-1
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ 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?) $\endgroup$ – rchilton1980 Mar 17 at 16:09
  • $\begingroup$ @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. $\endgroup$ – AsianInvasi0n Mar 17 at 17:47
  • $\begingroup$ 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? $\endgroup$ – Federico Poloni Mar 17 at 20:48
  • $\begingroup$ @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 $\endgroup$ – AsianInvasi0n Mar 18 at 12:30
  • $\begingroup$ @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). $\endgroup$ – Federico Poloni Mar 20 at 12:52
2
$\begingroup$

Yes, those Givens rotations do not seem correctly implemented to me. Since you are doing this as a learning project: learn from the masters http://www.netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga54d516b6e0497df179c9831f2554b1f9.html.

Also make sure that the givens rotations you are applying from the left are the same as the ones you are applying from the right. I haven't checked it myself, but it's not obvious to me that this is the case in your current implementation. The easiest way would be to store the rotations in some array or implement the more popular bulge chasing version of the QR algorithm instead.

Lastly, there are many other details in an implementation of the QR algorithm that need to be adressed: shifts to speed up convergence, deflations, over/underflow. Again, I recommend you learn from the masters: https://www.netlib.org/lapack/explore-html/d3/d03/dlahqr_8f_source.html#l00209

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the tips. I have checked the orthogonality of Q (i.e. the error is something like (+-)5e-16). Unfortunately, adapting the Givens code from LAPACK didn't work so I decided to move onto the bulge chasing version. I also decided to move on because implementation-wise, I believe it more important to be able to implement the "most relevant" options (I have no idea why this slipped my mind earlier). Currently I am working on understanding the components of the bulge chasing version of the QR algorithm and will get back to you as to whether or not it works. $\endgroup$ – AsianInvasi0n Mar 18 at 17:02
  • $\begingroup$ The bulge chasing algorithm is based on Givens rotations too, you need to get those functions right. There is a reference implementation in LAPACK, see slartg and dlartg. $\endgroup$ – rchilton1980 Mar 23 at 16:35
  • $\begingroup$ Well technically LAPACK uses small householder reflectors, but yes. The main difference I was pointing to is the order in which the transformations are applied. That way it's easier to make sure the left and right transformations are the same. $\endgroup$ – Thijs Steel Mar 24 at 7:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.