Finding matrix eigenvectors using QR decomposition

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

• 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
• 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
• @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

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

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

• 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
• Hi, I think I spotted an error in your code. I tried it and this is the results I got : gist.github.com/gbersac/4335bc17ca38fe07f7b3 The eigen vectors are ok, but the eigen values are wrong. Could you fix that please. – Moebius Apr 15 '15 at 12:15
• The error is in your code. Try with A <- Origin instead. – Elvis Apr 15 '15 at 12:25

Sorry it's a bit late, but I think there is still room for basic answers.

First, a general linear algebra question: Can a matrix have more than one set of (unit size) eigenvectors?

If the eigenvalues are distinct, the answer is "No". If there are duplicate eigenvalues, then for these eigenvalues, the eigenvectors are not distinct (but any eigenvectors corresponding to unique eigenvalues are still distinct).

For example, if you are looking at a matrix with iid Gaussian entries, then barring some floating point fluke, it will not have repeat eigenvalues, and hence the eigenvectors will be uniquely determined.

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

The QR does not give this to you, but the rank-revealing QR (RRQR) does, and there are known bounds on the error. You should not implement either the QR or the RRQR yourself since there is excellent code and this is a major topic of research in numerical linear algebra. (The big issues: stability and taking advantage of memory and communication to make it efficient)

For a reference on the RRQR giving bounds on the eigenvalues, try Some applications of the rank revealing QR factorization (1992), by T F Chan and P C Hansen.

Also, be careful with the distinction of the QR Factorization and the QR Algorithm. The QR Algorithm, which the other answer shows, uses QR factorizations at every step, hence the name, but otherwise they are different algorithms. You asked about the QR Factorization/Decomposition, and got an answer with the QR Algorithm.

On wikipedia, check out "QR_algorithm" compared to "QR_decomposition"

• Re the first answer: you should be more explicit about what you mean by two sets of eigenvectors to be the "same". After all, the $1$ by $1$ matrix $(1)$ has distinct eigenvalues $\{1\}$ for which the unit vectors $(1)$ and $(-1)$ are both eigenvectors, but few people would consider the sets $\{(1)\}$ and $\{(-1)\}$ to be the same. – whuber Jun 12 '13 at 11:51
• True. Any non-zero multiple of an eigenvector is still an eigenvector (and even with the SVD, there is still a +/- issue). So what I mean by "distinct" is that two vectors are distinct if they are linearly independent. Basically, every eigenvalue corresponds to an eigenspace, and the dimension of that eigenspace matches the multiplicity of the eigenvalue. So if the multiplicity is larger than 1, then there are infinitely many non-trivial choices of a basis (i.e. eigenvectors) for that space. – Stephen Jun 13 '13 at 9:00