Numerical Eigenvalue Problems and the SVD

An eigenvalue is a root of the characteristic polynomial, and for matrices of size $m\ge 5$ the Abel–Ruffini theorem forbids a closed-form root formula in radicals. Eigenvalues must therefore be iterated toward, not solved for. Two iterations do it: power iteration for a single eigenpair, the QR algorithm for the whole spectrum. The singular value decomposition follows from the same tools, again without squaring the condition number.

The characteristic polynomial versus Schur form

An eigenpair $(\lambda ,x)$ of $A\in {\mathbb{C}}^{m\times m}$ satisfies $Ax=\lambda x$ with $x\ne 0$; the eigenvalues form the spectrum $\sigma (A)$, the roots of the characteristic polynomial $\chi (\zeta )=\det (\zeta I-A)$. The obvious algorithm — form $\chi$, then find its roots — fails numerically, because polynomial rootfinding is severely ill-conditioned.

For example, take the diagonal matrix with eigenvalues $1,2,\dots ,22$, whose spectrum is read off instantly. Forming the characteristic polynomial and rooting it returns $15.4\pm 1.07i$ near the true $\lambda =15$: an absolute error of about $1$, where the well-conditioned eigenvalue itself permits error $O({\epsilon }_{\mathrm{mach}})$. The condition number of that root as a function of the polynomial coefficients is $\approx 6\times 10^{16}$. The eigenvalue problem and the rootfinding problem have very different conditioning, and going through the polynomial destroys the good conditioning of the former.¹

The stable target instead is the Schur decomposition, built from unitary similarities that preserve the spectrum.

Because unitary similarity leaves both the spectrum and the spectral norm unchanged, every algorithm here works by driving $A$ toward triangular form through a sequence of unitary changes of basis. This is the same reason unitary matrices dominated the stable solvers: they do not amplify perturbations.

Conditioning of eigenvalues

Backward error transfers directly from linear systems. For an approximate eigenpair $(\lambda ,x)$, the residual $r=Ax-\lambda x$ gives the backward error $\omega =\Vert r{\Vert }_2/(\Vert A{\Vert }_2\Vert x{\Vert }_2)$, and a computed eigenpair with $\omega =O({\epsilon }_{\mathrm{mach}})$ is backward stable. For a fixed approximate eigenvector, the eigenvalue minimizing that residual is the Rayleigh quotient

the least-squares fit of a scalar $\mu$ to $Ax\approx \mu x$. Two perturbation theorems govern the forward conditioning.²

• Eigenvalues of normal matrices are perfectly conditioned. By the Bauer–Fike theorem, if $A$ is normal and $\lambda \in \sigma (A+E)$, then $\mathrm{dist}(\lambda ,\sigma (A))\le \Vert E{\Vert }_2$: absolute condition number exactly $1$. A perturbation of size $\epsilon$ moves each eigenvalue by at most $\epsilon$.
• Eigenvectors are conditioned by the spectral gap. The Davis–Kahan theorem bounds the rotation of a simple eigenvector by $\Vert E{\Vert }_2$ divided by the distance from $\lambda$ to the rest of the spectrum. Two nearly equal eigenvalues make their eigenvectors ill-defined — in the limit, any direction in the degenerate eigenspace is an eigenvector.

An eigenvector is ill-conditioned exactly when its eigenvalue is close to another, and the same closeness slows every iteration below.

Power iteration

The simplest iteration repeatedly applies $A$ and renormalizes. Starting from $v_0$, for $k=1,2,\dots$:

Since $v_k$ is $A^k{v}_0$ normalized, the component along the largest eigenvalue grows fastest and dominates.

The eigenvector error decays geometrically at rate $\theta$, and the Rayleigh quotient converges twice as fast. Convergence is fast when the top two eigenvalues are well separated and slow when $\theta \approx 1$ — precisely the ill-conditioned eigenvector case. Backward stability is reached after roughly $15/|{\log }_{10}\theta |$ steps.³

The angle falls by $\theta =1/3$ per step and the eigenvalue error by ${\theta }^2=1/9$, matching the theorem's $O({\theta }^k)$ and $O({\theta }^{2k})$ rates.

Inverse iteration targets any eigenvalue, not just the dominant one, using the equivalence that $(\lambda ,x)$ is an eigenpair of $A$ iff $((\lambda -\mu {)}^{-1},x)$ is one of $(A-\mu I{)}^{-1}$. A shift $\mu$ near a chosen eigenvalue ${\lambda }_1$ makes $({\lambda }_1-\mu {)}^{-1}$ dominant, so power iteration on $(A-\mu I{)}^{-1}$ converges to $x_1$:

