Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.

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)?

share|improve this question
2  
First: Yes. Can you think of an example? (Hint: Try to find an example where every orthogonal matrix is a set of eigenvectors). Second (partial): Eigenvectors are only guaranteed to be orthogonal when the underlying matrix is symmetric. If you think about that, that answers a lot of your question right away. For diagonalizable matrices an eigendecomposition takes the form $X = P D P^{-1}$. – cardinal Jan 5 '12 at 21:18
1  
Is this the right place for this question (asking sincerely, not rhetorically)? Certainly linear algebra is important in statistics, but I don't see any statistical content here per se. – bnaul Jan 5 '12 at 21:33
Also, I will defer to @cardinal in terms of how much help is appropriate to give here, but I think it would be helpful to explain/link to the correct way to find eigenvalues via QR decomposition, which IMO is not at all obvious and would be hard to discover even with good hints. – bnaul Jan 5 '12 at 21:40
I believe this question might be of interest for people using multivariate statistics where matrix decomposition are extensively used, although it focus on linear algebra and not statistics per se (which probably explains the 2 votes to close). Do you have a possible statistical application in mind? – chl Jan 5 '12 at 22:50
1  
@Elvis, I think some elaboration on your comment in the other question would make a fine answer to this one. Probably no need to go into shifts and other ways of speeding up the algorithm you allude to, though the basic one has somewhat poor convergence properties. – cardinal Jan 6 '12 at 13:39
show 2 more comments

1 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.

share|improve this answer
Thanks much. I've implemented it, but it seems it takes a long time for this algorithm to converge. Could you please also briefly explain the SVD algorithm (similarly to the above)? I read about it, but did not find an explanation as simple and to-the-point as yours. – Bliss Jan 8 '12 at 21:03
@YanRaf 1) What convergence criterion do you use? You shouldn’t wait until non-diagonal coefficients are null, but stop for example as soon as they are < 0.001 in absolute value... this should be precise enough. 2), the SVD. This is much more complicated, there is a part to construct a tridiagonal matrix that I don’t remember at all. The brute force version is to compute the product $A \times A'$ or $A' \times A$ to have the smallest possible square symmetric matrix, and to compute its PC that are the singular vectors of A. – Elvis Jan 12 '12 at 21:59
Thanks. 1) I'm using this exact convergence criterion, i.e. non-diagonal coefficients (absolute values) should be epsilon distant from 0, at most. – Bliss Jan 15 '12 at 15:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.