How to get eigenvectors using QR algorithm? 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

 A: 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
