I'm trying to code up a matrix library (purely as a learning exercise). This question is about the math I'm trying to understand in order to implement it. I just want to make sure I have a firm grasp of the underlying algorithm.

For reference, I'm referring mainly to this book for the algorithms and mathematical understanding: http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter3.pdf

First, basic QR algorithm: Given an N x N non-symmetric real A, let QR(A) be the QR Decomposition into an orthogonal matrix Q and an upper triangular matrix R.

Then, iterate QR(Ak-1), Ak = Rk*Qk until the last element in the first sub diagonal is close (enough) to zero. Take the last element (in the main diagonal) as an eigenvalue and deflate the matrix by removing its last row and column and begin the process again, until the entire matrix has been consumed.

To handle complex eigenvalues, I've assumed that if this element does not converge to zero before a specified iteration_max, then the trailing 2x2 sub matrix contains two eigenvalues (which still may be real eigenvalues, I think, so I check for that). I apply a straight quadratic equation to that sub-matrix's characteristic polynomial to get its eigenvalues, and deflate the A matrix by two rows and columns, and iterate.

Everything before this point seems to work when I implement it in code.

Now, the problem is that the convergence is much too slow. I have to set the maximum iterations severely high in order to handle some of my test cases. So I went and read up on advanced QR tricks.

One is to reduce the input matrix to upper Hessenberg form, which I believe I've done correctly through the transformation H = Hessenberg(A). However, for one of my test cases, this causes the required iterations to actually increase. That makes me scratch my head.

More importantly, I tried to implement shifting, and this broke my algorithm. So I need to make sure I understand that.

Let the shift S be defined as the last element on the main diagonal. Let I be defined as an appropriately sized identity matrix.

Then, I do QR(Hk - I * Sk). And Hk+1 = R * Q + I * Sk.

This causes my matrix to rapidly converge. Unfortunately, it also effectively zeroes the last element in the main diagonal of H in my test case, causing it to incorrectly report an eigenvalue near zero.

What am I doing wrong?

Edit: (please don't laugh at my little 4 x 4 matrix) Test matrix:

double[][] Matrix = {
        {-2.476, -2.814, 4.29, -3.649},
        {2.839, -2.859, 1.623, -2.926},
        {-0.392, -3.206, -0.401, -2.174},
        {2.241, -4.435, -3.963, 4.102}};

Edit: link to repository with code. Relevant function is large or I would embed it. Located in SquareMatrix.eigenvalues() https://github.com/rwthompsonii/matrix-java

  • $\begingroup$ Can you post your test matrix and the $H$ in the iteration before it is zeroed? Also, the code would be helpful. Are you always choosing the last element of the main diagonal, or are you choosing the diagonal element above the unconverged part? $\endgroup$ – Victor Liu Mar 19 '15 at 5:38
  • $\begingroup$ I always take the very last element in the main diagonal if the last element in the main sub-diagonal has converged to ~0. $\endgroup$ – Ron Thompson Mar 19 '15 at 5:42
  • $\begingroup$ Apologies if this is what you're doing, but I'm not sure I understand what you're saying. If the subdiagonal element $h_{j+1,j}$ has converged to zero, then the next shift should be $h_{j,j}$. Also, I stand corrected... the code does not help :) Java is so verbose... $\endgroup$ – Victor Liu Mar 19 '15 at 5:45
  • $\begingroup$ In 1based array notafion if matrix[n][n-1] ~ 0 then matrix [n][n] is selected as an eigenvalue else the next shift is = matrix[n][n]. Is that right? $\endgroup$ – Ron Thompson Mar 19 '15 at 5:50
  • $\begingroup$ Yeah I know. I'm learning it or I would do it in C. $\endgroup$ – Ron Thompson Mar 19 '15 at 5:51

The problem is that a real upper Hessenberg matrix can have complex eigenvalues, which your code seemingly allows for. However, if you always choose your shift to be the last diagonal element above the converged part, it will always be real, since the QR decomposition of a real shifted upper Hessenberg matrix is still real. If you look at section 3.5,

Another problem occurs if real Hessenberg matrices have complex eigenvalues. We know that for reasonable convergence rates the shifts must be complex.

Hence the need for the double shift algorithm.

  • $\begingroup$ Oh! The single shift method ONLY works for real eigenvalues? I didn't get that from that paragraph, I just thought it meant it wouldn't converge, not that it would converge incorrectly. $\endgroup$ – Ron Thompson Mar 19 '15 at 6:08
  • $\begingroup$ Right now I'm choosing the very bottom right element as the shift if the element DIRECTLY to its left has not converged, is that properly done single shift? $\endgroup$ – Ron Thompson Mar 19 '15 at 6:09
  • $\begingroup$ @Ron Thompson: I believe so. $\endgroup$ – Victor Liu Mar 19 '15 at 16:03
  • $\begingroup$ Ok then. Thanks a bunch, I'll accept your answer and go work on double QR shift then. $\endgroup$ – Ron Thompson Mar 19 '15 at 17:48

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.