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.
