I am trying to perform subspace iteration in MatLab but I am not getting the results I expected. Let the matrices $A$ and $S$ be given by $$ A = \begin{pmatrix} 5 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 2 \end{pmatrix}, \quad S = \begin{pmatrix} 3 & 4 & 6 \\ 9 & 5 & 2 \\ 0 & 0 & 0 \end{pmatrix}. $$
$A$ clearly has the eigenvalues $5, 3$, and $2$, and $S$ defines a subspace of dimension $k = 2$. I then assumed that repeated application of $A$ to $S$, followed by normalization, should 'push' the subspace represented by $S$ to the subspace given by the 2 dominant eigenvectors (the ones associated with eigenvalues $5$ and $3$) of $A$. I would have expected to end up with a matrix very close to the following matrix $M$, after performing many iterations: $$ M = A^k S = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0.6 & 0 \\ 0 & 0 & 0 \end{pmatrix}, $$ as $k$ becomes large.
However, when I perform this iterative process in MatLab I end up with a one-dimensional matrix $$ \hat{M} \approx \begin{pmatrix} 0.38 & 0.51 & 0.76 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$
I thought iteration with the matrix $A$ should take matrix $S$ to the $2$ dominant eigenvectors of $A$. So why hasn't it?
Edit
Here is my code:
A = [ 5, 0, 0 ; 0, 3, 0; 0, 0, 2 ];
S = [ 3, 4, 6; 9, 5, 2;, 0, 0, 0];
jMax = 20;
for j=1:jMax
S = A*S/norm(A*S);
disp(S);
end
Normalization for dimensions greater than $1$ requires to find a basis with some fixed properties like orthonormality for the subspace in each step. This can be achieved by the QR decomposition,
Then $$ S_kR_kR_{k−1}…R_1=A^kS_0, $$ so $S_k$ is an orthonormal basis of the space spanned by the rows of $M_k=A^kS_0$.
Executed (here in scilab) it results in the expected relatively slow linear convergence towards a diagonal matrix with a factor of $0.6$ or $1$ digit per $4.5$ steps in the off-diagonal elements.
.