# Why SVD on $X$ is preferred to eigendecomposition of $XX^\top$ in PCA

In this post J.M. has mentioned that ...

In fact, using the SVD to perform PCA makes much better sense numerically than forming the covariance matrix to begin with, since the formation of $XX^\top$ can cause loss of precision. This is detailed in books on numerical linear algebra, but I'll leave you with an example of a matrix that can be stable SVD'd, but forming $XX^\top$ can be disastrous ...

Could you elaborate on why the calculation of $XX^\top$ is disastrous for the matrix given in that post? I calculated $XX^\top$ for the numbers it specifies. I got a $3\times3$ matrix with all $1$'s.

• This post apparently refers to this question. – Andrew Uzzell Apr 12 '13 at 13:05
• I changed your title to reflect the question more precisely. Please feel free to roll back if you think the change is inappropriate. – user1551 Apr 13 '13 at 4:50

Alright, to further explain Läuchli's example in detail:

The precise problem with forming the cross-product matrix in this case is that you have effectively made a singular matrix out of your starting matrix.

In exact arithmetic, the cross-product matrix of the Läuchli matrix should have looked like this:

$$\begin{pmatrix}1+\varepsilon^2&1&1\\1&1+\varepsilon^2&1\\1&1&1+\varepsilon^2\end{pmatrix}$$

Its exact eigenvalues are $(3+\varepsilon^2,\varepsilon^2,\varepsilon^2)$. The trouble with inexact arithmetic is that if $\varepsilon^2$ is tiny, then $1+\varepsilon^2$ is the same as $1$, and you now have the spectrum $(3,0,0)$. You do think there's a difference between $\bf \varepsilon^2$ and $\bf 0$, don't you?

The good thing, then, about direct SVD algorithms is that it is able to compute these tiny singular values accurately, much better than taking the eigensystem of the cross-product matrix. Here's a Mathematica comparison:

lauchli = With[{ε = 1*^-20},
N[{{1, 1, 1}, {ε, 0, 0}, {0, ε, 0}, {0, 0, ε}}]];

SingularValueList[lauchli, Tolerance -> 0]
{1.73205, 1.*10^-20, 1.*10^-20}

Sqrt[Eigenvalues[Transpose[lauchli].lauchli]]
{1.73205, 0. + 1.82501*10^-8 I, 0.}

(A similar demonstration could of course be done in MATLAB; luckily this example is built-in there, as gallery('lauchli', n))

Now, we even have a spurious (but tiny) imaginary component in the second answer; neither of the last two results matches the first.

That's why the SVD is vastly preferred.

• +1. I used (and referenced) your example in my answer to the almost identical question on CV: stats.stackexchange.com/questions/79043. You might be interested in adding something to that discussion. – amoeba Jun 30 '16 at 18:09