---
title: Numerical Eigenvalue Problems and the SVD
module: Numerical Linear Algebra
moduleNumber: 8
lessonNumber: 6
order: 806
summary: >
  Eigenvalues cannot be found by a formula for large matrices, so they are found
  by iteration. Power and inverse iteration converge to one eigenvector at a rate
  set by the eigenvalue gap; the QR algorithm sweeps a matrix to Schur form and,
  with a good shift and a Hessenberg reduction, computes the whole spectrum in
  cubic time. Singular values follow from the same machinery applied without ever
  forming AᵀA.
topics: [Numerical Linear Algebra]
sources:
  - book: Bornemann
    ref: "Part V — Eigenvalue Problems; §18 Basic Concepts, §19 Perturbation Theory"
  - book: Bornemann
    ref: "§20 Power Iteration, §21 QR Algorithm"
draft: false
---

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(\varepsilon_{\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.[^bor-basic]

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

> **Definition (Schur decomposition).** Every $A \in \mathbb{C}^{m\times m}$ has a
> factorization $Q^\ast A Q = T$ with $Q$ unitary and $T$ upper triangular; the
> eigenvalues of $A$ appear on the diagonal of $T$. If $A$ is **normal** ($A^\ast A
> = A A^\ast$), then $T$ is diagonal, so $A = Q D Q^\ast$ and the columns of $Q$
> form an orthonormal basis of eigenvectors.[^bor-basic]

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](/linear-algebra/numerical-linear-algebra/qr-and-numerical-least-squares):
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 = \lVert r\rVert_2 / (\lVert A\rVert_2\lVert x\rVert_2)$, and a computed
eigenpair with $\omega = O(\varepsilon_{\mathrm{mach}})$ is backward stable. For a
fixed approximate eigenvector, the eigenvalue minimizing that residual is the
**Rayleigh quotient**

$$
\lambda = \frac{x^\ast A x}{x^\ast x},
$$

the least-squares fit of a scalar $\mu$ to $Ax \approx \mu x$. Two perturbation
theorems govern the forward conditioning.[^bor-pert]

- **Eigenvalues of normal matrices are perfectly conditioned.** By the Bauer–Fike
  theorem, if $A$ is normal and $\lambda \in \sigma(A + E)$, then
  $\operatorname{dist}(\lambda, \sigma(A)) \le \lVert E\rVert_2$: absolute
  condition number exactly $1$. A perturbation of size $\varepsilon$ moves each
  eigenvalue by at most $\varepsilon$.
- **Eigenvectors are conditioned by the spectral gap.** The Davis–Kahan theorem
  bounds the rotation of a simple eigenvector by $\lVert E\rVert_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$:

$$
w_k = A v_{k-1},
\qquad
\mu_k = v_{k-1}^\ast A v_{k-1},
\qquad
v_k = \frac{w_k}{\lVert w_k\rVert_2}.
$$

Since $v_k$ is $A^k v_0$ normalized, the component along the largest eigenvalue
grows fastest and dominates.

> **Theorem (Power iteration).** Let $A$ be normal with a **dominant** eigenvalue,
> $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert \ge \cdots$, and set $\theta =
> \lvert\lambda_2/\lambda_1\rvert < 1$. If $v_0$ has a nonzero component along the
> dominant eigenvector $x_1$, then
>
> $$
> \sin\angle(v_k, x_1) = O(\theta^k),
> \qquad
> \lvert\mu_k - \lambda_1\rvert = O(\theta^{2k}).
> $$

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/\lvert\log_{10}
\theta\rvert$ steps.[^bor-power]

$$
% caption: Power iteration rotates each iterate $v_k$ toward the dominant
% eigenvector $x_1$; the angle to $x_1$ shrinks by a factor $\theta =
% \lvert\lambda_2/\lambda_1\rvert$ each step.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\coordinate (O) at (0,0);
% dominant eigenvector direction
\draw[->, acc, thick] (O) -- (3.4,0) node[anchor=west] {$x_1$};
% successive iterates closing in
\draw[->, thick] (O) -- (60:2.6) node[anchor=south] {$v_0$};
\draw[->, thick] (O) -- (36:2.6) node[anchor=south] {$v_1$};
\draw[->, thick] (O) -- (19:2.6) node[anchor=south] {$v_2$};
\draw[->, thick] (O) -- (8:2.6) node[anchor=south west] {$v_3$};
% angle arcs
\draw[black] (1.2,0) arc (0:60:1.2);
\node[black] at (40:1.6) {angle};
\end{tikzpicture}
$$

> **Worked example.** Run power iteration on
>
> $$
> A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix},
> $$
>
> with eigenvalues $\lambda_1 = 3$ (eigenvector $x_1 = (1,1)^\top/\sqrt2$) and
> $\lambda_2 = 1$, so $\theta = \lvert\lambda_2/\lambda_1\rvert = 1/3$. Start from
> $v_0 = (1,0)^\top$. The angle to $x_1$ shrinks by a factor near $\theta$ each
> step, and the Rayleigh quotient $\mu_k$ approaches $3$ twice as fast.

