How to get eigenvectors using QR algorithm?

Question

From everything I've heard, this matlab code ought to spit out a matrix where each row is the same. So why doesn't it?

A = [ 5  2  0  0;
      3  9  4  0;
      0  9  5 -2;
      0  0 -3  4 ];
B=A;
QQQ=eye(4);

%QR algorithm
for i=1:100
    [Q,R] = qr(B);
    B=R*Q;
    QQQ = QQQ*Q;
end

(A*QQQ)./QQQ %should have constant rows, but doesn't

Answer 1

First, this is not true

(A*QQQ)./QQQ %should have constant rows, but doesn't

You basically state that $$ A Q = Q \Lambda, $$ where $Q$ is an orthogonal matrix and $\Lambda = \operatorname{diag} (\lambda_1, \dots, \lambda_n)$. This means that $A$ is ortogonaly similar to a diagonal matrix, but it is true only for normal matrices which is not the case.

After each iteration the QR algorithm the following relation holds: $$ A Q_k = Q_k B_k $$ But $B_k$ converge to an upper triangular matrix, not a diagonal one, as you might expect. Thus QR algorithm computes the Schur decomposition of the matrix, not its eigendecomposition.

Consider now the Schur decomposition of the original matrix: $$ A Q = Q R $$ It is not hard to obtain eigenvectors when the Schur decomposition is known. We need to find eigenvectors of $R$ (and the eigenvalues are already known - they are on the main diagonal of $R$).

Let $\mathbf v_i$ be the $i$-th eigenvector of $R$: $$ R \mathbf v_i = r_{ii} \mathbf v_i \implies (R - r_{ii} I) \mathbf v_i = \mathbf 0. $$ The corresponding $\mathbf v_i$ may be sought in form $$ \mathbf v_i = \begin{pmatrix} \ast\\ \ast\\ \vdots\\ \ast\\ 1\\ 0\\ \vdots\\ 0 \end{pmatrix} $$ with $1$ on the $i$-th row. Plugging this to the equation we obtain a triangular system for $\ast$ unknown values: $$ R(1:i-1, 1:i-1) v_i(1:i-1) + R(1:i-1, i) - r_{ii} v_i(1:i-1) = 0 $$ Here I've used the Matlab notation for denoting submatrices. The solution of the system is $$ v_i(1:i-1) = \left[ r_{ii} I - R(1:i-1, 1:i-1) \right]^{-1} R(1:i-1, i) $$

After we have found the eigendecomposition of $R$ with $V = [\mathbf v_1, \dots, \mathbf v_n]$ $$ R V = V \operatorname{diag}(R) $$ we can plug it into the Schur decomposition of $A$: $$ A Q V = Q R V = Q V \operatorname{diag}(R). $$ It is now clear that $QV$ form the eigenvectors of $A$.

Here's an updated version of the code: https://gist.github.com/uranix/2b4bb821a0e3ffc4531bec547ea67727