Working with Givens rotations

Question

If we consider a matrix R of size pxp. If we want to multiply A'RA where A is equal to (I+Givens rotation). Here I is an identity matrix and ' denotes the transpose operator.

We know that a Givens rotation is a sparse matrix written as:

To perform the multiplication A'RA in matlab, we can do this fast implementation:

%Fast implementation
  ci = R(:,ik)*(cos(theta))+R(:,jk)*(sin(theta)); % R*A
  cj = R(:,jk)*(cos(theta)) - R(:,ik)*(sin(theta));
  R(:,ik) = ci;
  R(:,jk) = cj;

  ri = R(ik,:)*(cos(theta))+R(jk,:)*(sin(theta)); % A'*R*A
  rj = R(jk,:)*(cos(theta)) - R(ik,:)*(sin(theta));
  R(ik,:) = ri;  
  R(jk,:) = rj;

But I didn't understand how they wrote this Matlab code. In other terms, I am not understanding how this Matlab code apply the multiplication A'RA. Kindly, can someone help me to understand this code?

Answer 1

One possible source of confusion is that either the signs in the Givens rotation matrix, or the side on which we need to transpose, is wrong in your example. I'll assume the latter: I'll use the same A matrix as you defined, but transform with A*R*A' (changing the A to transpose is equivalent to taking the rotation angle with opposite sign).

The algorithm is relatively straightforward. For starters, as the comments in the code suggest, the transformation is performed in two steps:

Rnew = A * R * A' = A * (R * A')

First, we compute R*A'. For this, imagine the transformation matrix A = I + M with the Givens rotation matrix M. The formula which you showed basically says "Take a unit matrix, except for 2 specified dimensions in which you rotate by a given angle". Here's how the full A matrix looks like for a small problem (6d matrix, ik=2, jk=4, both in full and sparse form):

You can see that except for the (ik,jk) 2d subspace, this matrix is a unit matrix, leaving every other dimension intact. So the action of R*A' will result in R for every dimension except for columns ik and jk.

In these two columns the result of R*A' is the linear combination of R(:,ik) and R(:,jk) with these trigonometric coefficients:

[R*A'](:,ik) =  R(:,ik)*cos(theta) + R(:,jk)*sin(theta)
[R*A'](:,jk) = -R(:,ik)*sin(theta) + R(:,jk)*cos(theta)

while the rest of the columns are left unchanged. If you look at the code you cited: this is exactly what it's doing. This is, by definition, what R*A' means with the A matrix shown above. All of this is the implication of that the A matrix is a unit matrix except for a 2d subspace.

The next step is then quite similar: using this new R*A' matrix we multiply from the left with A. Again, the effect along most of the dimensions (rows, this time) will be identity, but in rows ik and jk we again get a linear combination:

[A*[R*A']](ik,:) =  cos(theta)*[R*A'](ik,:) + sin(theta)*[R*A'](jk,:)
[A*[R*A']](jk,:) = -sin(theta)*[R*A'](ik,:) + cos(theta)*[R*A'](jk,:)

By noting that the code overwrites the R matrix with R*A' after the first step, it's again clear that the same is performed in the "fast implementation" code.

Disclaimer: A' is the adjoint (conjugate transpose) in matlab, so you should use A.' to refer to the transpose. For complex matrices there's a huge difference, and people often forget to use the proper transpose when eventually encountering complex matrices.