| $k$ | $\mu_k$ | $\lvert 3 - \mu_k\rvert$ | $\angle(v_k, x_1)$ |
| --- | --- | --- | --- |
| 0 | $2.0000$ | $1.0\times 10^{0}$ | $45.0^\circ$ |
| 1 | $2.8000$ | $2.0\times 10^{-1}$ | $18.4^\circ$ |
| 2 | $2.9756$ | $2.4\times 10^{-2}$ | $6.3^\circ$ |
| 3 | $2.9973$ | $2.7\times 10^{-3}$ | $2.1^\circ$ |
| 4 | $2.99970$ | $3.0\times 10^{-4}$ | $0.7^\circ$ |
| 5 | $2.999966$ | $3.4\times 10^{-5}$ | $0.2^\circ$ |

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$:

$$
(A - \mu I) w_k = v_{k-1},
\qquad
v_k = \frac{w_k}{\lVert w_k\rVert_2},
\qquad
\mu_k = v_k^\ast A v_k.
$$

The matrix $A - \mu I$ is fixed, so a single [factorization](/linear-algebra/numerical-linear-algebra/lu-and-cholesky)
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.[^bor-power]

## 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$:

$$
A_k - \mu_k I = Q_k R_k \quad(\text{QR factorization}),
\qquad
A_{k+1} = R_k Q_k + \mu_k I \quad(\text{reverse the product}).
$$

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.[^bor-qr]

$$
% caption: The QR algorithm drives the subdiagonal of $A_k$ to zero; once the
% bottom entry is negligible the eigenvalue is deflated and the algorithm recurses
% on the leading block.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\foreach \i/\x in {0/0, 1/3.4, 2/6.8} {
  \begin{scope}[xshift=\x cm]
  \draw[black] (0,0) rectangle (2.2,2.2);
  \draw[black] (0,2.2)--(2.2,0);
  \end{scope}
}
% k=0: full subdiagonal
\node at (0.55,1.5) {$\ast$}; \node at (1.1,0.95) {$\ast$}; \node at (1.65,0.4) {$\ast$};
\node[acc] at (0.75,0.35) {$\ast$};
\node[anchor=north] at (1.1,-0.15) {$A_0$};
% k=mid
\begin{scope}[xshift=3.4cm]
\node[acc] at (0.75,0.35) {small};
\node[anchor=north] at (1.1,-0.15) {$A_k$: subdiag shrinking};
\end{scope}
% converged
\begin{scope}[xshift=6.8cm]
\node[acc] at (0.75,0.35) {$0$};
\node[acc] at (1.7,0.35) {eig};
\draw[acc, thick] (0,0.7)--(2.2,0.7);
\node[anchor=north] at (1.1,-0.15) {deflate eig};
\end{scope}
\draw[->, black, thick] (2.4,1.1) -- (3.2,1.1);
\draw[->, black, thick] (5.8,1.1) -- (6.6,1.1);
\end{tikzpicture}
$$

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.[^bor-qr]

> **Worked example.** Run the Wilkinson-shifted QR algorithm on the symmetric
>
> $$
> A = \begin{bmatrix} 4 & 1 & 0 \\ 1 & 3 & 1 \\ 0 & 1 & 2 \end{bmatrix},
> $$
>
> whose eigenvalues are $\approx 1.27,\ 3,\ 4.73$. Track the trailing subdiagonal
> entry $(A_k)_{32}$, whose decay to zero deflates the bottom eigenvalue.

| $k$ | trailing subdiagonal $\lvert (A_k)_{32}\rvert$ |
| --- | --- |
| 0 | $1.0\times 10^{0}$ |
| 1 | $9.5\times 10^{-2}$ |
| 2 | $1.2\times 10^{-5}$ |
| 3 | $8.3\times 10^{-18}$ |

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 $\lvert\lambda_2/\lambda_1\rvert$ per step.

$$
% caption: Convergence of the two iterations on a log-error axis. Power iteration
% descends along a straight line at rate theta; the shifted QR algorithm plunges as
% the correct digits triple each step.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (0,-3.4) -- (0,0.6) node[anchor=south] {log error};
\draw[->, black] (0,0) -- (6.9,0) node[anchor=west] {iteration};
\draw[acc, thick] (0.3,-0.15) -- (6.2,-3.0);
\foreach \x/\y in {0.3/-0.15, 1.5/-0.72, 2.7/-1.30, 3.9/-1.87, 5.1/-2.44} {\fill[acc] (\x,\y) circle (1.6pt);}
\node[acc, anchor=west] at (3.5,-1.35) {power iteration};
\draw[black, thick] (0.3,-0.15) .. controls (1.4,-0.4) and (1.9,-0.9) .. (2.4,-3.2);
\foreach \x/\y in {0.3/-0.15, 1.3/-0.55, 2.1/-1.4} {\fill[black] (\x,\y) circle (1.6pt);}
\node[anchor=west] at (2.7,-2.75) {shifted QR};
\end{tikzpicture}
$$

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

