Solved – Finding matrix eigenvectors using QR decomposition

linear algebramatrix decompositionsvd

First, a general linear algebra question: Can a matrix have more than one set of (unit size) eigenvectors? From a different angle: Is it possible that different decomposition methods/algorithms (QR, NIPALS, SVD, Householder etc.) give different sets of eigenvectors for the same matrix?

Second, regarding QR decomposition: Are the columns of the Q matrix the eigenvectors? How can their eigenvalues be easily found (post the QR decomposition)?

Best Answer

The (basic) algorithm with QR decomposition is as follows.

  • Let $X$ by a symmetric matrix.

  • Let $X_1 = X$, and iterate the following:

  • Given $X_k$, write a QR decomposition $X_k = Q_k R_k$, and let $X_{k+1} = R_k Q_k$;

  • The matrices sequence $X_n$ converges to some diagonal matrix $D$ with the eigenvalues on the diagonal; you retrieve the corresponding eigenvectors as the columns of $\prod_i Q_i$.

Here is an example code in R.

# some symmetric matrix
A <- matrix( sample(1:30,16), ncol=4)
A <- A + t(A);

# initialize
X <- A;
pQ <- diag(1, dim(A)[1]);

# iterate 
for(i in 1:30)
{
  d <- qr(X);
  Q <- qr.Q(d);
  pQ <- pQ %*% Q;
  X <- qr.R(d) %*% Q;
}

Now we have a look on the result

> A
     [,1] [,2] [,3] [,4]
[1,]   52   30   49   28
[2,]   30   50    8   44
[3,]   49    8   46   16
[4,]   28   44   16   22

The matrix X contains the eigenvalues on the diagonal:

> round(X,5)
         [,1]    [,2]      [,3]     [,4]
[1,] 132.6279  0.0000   0.00000  0.00000
[2,]   0.0000 52.4423   0.00000  0.00000
[3,]   0.0000  0.0000 -11.54113  0.00000
[4,]   0.0000  0.0000   0.00000 -3.52904

And the product of all Q contains the eigenvectors:

> round(pQ,5)
        [,1]     [,2]     [,3]     [,4]
[1,] 0.60946 -0.29992 -0.09988 -0.72707
[2,] 0.48785  0.65200  0.57725  0.06069
[3,] 0.46658 -0.60196  0.22156  0.60898
[4,] 0.41577  0.35013 -0.77956  0.31117

We can compare to the result of eigen(A) :

> eigen(A)
$values
[1] 132.627875  52.442300  -3.529045 -11.541131

$vectors
           [,1]       [,2]        [,3]        [,4]
[1,] -0.6094595 -0.2999194  0.72707077  0.09987744
[2,] -0.4878528  0.6519967 -0.06068999 -0.57724915
[3,] -0.4665778 -0.6019623 -0.60897966 -0.22156327
[4,] -0.4157690  0.3501285 -0.31117293  0.77956246

Of there is room for lots of improvements, but basically here it is. I once read lots of papers on the subject but my memory is leacking :(

Note that, as your problem is to perform PCA, you will find easily many PCA programs on the internet, you may prefer to do than rather than program it yourself.