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.