> **Two-phase Schur.**
> 1. Reduce $A$ to upper Hessenberg $H = Q^\ast A Q$ by a direct sequence of
>    Householder reflections, at cost $O(m^3)$.
> 2. Run the shifted QR algorithm on $H$, at $O(m^2)$ per step, $O(m^2)$ total
>    after deflations.
>
> The whole computation is $O(m^3)$, comparable to a single QR factorization. For
> self-adjoint $A$, the Hessenberg form is **tridiagonal**, and phase 2 drops to
> $O(m^2)$.[^bor-qr]

$$
% caption: Two phases of a Schur decomposition: a direct $O(m^3)$ Householder
% reduction to Hessenberg (tridiagonal if symmetric), then cheap $O(m^2)$ QR
% sweeps that finish the triangularization.
\begin{tikzpicture}[scale=1.0, font=\footnotesize,
  box/.style={draw, minimum width=24mm, minimum height=13mm, align=center}]
\definecolor{acc}{HTML}{4A6FA5}
\node[box] (a) at (0,0) {full $A$};
\node[box] (h) at (4,0) {Hessenberg $H$};
\node[box] (t) at (8.4,0) {Schur $T$};
\draw[->, black, thick] (a) -- (h) node[midway, above, align=center] {phase 1\\Householder $O(m^3)$};
\draw[->, black, thick] (h) -- (t) node[midway, above, align=center] {phase 2\\QR sweeps $O(m^2)$};
\end{tikzpicture}
$$

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$,

$$
\sigma_i(A) = \sqrt{\lambda_i(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](/linear-algebra/numerical-linear-algebra/qr-and-numerical-least-squares):
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$.[^golub-kahan]
- **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](/linear-algebra/symmetric-quadratic-svd/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(\varepsilon_{\mathrm{mach}})\lVert A\rVert$ of $A$.

$$
% caption: Stable SVD: two-sided Householder bidiagonalization of $A$, then
% implicit QR on the bidiagonal $B$ — the singular values emerge without ever
% forming the condition-squaring product $A^\ast A$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize,
  box/.style={draw, minimum width=26mm, minimum height=12mm, align=center}]
\definecolor{acc}{HTML}{4A6FA5}
\node[box] (a) at (0,0) {$A$};
\node[box] (b) at (3.8,0) {bidiagonal $B$};
\node[box] (s) at (8.2,0) {singular values};
\draw[->, black, thick] (a) -- (b) node[midway, above, align=center] {Householder\\two-sided};
\draw[->, black, thick] (b) -- (s) node[midway, above, align=center] {implicit\\QR};
\end{tikzpicture}
$$

## 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](/linear-algebra/numerical-linear-algebra/lu-and-cholesky)
factorizations triangularize for solving systems; [QR and
Householder](/linear-algebra/numerical-linear-algebra/qr-and-numerical-least-squares)
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](/linear-algebra/numerical-linear-algebra/stability-and-error-analysis) —
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.

[^bor-basic]: **Bornemann**, _Numerical Linear Algebra_, §18 — Basic Concepts: the eigenvalue problem, the instability of the characteristic-polynomial route (the $\operatorname{diag}(1{:}22)$ example), deflation, the Schur decomposition, and unitary diagonalization of normal matrices.
[^bor-pert]: **Bornemann**, §19 — Perturbation Theory: the residual backward error, the Rayleigh quotient as the optimal eigenvalue, the Bauer–Fike bound (normal eigenvalues have $\kappa_{\mathrm{abs}} = 1$), and the Davis–Kahan spectral-gap bound on eigenvector conditioning.
[^bor-power]: **Bornemann**, §20 — Power Iteration: the iteration and its $O(\theta^k)$ / $O(\theta^{2k})$ convergence rates, inverse iteration with a shift, and the resolution of the "ill-conditioned solve" objection via backward stability.
[^bor-qr]: **Bornemann**, §21 — QR Algorithm: the $A_k - \mu_k I = Q_k R_k$, $A_{k+1} = R_k Q_k + \mu_k I$ iteration as disguised inverse iteration, Rayleigh and Wilkinson shifts with quadratic/cubic convergence, and the Hessenberg (tridiagonal) reduction bringing the cost to $O(m^3)$.
[^golub-kahan]: The two-sided bidiagonalization plus implicit-shift QR is the Golub–Kahan SVD algorithm: G. Golub and W. Kahan, "Calculating the Singular Values and Pseudo-Inverse of a Matrix," _SIAM J. Numer. Anal._ 2 (1965), 205–224, https://doi.org/10.1137/0702016.