Trying to perform subspace iteration by applying a matrix $A\in\mathbb{R}^3$ repeately to a matrix $S$, but its not working?

109 Views Asked by At

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
1

There are 1 best solutions below

0
On BEST ANSWER

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,

[S,R]=qr(A*S). 

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.

for J=1:20
    [S,R] = qr(A*S); S
end

.

 S  =
  - 0.4856429  - 0.8741573    0.  
  - 0.8741573    0.4856429    0.  
    0.           0.           1.  
 S  =
  - 0.6794080  - 0.7337607    0.  
  - 0.7337607    0.6794080    0.  
    0.           0.           1.  
...

 S  =
  - 0.9999996  - 0.0008463    0.  
  - 0.0008463    0.9999996    0.  
    0.           0.           1.  
 S  =     
  - 0.9999999  - 0.0005078    0.  
  - 0.0005078    0.9999999    0.  
    0.           0.           1.  
 S  =
  - 1.0000000  - 0.0003047    0.  
  - 0.0003047    1.0000000    0.  
    0.           0.           1.  
...