# How to implement the Hessenberg QR Algorithm?

For context, I'm creating a linear algebra library from scratch for learning purposes in C. Right now I'm working on calculating eigenvalues but my implementation of the QR Algorithm is diverging. Below is my understanding of the QR Algorithm, and after the first block of code is (what I believe to be) the meat of the problem: Givens rotations.

Note that the code given is python-like pseudocode, not C (though the code given doesn't look much different than the source).

I'm having trouble understanding what's going wrong here and why. Details are below:

My implementation diverges right off the bat, and is definitely noticeable after 3 iterations. My current understanding of the algorithm (using a Hessenberg) is this:

Given $$n \times n$$ Hessenberg matrix $$\mathbf H$$, I want to calculate its eigenvalues with the QR Algorithm. To do so, I must apply $$n - 1$$ transposed Givens rotations from the left (i.e. $$\mathbf G^T_{n-1}\mathbf G^T_{n-2}...\mathbf G^T_1\mathbf H$$). The givens rotations are equivalent to $$\mathbf Q^T$$. As such, we can multiply in reverse order to fulfill the QR algorithm (i.e. $$\mathbf {\overline H} = \mathbf H\mathbf G_1\mathbf G_2...\mathbf G_{n-1}$$). This eventually converges into an upper triangular, on whose diagonal the corresponding eigenvalues are to be found.

In code:

for i = 0 until convergence:
for k = 0 to n-1:       // rotation to make H upper triangular
givensLeft(H, k)
for k = 0 to n-1:       // back to Hessenberg
givensRight(H, k)


The Givens rotations are my primary worry as I'm not entirely sure that they're working correctly. I understand that only rows/columns $$k$$ and $$k+1$$ are affected. As such, I can perform a Givens rotation from the left by multiplying the submatrix $$\mathbf H_{k:k+1, \ k:n}$$ by the transposed Givens matrix $$\mathbf G^T$$: $$\begin{bmatrix} c_k &-s_k \\ s_k & c_k \end{bmatrix}$$

To perform a Givens rotation from the right (in the QR algorithm this would be retruning the Hessenberg back to its form from the upper triangle caused by the left Givens rotation), I would multiply submatrix $$\mathbf H_{1:k+1, \ k:k+1}$$ by the (not transposed) Givens matrix $$\mathbf G$$: $$\begin{bmatrix} c_k & s_k \\ -s_k & c_k \end{bmatrix}$$

In code:

givensLeft(H, k):
a = H[k][k]
b = H[k+1][k]
r = hypot(a, b)
c = a/r
s = -1*b/r

// these could be more efficient, I know
for i = k+1 to H.columns:
a = H[k][i]
b = H[k+1][i]
H[k][i] = c*a - s*b

for i = k+1 to H.columns:
a = H[k][i]
b = H[k+1][i]
H[k][i] = s*a + c*b

givensRight(H, k):
a = H[0][k]
b = H[0][k+1]
r = hypot(a, b)
c = a/r
s = -1*b/r

for i = 0 to k+1:
a = H[i][k]
b = H[i][k+1]
H[i][k] = c*a - s*b

for i = 0 to k+1:
a = H[i][k]
b = H[i][k+1]
H[i][k+1] = s*a + c*b


The only problem I can think of is that of the Givens rotations, as they seem to make up the whole QR algorithm, but I have no idea what the problem would be. I've looked at several things online outlining the QR algorithm with a Hessenberg (e.g. this and this) and this implementation is written directly from the former (using Algorithms 4.1 and 4.2). I have not been able to find anything on here or elsewhere with people having a similar problem.

Is my problem conceptual? or is it just faulty programming? Many thanks in advance.

• Hard to tell what language you're using, but I would suggest you implement this first in a language that puts a rich linear algebra library right at your fingertips. (matlab/octave or python+numpy spring to mind). Then you'll be able to check your work as you go. In particular, unitary transformations mean that H should have the same eigenvalues after every iteration, you can check this easily in matlab using eig(). You could also test that your givens matrix does the right thing on the 2x2 case (for example, does it preserve eigenvalues and truly yield a zero where you expect it to?) – rchilton1980 Mar 17 at 16:09
• @rchilton1980 It's in C, sorry I should have been more clear. As for using another library, I'm trying to avoid doing so because this is primarily a programming and math learning exercise. With the Givens functions, I have stepped through them and they are zeroing and filling the right spots, so I really don't know what's going on. If you're curious as to the source (which the above doesn't look much different at all than the source), here it is. – AsianInvasi0n Mar 17 at 17:47
• What is the question here? You mention a "problem", but you never describe it. You write that "The Givens rotations are my primary worry as I'm not entirely sure that they're working correctly". Have you tested your implementation? Does it pass tests? – Federico Poloni Mar 17 at 20:48
• @Frederico Poloni The main problem is that I'm unsure how to verify that my work is working correctly aside from the algorithm converging i.e I lack understanding of all of this, am having trouble finding and understanding resources, and am looking for guidance on where to proceed. With the Givens rotations, I have been unable to find enough to allow me to verify that it is working beyond the change in shape. The change of shape is working, and that is what led to my confusion – AsianInvasi0n Mar 18 at 12:30
• @AsianInvasi0n These Givens transformations construct an orthogonal similarity of your matrix $A$. A way to test is accumulating the product of all rotations $Q$, and check that at each step $QHQ^T=A$ (and $Q^TQ=I$, for good measure). (Also, please use autocompletion or copypaste to avoid mistyping names in at-notifications). – Federico Poloni Mar 20 at 12:52

An additional hint is that you can also compute the matrix $$Q$$. You can then verify whether it is orthogonal ($$Q^TQ-I$$ should be small) and you can also check the backwards error ($$Q^THQ - \overline{H}$$ should be small). If any of these become large you can even check this line by line to find exactly where the error occurs.