The matrix $A-\mu I$ is fixed, so a single factorization serves every step. A tempting objection is that a shift near the spectrum makes $A-\mu I$ ill-conditioned, so $w_k$ is computed inaccurately. The objection is misguided: only the direction of $w_k$ matters, and a backward stable solve returns the exact solution for a slightly perturbed $A+E$, whose eigenvector is what the iteration converges to. In fact, using a backward stable eigenvalue as the shift, one step of inverse iteration from a random start computes a backward stable eigenpair.³

The QR algorithm

Inverse iteration finds one eigenvector; the QR algorithm finds the whole Schur form by folding inverse iteration into a similarity transform. From $A_0=A$, with shifts ${\mu }_k$:

Each $A_{k+1}=Q_k^{\ast }{A}_k{Q}_k$ is unitarily similar to $A_k$, so the spectrum is preserved throughout. With a constant shift the recurrence reduces to inverse iteration keeping the eigenvector approximation pinned at the last basis vector, and the iterates converge to Schur form: the subdiagonal entries shrink to zero and the eigenvalues surface on the diagonal, deflated one at a time from the bottom.⁴

The shift determines the speed. The Rayleigh shift ${\mu }_k=(A_k{)}_{mm}$ gives locally quadratic convergence in general and locally cubic convergence for normal matrices, so the number of correct digits doubles or triples each step. It fails on real matrices with complex-conjugate eigenvalue pairs, which the Wilkinson shift — the eigenvalue of the trailing $2\times 2$ block nearest $(A_k{)}_{mm}$ — resolves. In practice, only two to four iterations per eigenvalue are needed.⁴

The correct digits roughly triple each step: the shift makes convergence cubic for a symmetric matrix, so three iterations drive the entry below machine precision and one eigenvalue deflates. Unshifted, the same entry would fall only linearly, by a constant factor $|{\lambda }_2/{\lambda }_1|$ per step.

Making it cubic

A naive QR step costs $O(m^3)$, and $O(m)$ eigenvalues would give $O(m^4)$ overall. The fix is to precondition $A$ into a shape a QR step preserves and can sweep cheaply. That shape is upper Hessenberg (zero below the first subdiagonal): the space of Hessenberg matrices is invariant under a QR step, and a QR step on a Hessenberg matrix costs only $O(m^2)$ using Givens rotations. So a Schur decomposition proceeds in two phases.

Individual eigenvectors are then recovered by inverse iteration with the computed eigenvalues as shifts. Computing only the eigenvalues is far cheaper than a full Schur decomposition, since the unitary transformations need not be accumulated.

Computing the SVD

The singular values of $A\in {\mathbb{K}}^{m\times n}$ are the square roots of the eigenvalues of the s.p.d. matrix $A^{\ast }A$,

which suggests forming $A^{\ast }A$ and running the symmetric eigenvalue algorithm on it. That route repeats the mistake of the normal equations: forming $A^{\ast }A$ squares the condition number, and the small singular values — usually the ones of interest — are computed with squared error. The stable method never forms $A^{\ast }A$. Instead it works on $A$ directly, in two phases mirroring the eigenvalue algorithm:

• Bidiagonalize. Apply Householder reflectors on the left and right to reduce $A$ to an upper bidiagonal matrix $B$, with the same singular values as $A$.⁵
• Diagonalize. Run an implicit-shift QR iteration on $B$ that operates as if on the tridiagonal $B^{\ast }B$ without ever forming it, driving $B$ to a diagonal of singular values.

The result is the singular value decomposition$A=U\Sigma V^{\ast }$, with the singular values accurate to full precision. Because the whole computation is a chain of unitary transformations, it inherits their backward stability: the computed factorization is the exact SVD of a matrix within $O({\epsilon }_{\mathrm{mach}})\Vert A\Vert$ of $A$.

Summary

Every algorithm here reduces to two ideas. Hard matrices are pushed toward triangular or diagonal form, and the pushing is done with unitary transformations because they preserve length and so do not amplify error. The $LU$ and Cholesky factorizations triangularize for solving systems; QR and Householder orthogonalize for least squares; the QR algorithm iterates orthogonalization into the Schur form; and the SVD does it from both sides at once. Underneath all of it sits the rule of thumb — forward error at most condition times backward error — which decides, for every one of these methods, exactly how much accuracy the problem allows and the algorithm delivers.