Eigenvalue and SVD algorithms (QR algorithm, Golub–Kahan) and Krylov methods

Before this lesson: lesson 05 (eigenvalues, power iteration), lesson 06 (spectral theorem, Rayleigh quotient), lesson 07 (positive definite matrices), lesson 12 (QR) and lesson 13 (floating point and conditioning). Time: about 90 minutes, plus the notebook.

---

Why this matters #

Nobody computes eigenvalues as the roots of $\det(A - \lambda I)$. Polynomial roots are very sensitive to rounding in the coefficients, and for $n \ge 5$ there is no formula for them in radicals (Abel–Ruffini). So every eigenvalue algorithm iterates. Two families cover practice. Dense methods (the QR algorithm, Golub–Kahan for the SVD) find all eigenvalues of a matrix you can store, in $O(n^3)$ operations; they run inside np.linalg.eig, eigh and svd. Krylov methods only need the product $v \mapsto Av$. They find a few eigenvalues or solve $Ax = b$ for sparse or implicit matrices with millions of unknowns: the pressure-Poisson step of a fluid solver, the top eigenvalues of a network Hessian. CG also explains why the condition number decides how long such a solve takes.

1. Predict first #

2. The picture #

Power iteration, two upgrades. Power iteration (lesson 05) repeats $x \leftarrow Ax/\|Ax\|$ and finds the dominant eigenvector at rate $|\lambda_2/\lambda_1|$ per step. Two ideas improve it.

• Change the matrix. $(A - \sigma I)^{-1}$ has the same eigenvectors as $A$ and eigenvalues $1/(\lambda_i - \sigma)$. The eigenvalue closest to the shift $\sigma$ becomes the dominant one. This is inverse iteration. If you update $\sigma$ at every step with the current estimate, you get Rayleigh quotient iteration.
• Iterate many vectors at once. Multiply an orthonormal block of $j$ vectors by $A$ and re-orthonormalize it (QR). The block converges to the span of the top $j$ eigenvectors. With $j = n$ this is the QR algorithm.

Interactive figure: loading…

Keep every vector. Power iteration keeps only the last vector $A^{k-1}b$ and throws away $b, Ab, \dots, A^{k-2}b$. A Krylov method keeps all of them and picks the best combination. Every vector in their span is $p(A)b$ for a polynomial $p$ of degree at most $k - 1$. Power iteration always uses $p(\lambda) = \lambda^{k-1}$; a Krylov method chooses the best $p$ for the task. "Best" can be made precise, and that gives the convergence theory below.

3. Definitions and main results #

Inverse iteration and Rayleigh quotient iteration #

Inverse iteration with shift $\sigma$: solve $(A - \sigma I)y = x_k$, set $x_{k+1} = y/\|y\|$. Factor $A - \sigma I$ once (LU), then each step costs two triangular solves, $O(n^2)$; never form the inverse. Let $\lambda_J$ be the eigenvalue closest to $\sigma$ and $\lambda_K$ the next closest. The vector converges to $\pm v_J$ with error $O(|\lambda_J - \sigma|^k / |\lambda_K - \sigma|^k)$: a good shift gives a tiny ratio. The system is nearly singular when $\sigma \approx \lambda_J$, but this is harmless: the large error in $y$ points along $v_J$, which is the direction we want.

Rayleigh quotient iteration (RQI): use the shift $\sigma_k = r(x_k) = x_k^\top A x_k$ (with $\|x_k\| = 1$), updated every step. For symmetric $A$ it converges cubically:

1. If $\|x_k - v\| = O(\epsilon)$, the Rayleigh quotient is accurate to $O(\epsilon^2)$ (lesson 06).
2. One inverse-iteration step with a shift error $O(\epsilon^2)$ multiplies the error by $O(\epsilon^2)$.
3. So $\|x_{k+1} - v\| = O(\epsilon^3)$: the number of correct digits roughly triples per step.

The price is a new factorization at every step. For non-symmetric $A$ the rate is quadratic.

The QR algorithm #

Two vectors at a time, in words. Take $A = \begin{pmatrix} 2 & 1 \\ 1 & 2\end{pmatrix}$. Its eigenvectors are $(1, 1)/\sqrt2$ (eigenvalue 3) and $(-1, 1)/\sqrt2$ (eigenvalue 1). Start with the two vectors $e_1, e_2$.

1. Power step. Multiply both by $A$: $Ae_1 = (2, 1)$ and $Ae_2 = (1, 2)$.
2. Re-orthonormalize with Gram–Schmidt: $q_1 = (2, 1)/\sqrt5$ and $q_2 = (-1, 2)/\sqrt5$. Steps 1 and 2 together are one QR factorization of $A[e_1\; e_2] = A$.
3. Repeat. The first vector does plain power iteration, so it turns toward the dominant eigenvector $(1, 1)/\sqrt2$; the angle error shrinks by about $1/3$ per step.
4. The second vector must stay orthogonal to the first. In $\mathbb{R}^2$ it has only one place to go: the other eigenvector $(-1, 1)/\sqrt2$, up to sign.
5. So $\underline Q = [q_1\; q_2]$ tends to the eigenvector matrix, and $\underline Q^\top A\underline Q$ tends to $\operatorname{diag}(3, 1)$.

In $n$ dimensions the first $j$ vectors do power iteration on a $j$-dimensional subspace, so their span tends to the span of the top $j$ eigenvectors. The QR algorithm does not store these vectors. It updates the matrix $A_k = \underline Q_k^\top A\,\underline Q_k$ directly, as the theorem shows.

Theorem (the QR algorithm is simultaneous power iteration) #

Let $\underline{Q}_k = Q_0Q_1\cdots Q_{k-1}$ and $\underline{R}_k = R_{k-1}\cdots R_1R_0$. Then (a) $A_k = \underline{Q}_k^\top A\,\underline{Q}_k$, and (b) $A^k = \underline{Q}_k\underline{R}_k$.

Proof.

1. From $A_k = Q_kR_k$: $R_k = Q_k^\top A_k$, so $A_{k+1} = R_kQ_k = Q_k^\top A_kQ_k$. Repeating this from $A_0 = A$ gives (a). Each step is an orthogonal similarity: same eigenvalues, and a symmetric $A$ stays symmetric.
2. For (b), use induction. $k = 1$: $A = Q_0R_0$. Assume $A^k = \underline{Q}_k\underline{R}_k$. By (a), $A\underline{Q}_k = \underline{Q}_kA_k$. So $A^{k+1} = A\underline{Q}_k\underline{R}_k = \underline{Q}_kA_k\underline{R}_k = \underline{Q}_kQ_kR_k\underline{R}_k = \underline{Q}_{k+1}\underline{R}_{k+1}$.
3. $\underline{R}_k$ is upper triangular (with a nonzero diagonal if $A$ is nonsingular), so (b) is a QR factorization of $A^k$. The first $j$ columns of $\underline{Q}_k$ are an orthonormal basis of the span of the first $j$ columns of $A^k$, that is, of $A^k\operatorname{span}(e_1, \dots, e_j)$: power iteration on a $j$-dimensional subspace. $\blacksquare$

Convergence. Let $A$ be symmetric with $|\lambda_1| > |\lambda_2| > \dots > |\lambda_n|$, and let all leading principal submatrices of its eigenvector matrix be nonsingular (true for almost every $A$). Then $A_k$ tends to $\operatorname{diag}(\lambda_1, \dots, \lambda_n)$, ordered by decreasing $|\lambda|$, and $\underline{Q}_k$ tends to the eigenvectors (up to signs). The entry $(j+1, j)$ decays like $|\lambda_{j+1}/\lambda_j|^k$: linear convergence, slow when two eigenvalues are close in modulus. For a real non-symmetric $A$ the limit is the real Schur form: upper triangular, with $2\times2$ blocks for complex-conjugate pairs.

The hidden inverse iteration. For symmetric $A$, invert (b): $A^{-k} = \underline{R}_k^{-1}\underline{Q}_k^\top$, and since $A^{-k}$ is symmetric, $A^{-k} = \underline{Q}_k\underline{R}_k^{-\top}$ with $\underline{R}_k^{-\top}$ lower triangular. Its last column gives $A^{-k}e_n \propto \underline{Q}_k e_n$. So the last column of $\underline{Q}_k$ runs inverse iteration. This is why a shift taken from the bottom-right corner works so well.

Answers to section 1, question 2. The off-diagonal entry shrinks by $\lambda_2/\lambda_1 = 1/3$ per step, asymptotically (section 4 shows $1 \to 0.6 \to 0.22 \to 0.074$). For $\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}$ nothing happens: $Q_0 = A$, $R_0 = I$, so $A_1 = A$. The eigenvalues $\pm1$ have equal modulus, and the ratio is 1.

The practical eigenvalue pipeline #

1. Reduce once to Hessenberg form. Householder reflections give $A = QHQ^\top$ with $H$ zero below the first subdiagonal, in about $\tfrac{10}{3}n^3$ flops. A symmetric $A$ becomes tridiagonal (about $\tfrac43 n^3$ flops). Triangular form cannot be reached in a finite number of steps (for $n \ge 5$ it would give the eigenvalues in closed form, see Abel–Ruffini); Hessenberg is the closest form one can reach.
2. QR steps keep the form and are cheap. If $H = QR$ is Hessenberg, then $Q = HR^{-1}$ and $RQ$ are Hessenberg too. One step costs $O(n^2)$ with $n - 1$ Givens rotations instead of $O(n^3)$; on a tridiagonal matrix, $O(n)$.
3. Shifts. $A_k - \mu_kI = Q_kR_k$, $A_{k+1} = R_kQ_k + \mu_kI$. This is still a similarity: $A_{k+1} = Q_k^\top(A_k - \mu_kI)Q_k + \mu_kI = Q_k^\top A_kQ_k$. The Rayleigh shift $\mu_k = (A_k)_{nn}$ is, for symmetric $A$, RQI on the last column (hidden inverse iteration). It fails on $\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}$ ($\mu = 0$, no progress). The Wilkinson shift is the eigenvalue of the trailing $2\times2$ block that is closer to $(A_k)_{nn}$. For symmetric tridiagonal matrices it always converges (Wilkinson, 1968), at least quadratically and in practice cubically.
4. Deflation. When $|h_{n,n-1}| \le \varepsilon\,(|h_{n-1,n-1}| + |h_{nn}|)$, set it to 0. Then $h_{nn}$ is an eigenvalue; continue on the leading $(n-1)\times(n-1)$ block.
5. Real non-symmetric matrices have complex pairs. The Francis double-shift step uses the two eigenvalues of the trailing $2\times2$ block as shifts in one real step, applied implicitly ("bulge chasing").

Interactive figure: loading…

The SVD pipeline: Golub–Kahan #

We want orthogonal $U, V$ with $U^\top AV$ diagonal. The SVD is not a similarity ($U$ and $V$ can differ), so a finite procedure can reach more than Hessenberg form: a bidiagonal matrix $B$, with nonzeros only on the diagonal and the first superdiagonal. Golub–Kahan bidiagonalization does this with Householder reflections from the left and from the right: $A = UBV^\top$.

• Why this is the right form. $B^\top B = V^\top A^\top AV$ is tridiagonal: Golub–Kahan tridiagonalizes $A^\top A$ without forming it.
• Phase 2. The singular values of $B$ are computed by an implicitly shifted QR iteration on $B$ itself, or by divide and conquer.
• Why not $\operatorname{eig}(A^\top A)$. Forming $A^\top A$ squares the condition number (lesson 13). Its computed eigenvalues have absolute errors of about $\varepsilon\sigma_1^2$, so every singular value below about $\sqrt{\varepsilon}\,\sigma_1 \approx 1.5\times10^{-8}\sigma_1$ is lost (section 5 shows $10^{-10}$ coming out as $10^{-8}$).

Golub–Kahan step by step, and which LAPACK routine NumPy calls

For $m \ge n$ the reflections alternate: from the left, zero column $j$ below the diagonal; from the right, zero row $j$ to the right of the superdiagonal. For a $4 \times 3$ matrix ($\times$ = nonzero):

$$\begin{pmatrix}\times&\times&\times\\\times&\times&\times\\\times&\times&\times\\\times&\times&\times\end{pmatrix} \xrightarrow{U_1^\top} \begin{pmatrix}\times&\times&\times\\0&\times&\times\\0&\times&\times\\0&\times&\times\end{pmatrix} \xrightarrow{V_1} \begin{pmatrix}\times&\times&0\\0&\times&\times\\0&\times&\times\\0&\times&\times\end{pmatrix} \xrightarrow{U_2^\top} \begin{pmatrix}\times&\times&0\\0&\times&\times\\0&0&\times\\0&0&\times\end{pmatrix} \xrightarrow{U_3^\top} \begin{pmatrix}\times&\times&0\\0&\times&\times\\0&0&\times\\0&0&0\end{pmatrix}$$

The right reflection $V_1$ acts on columns 2–3 only, so it does not refill column 1. The cost is about $4mn^2 - \tfrac43 n^3$ flops. For $m \gg n$ it is cheaper to compute $A = QR$ first and bidiagonalize the $n \times n$ matrix $R$ (about $2mn^2 + 2n^3$ flops).

Phase 2 is either the implicitly shifted QR iteration on $B$ (Golub–Kahan–Reinsch; LAPACK driver gesvd) or divide and conquer (gesdd, the NumPy and SciPy default). So NumPy's SVD is Golub–Kahan bidiagonalization followed by divide and conquer, not a QR iteration; the QR driver is available as scipy.linalg.svd(..., lapack_driver="gesvd").

For eigenvalues: np.linalg.eig calls LAPACK geev (Hessenberg reduction, then Francis QR); np.linalg.eigh calls syevd (tridiagonal reduction, then divide and conquer); scipy.linalg.eigh uses syevr (tridiagonal reduction, then the MRRR method) by default. The QR algorithm is the core for non-symmetric matrices; for symmetric ones it is one of several choices after the reduction.

Krylov subspaces and the Arnoldi relation #

The Krylov subspace of $A \in \mathbb{R}^{n\times n}$ and $b \in \mathbb{R}^n$ is

Building it needs only matrix-vector products.

Why the plain basis is bad. The vectors $b, Ab, A^2b, \dots$ are the iterates of power iteration. So they all turn toward the dominant eigenvector, and after a few steps they are almost parallel. A basis of almost parallel vectors is extremely ill-conditioned: small rounding errors in the vectors change the coefficients completely.

Interactive figure: loading…

Arnoldi builds an orthonormal basis of the same space, one vector at a time. It orthogonalizes each new vector $Aq_j$ against the old ones (modified Gram–Schmidt, lesson 12):

$q_1 = b/\|b\|$; for $j = 1, \dots, k$: $w = Aq_j$; for $i = 1, \dots, j$: $h_{ij} = q_i^\top w$, $w \leftarrow w - h_{ij}q_i$; then $h_{j+1,j} = \|w\|$, $q_{j+1} = w/h_{j+1,j}$.

The loop as equations, for $k = 2$. Rearrange each step: what was subtracted plus what is left gives back $Aq_j$.

• Step $j = 1$: $Aq_1 = h_{11}q_1 + h_{21}q_2$.
• Step $j = 2$: $Aq_2 = h_{12}q_1 + h_{22}q_2 + h_{32}q_3$.

Put the two equations side by side as columns:

The 0 in the bottom-left corner is there because $Aq_1$ needs only $q_1$ and $q_2$. In general $Aq_j$ needs only $q_1, \dots, q_{j+1}$, so $\tilde H_k$ is zero below its first subdiagonal: it is upper Hessenberg.

Theorem (Arnoldi relation) #

If no $h_{j+1,j}$ is zero, the columns of $Q_k = [q_1 \cdots q_k]$ ($n \times k$) are an orthonormal basis of $\mathcal K_k(A, b)$, and

with $Q_{k+1}$ of size $n\times(k+1)$, $\tilde H_k$ of size $(k+1)\times k$ upper Hessenberg, and $H_k = Q_k^\top AQ_k$ its top $k \times k$ block.

Proof.

1. Step $j$ of the loop, rearranged, reads $Aq_j = \sum_{i=1}^{j+1} h_{ij}q_i$. This is column $j$ of $AQ_k = Q_{k+1}\tilde H_k$, with $h_{ij} = 0$ for $i > j + 1$: upper Hessenberg.
2. By construction $q_{j+1}$ is orthogonal to $q_1, \dots, q_j$ and has norm 1.
3. By induction, $q_{j+1} \in \operatorname{span}(q_1, \dots, q_j, Aq_j) \subseteq \mathcal K_{j+1}$. The $k$ orthonormal vectors lie in $\mathcal K_k$, whose dimension is at most $k$, so they are a basis.
4. Multiply $AQ_k = Q_kH_k + h_{k+1,k}q_{k+1}e_k^\top$ by $Q_k^\top$. Since $Q_k^\top q_{k+1} = 0$: $Q_k^\top AQ_k = H_k$. $\blacksquare$

• Ritz values. $H_k$ is $A$ compressed onto $\mathcal K_k$. Its eigenvalues $\theta$ (Ritz values) approximate eigenvalues of $A$, the extreme ones first. For an eigenpair $H_ky = \theta y$ with $\|y\| = 1$, the residual is $\|AQ_ky - \theta Q_ky\| = h_{k+1,k}|e_k^\top y|$: a free error estimate.
• Breakdown $h_{j+1,j} = 0$ means $\mathcal K_j$ is invariant under $A$: the Ritz values are exact eigenvalues ("lucky" breakdown).
• Cost. $k$ products with $A$, about $2nk^2$ flops of orthogonalization, $nk$ numbers of memory. ARPACK (eigs, eigsh) restarts to bound $k$.
• Lanczos. If $A$ is symmetric, $H_k = Q_k^\top AQ_k$ is symmetric and Hessenberg, so it is tridiagonal: $T_k$ with diagonal $\alpha_j$ and off-diagonal $\beta_j$. The loop collapses to the three-term recurrence $\beta_jq_{j+1} = Aq_j - \alpha_jq_j - \beta_{j-1}q_{j-1}$, with $O(n)$ work per step beyond the product. In floating point the $q_j$ lose orthogonality as soon as a Ritz value converges (Paige), and copies ("ghosts") of converged eigenvalues appear. Full or selective reorthogonalization fixes this.

Interactive figure: loading…

Conjugate gradient #

Let $A$ be SPD, $x^* = A^{-1}b$, error $e = x - x^*$, residual $r = b - Ax = -Ae$, and $\|e\|_A = \sqrt{e^\top Ae}$ (the $A$-norm, or energy norm).

• Energy. $\varphi(x) = \tfrac12 x^\top Ax - b^\top x = \tfrac12\|x - x^*\|_A^2 - \tfrac12 x^{*\top}Ax^*$. Minimizing $\varphi$ is minimizing the $A$-norm error, and $-\nabla\varphi(x) = r$.
• Definition. The CG iterate is $x_k = \operatorname{argmin}\{\|x - x^*\|_A : x \in x_0 + \mathcal K_k(A, r_0)\}$.
• Optimality = orthogonal residual. $x_k$ is this minimizer exactly when $(x_k - x^*)^\top Av = 0$ for all $v \in \mathcal K_k$, that is, $r_k \perp \mathcal K_k$.
• The algorithm finds $x_k$ with short recurrences and one product $Ap$ per step:

The residuals are mutually orthogonal and the directions are $A$-conjugate ($p_i^\top Ap_j = 0$ for $i \ne j$). In an $A$-orthogonal basis the minimization splits into one-dimensional problems, so adding a new direction never changes the old coefficients: progress is never undone. CG is Lanczos applied to $Ax = b$.

• Why SPD. Positive definiteness makes $\|\cdot\|_A$ a norm and $\varphi$ strictly convex; symmetry gives the short recurrence.

CG as a polynomial problem. Every $x_k$ in $x_0 + \mathcal K_k$ has an error of the form $e_k = p(A)e_0$, where $p$ has degree at most $k$ and $p(0) = 1$ (step 1 of the proof below). In the eigenvector basis, $p(A)$ multiplies the component of $e_0$ along $v_i$ by $p(\lambda_i)$. So CG is fast when some polynomial with $p(0) = 1$ is small at every eigenvalue of $A$:

• eigenvalues spread over $[\lambda_{\min}, \lambda_{\max}]$: $p$ must be small on the whole interval, and the best one can do depends on $\kappa = \lambda_{\max}/\lambda_{\min}$ through $\sqrt\kappa$;
• eigenvalues in a few clusters: put one root of $p$ in each cluster;
• only $d$ distinct eigenvalues: put a root on each, and $p$ is exactly 0 on the spectrum.

Interactive figure: loading…

Theorem (CG convergence) #

Proof.

1. $x_k - x_0 \in \mathcal K_k(A, r_0)$ and $r_0 = -Ae_0$, so $x_k - x_0 = -q(A)Ae_0$ with $\deg q \le k - 1$. Then $e_k = (I - Aq(A))e_0 = p(A)e_0$ with $p(\lambda) = 1 - \lambda q(\lambda)$: degree $\le k$, $p(0) = 1$. Every such $p$ arises from some $q$.
2. CG minimizes $\|e_k\|_A$ over this set, which gives the equality.
3. Write $e_0 = \sum_i c_iv_i$ in orthonormal eigenvectors of $A$. Then $\|p(A)e_0\|_A^2 = \sum_i c_i^2\lambda_ip(\lambda_i)^2 \le \max_i p(\lambda_i)^2\sum_i c_i^2\lambda_i = \max_i p(\lambda_i)^2\,\|e_0\|_A^2$.
4. The last inequality: choose one good polynomial (the scaled Chebyshev polynomial) and bound its maximum on $[\lambda_{\min}, \lambda_{\max}]$. Details in the folded proof. $\blacksquare$

Step 4: the Chebyshev polynomial bound

The Chebyshev polynomial $T_k$ has degree $k$, satisfies $|T_k(z)| \le 1$ on $[-1, 1]$, and grows fastest outside $[-1, 1]$ among such polynomials. Map $[\lambda_{\min}, \lambda_{\max}]$ onto $[-1, 1]$ and normalize so that $p(0) = 1$:

$$p(\lambda) = T_k\Big(\tfrac{\lambda_{\max} + \lambda_{\min} - 2\lambda}{\lambda_{\max} - \lambda_{\min}}\Big)\Big/T_k\Big(\tfrac{\kappa + 1}{\kappa - 1}\Big).$$

1. At $\lambda = 0$ the argument is $\tfrac{\lambda_{\max} + \lambda_{\min}}{\lambda_{\max} - \lambda_{\min}} = \tfrac{\kappa + 1}{\kappa - 1}$, so $p(0) = 1$.
2. For $\lambda \in [\lambda_{\min}, \lambda_{\max}]$ the argument lies in $[-1, 1]$, so the numerator is at most 1 in absolute value, and $\max_{\lambda \in \sigma(A)} |p(\lambda)| \le 1/T_k\big(\tfrac{\kappa + 1}{\kappa - 1}\big)$.
3. For $z \ge 1$, $T_k(z) = \tfrac12\big[(z + \sqrt{z^2 - 1})^k + (z - \sqrt{z^2 - 1})^k\big] \ge \tfrac12(z + \sqrt{z^2 - 1})^k$.
4. For $z = \tfrac{\kappa + 1}{\kappa - 1}$: $z^2 - 1 = \tfrac{4\kappa}{(\kappa - 1)^2}$, so $z + \sqrt{z^2 - 1} = \tfrac{\kappa + 1 + 2\sqrt\kappa}{\kappa - 1} = \tfrac{(\sqrt\kappa + 1)^2}{(\sqrt\kappa + 1)(\sqrt\kappa - 1)} = \tfrac{\sqrt\kappa + 1}{\sqrt\kappa - 1}$.
5. Steps 2–4 give $\max|p| \le 2\big(\tfrac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\big)^k$. $\blacksquare$

Consequences:

• Iterations. About $\tfrac12\sqrt\kappa\,\ln(2/\text{tol})$, so $O(\sqrt\kappa)$.
• Steepest descent (exact line search) only reaches $\|e_k\|_A \le \big(\tfrac{\kappa - 1}{\kappa + 1}\big)^k\|e_0\|_A$: about $\tfrac12\kappa\ln(1/\text{tol})$ iterations, $O(\kappa)$. CG replaces $\kappa$ by $\sqrt\kappa$.
• Few distinct eigenvalues. If $A$ has $d$ distinct eigenvalues, $p(\lambda) = \prod_i(1 - \lambda/\lambda_i)$ vanishes on the spectrum, so CG is exact after at most $d$ steps (exact arithmetic). Clustered eigenvalues give fast convergence; $I + (\text{rank } 5)$ needs about 6 steps.
• Preconditioning. With an SPD $M \approx A$ that is cheap to invert (Jacobi, incomplete Cholesky, multigrid), CG on the preconditioned system converges with $\kappa(M^{-1}A)$ instead.

Interactive figure: loading…

Answers to section 1, question 1. CG: $\tfrac12\sqrt{10^4}\ln(2\times10^6) \approx 725$ iterations (the bound; in practice often fewer). Steepest descent: $\tfrac12 \cdot 10^4\ln(10^6) \approx 69{,}000$.

GMRES #

For any nonsingular $A$, GMRES chooses $x_k \in x_0 + \mathcal K_k(A, r_0)$ with the smallest residual $\|b - Ax_k\|_2$. Write $x_k = x_0 + Q_ky$ with the Arnoldi basis ($q_1 = r_0/\beta$, $\beta = \|r_0\|$). Then $b - Ax_k = r_0 - AQ_ky = Q_{k+1}(\beta e_1 - \tilde H_ky)$, and since $Q_{k+1}$ has orthonormal columns, $\|b - Ax_k\| = \|\beta e_1 - \tilde H_ky\|$: a small $(k+1)\times k$ least-squares problem, updated by one Givens rotation per step.

• The residual norms never increase (the spaces are nested).
• No short recurrence: step $k$ costs $O(nk)$ and all $k$ vectors are stored. In practice one restarts every $m$ steps (GMRES($m$)), which can stagnate.
• Symmetric indefinite systems use MINRES (the same minimization with a Lanczos short recurrence).

method	matrix	minimizes over $x_0 + \mathcal K_k$	recurrence	work per step
CG	SPD	$\lVert x - x^*\rVert_A$	short (Lanczos)	1 product + $O(n)$
MINRES	symmetric	$\lVert b - Ax\rVert_2$	short (Lanczos)	1 product + $O(n)$
GMRES	any nonsingular	$\lVert b - Ax\rVert_2$	full (Arnoldi)	1 product + $O(nk)$

The GMRES bound and non-normal matrices

For diagonalizable $A = X\Lambda X^{-1}$, the same polynomial argument as for CG gives

$$\frac{\|r_k\|}{\|r_0\|} \le \kappa(X)\min_{\substack{\deg p \le k\\ p(0) = 1}}\ \max_{\lambda\in\sigma(A)}|p(\lambda)|.$$

The rate depends on where the eigenvalues lie in the complex plane and on the non-normality $\kappa(X)$, not on $\kappa(A)$ alone. Eigenvalues in a disk that does not contain 0 give linear convergence (section 5 shows a disk of radius 0.5 around 1: the residual halves each step). For strongly non-normal matrices ($\kappa(X)$ huge) the eigenvalues can be misleading, and GMRES can stagnate for many steps before it converges.

4. Worked example #

A. QR steps by hand. $A = \begin{pmatrix} 2 & 1 \\ 1 & 2\end{pmatrix}$, eigenvalues 3 and 1.

Gram–Schmidt on the columns: $q_1 = (2, 1)/\sqrt5$, $r_{11} = \sqrt5$, $r_{12} = q_1^\top(1, 2) = 4/\sqrt5$. The remainder $(1, 2) - \tfrac45(2, 1) = (-\tfrac35, \tfrac65)$ has norm $3/\sqrt5$, so $q_2 = (-1, 2)/\sqrt5$:

Check: trace $4$ and determinant $3.36 - 0.36 = 3$ are unchanged. The diagonal moves toward $(3, 1)$. Further unshifted steps give off-diagonal magnitudes $0.6 \to 0.220 \to 0.0740 \to 0.0247$: the ratio approaches $1/3$.

Shifted steps. For $\begin{pmatrix} a & b \\ b & c\end{pmatrix}$ with the Rayleigh shift $\mu = c$: with $d = a - c$ and $\rho = \sqrt{d^2 + b^2}$, the QR factors of $\begin{pmatrix} d & b \\ b & 0\end{pmatrix}$ are $Q = \tfrac1\rho\begin{pmatrix} d & -b \\ b & d\end{pmatrix}$ and $R = \begin{pmatrix}\rho & db/\rho \\ 0 & -b^2/\rho\end{pmatrix}$. The $(2,1)$ entry of $RQ$ is $-b^3/(d^2 + b^2)$. When $b$ is small, the new $|b|$ is about $|b|^3/d^2$: cubic. From $A_1$ ($d = 1.6$, $b = 0.6$): $|b| = 0.6 \to 0.216/2.92 = 0.0740 \to 1.0\times10^{-4} \to 2.6\times10^{-13}$. Three shifted steps reach this level; unshifted steps from $A_1$ would need 26. The exact second iterate is $A_2 = \tfrac{1}{365}\begin{pmatrix}1094 & 27\\27 & 366\end{pmatrix}$. (The Wilkinson shift of a $2\times2$ matrix is an exact eigenvalue, so it deflates in one step.)

B. CG by hand. $A = \begin{pmatrix} 4 & 1 \\ 1 & 3\end{pmatrix}$ (SPD, $\kappa \approx 1.94$), $b = (1, 2)$, $x_0 = 0$. The solution is $x^* = (1/11, 7/11)$, and $\|e_0\|_A = \sqrt{165}/11 \approx 1.168$.

Step 0.

1. $r_0 = b - Ax_0 = (1, 2)$ and $p_0 = r_0 = (1, 2)$.
2. $Ap_0 = (4 + 2,\; 1 + 6) = (6, 7)$.
3. $\alpha_0 = r_0^\top r_0 / p_0^\top Ap_0 = (1 + 4)/(6 + 14) = 5/20 = 1/4$.
4. $x_1 = x_0 + \alpha_0p_0 = (1/4, 1/2)$.
5. $r_1 = r_0 - \alpha_0Ap_0 = (1 - 6/4,\; 2 - 7/4) = (-1/2, 1/4)$. Check: $r_0^\top r_1 = -1/2 + 1/2 = 0$.
6. $\beta_0 = r_1^\top r_1 / r_0^\top r_0 = (1/4 + 1/16)/5 = (5/16)/5 = 1/16$.
7. $p_1 = r_1 + \beta_0p_0 = (-1/2 + 1/16,\; 1/4 + 2/16) = (-7/16, 3/8)$.

Step 1.

1. $Ap_1 = (-28/16 + 6/16,\; -7/16 + 18/16) = (-11/8, 11/16)$.
2. Conjugacy check: $p_0^\top Ap_1 = 1\cdot(-11/8) + 2\cdot(11/16) = 0$.
3. $p_1^\top Ap_1 = (-7/16)(-11/8) + (3/8)(11/16) = 77/128 + 33/128 = 110/128$.
4. $\alpha_1 = r_1^\top r_1 / p_1^\top Ap_1 = (5/16)/(110/128) = 4/11$.
5. $x_2 = x_1 + \alpha_1p_1 = (1/4 - 7/44,\; 1/2 + 3/22) = (1/11, 7/11) = x^*$.
6. $r_2 = r_1 - \alpha_1Ap_1 = (-1/2 + 1/2,\; 1/4 - 1/4) = (0, 0)$.

What to notice:

• Step 0 is the best point on the line $\operatorname{span}\{b\} = \mathcal K_1$: $\alpha_0 = b^\top b/b^\top Ab$. Its error is $\|e_1\|_A \approx 0.337$.
• After $n = 2$ steps CG is exact, because $\mathcal K_2 = \mathbb{R}^2$.
• Steepest descent from the same start also has $x_1 = (0.25, 0.5)$, but then $x_2 = (1/12, 7/12)$ with $\|e_2\|_A \approx 0.097$, and $x_3 = (5/48, 5/8)$ with $\|e_3\|_A \approx 0.028$: it zigzags and is never exact.

5. In code #

SciPy checks the worked example and shows the facts of this lesson:

import numpy as np
from scipy.linalg import hessenberg, schur, svd

# The worked example: one QR step is an orthogonal similarity
A = np.array([[2., 1.], [1., 2.]])
Q0 = np.array([[2., -1.], [1., 2.]]) / np.sqrt(5)       # Gram-Schmidt Q of A
R0 = np.array([[5., 4.], [0., 3.]]) / np.sqrt(5)
A1 = R0 @ Q0
print(np.allclose(Q0 @ R0, A), A1.round(12).tolist())    # True [[2.8, 0.6], [0.6, 1.2]]
print(np.allclose(A1, Q0.T @ A @ Q0), np.linalg.eigvalsh(A1))  # True [1. 3.]

# Phase 1 of eig: Hessenberg form (symmetric input: tridiagonal form)
rng = np.random.default_rng(0)
M = rng.standard_normal((6, 6))
H, Qh = hessenberg(M, calc_q=True)
print(np.abs(np.tril(H, -2)).max(), np.allclose(Qh @ H @ Qh.T, M))   # 0.0 True
T = hessenberg(M + M.T)
print(np.abs(np.triu(T, 2)).max() < 1e-12)                          # True: tridiagonal

# The end point: real Schur form M = Z T Zᵀ, quasi-triangular (2x2 blocks = complex pairs)
Tm, Z = schur(M)
print(np.allclose(Z @ Tm @ Z.T, M), np.round(np.diag(Tm, -1), 3))   # True, nonzeros only in 2x2 blocks

# SVD: two LAPACK drivers agree; eig(BᵀB) loses the small singular value
Qa, _ = np.linalg.qr(rng.standard_normal((50, 4)))
Qb, _ = np.linalg.qr(rng.standard_normal((4, 4)))
B = Qa @ np.diag([1, 1e-3, 1e-6, 1e-10]) @ Qb.T
s_dd = svd(B, compute_uv=False, lapack_driver="gesdd")   # bidiagonalize + divide and conquer
s_qr = svd(B, compute_uv=False, lapack_driver="gesvd")   # bidiagonalize + implicit QR iteration
s_eig = np.sqrt(np.abs(np.linalg.eigvalsh(B.T @ B)))[::-1]
print(s_dd)          # [1.e+00 1.e-03 1.e-06 1.e-10] (the last one to about 7 digits)
print(s_qr)          # the same
print(s_eig)         # last entry 1.0e-08: rounding noise, not 1e-10

Krylov methods, with only products $v \mapsto Av$:

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import cg, gmres, eigsh, LinearOperator

def poisson_2d(n):                       # 5-point Laplacian on an n x n grid, size n² x n²
    T = sp.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(n, n))
    I = sp.identity(n)
    return (sp.kron(I, T) + sp.kron(T, I)).tocsr()

rng = np.random.default_rng(0)

# 1) CG iterations grow like sqrt(kappa): kappa = cot²(pi h / 2) with h = 1/(n+1)
for n in (16, 32, 64, 128):
    A, b, it = poisson_2d(n), rng.standard_normal(n * n), [0]
    x, info = cg(A, b, rtol=1e-8, callback=lambda xk: it.__setitem__(0, it[0] + 1))
    kappa = 1 / np.tan(np.pi / (2 * (n + 1))) ** 2
    print(f"size {n*n:6d}  sqrt(kappa) {np.sqrt(kappa):5.1f}  CG iterations {it[0]:4d}")
# size 256: 10.8, 52 · 1024: 21.0, 101 · 4096: 41.4, 198 · 16384: 82.1, 385  (about 4.8 sqrt(kappa))

# 2) Few distinct eigenvalues: I + (rank 5) has at most 6, so CG needs about 6 steps
N = 2000
U = rng.standard_normal((N, 5))
b, it = rng.standard_normal(N), [0]
x, info = cg(np.eye(N) + U @ U.T, b, rtol=1e-10, callback=lambda xk: it.__setitem__(0, it[0] + 1))
print("I + rank 5:", it[0], "iterations")                     # 6

# 3) Lanczos (ARPACK eigsh) with only products v -> Xᵀ(Xv): PCA without forming XᵀX
m, n = 20000, 1000
X = rng.standard_normal((m, n))
X += rng.standard_normal((m, 3)) @ np.diag([6.0, 4.0, 3.0]) @ np.linalg.qr(rng.standard_normal((n, 3)))[0].T
calls = [0]
def gram(v):
    calls[0] += 1
    return X.T @ (X @ v)
top = eigsh(LinearOperator((n, n), matvec=gram, dtype=float), k=3, which="LA", return_eigenvectors=False)
print(np.allclose(np.sort(top), np.linalg.eigvalsh(X.T @ X)[-3:]), calls[0], "products")   # True 21 products

# 4) GMRES: residuals never increase; rate set by the eigenvalues (here a disk of radius 0.5 around 1)
n = 1000
A = np.eye(n) + 0.5 * rng.standard_normal((n, n)) / np.sqrt(n)
res = []
x, info = gmres(A, rng.standard_normal(n), rtol=1e-10, restart=200, maxiter=1,
                callback=lambda r: res.append(r), callback_type="pr_norm")
print(len(res), "iterations", np.round(res[:4], 3), bool(np.all(np.diff(res) <= 0)))
# 33 iterations [0.447 0.226 0.109 0.055] True: each step about halves the residual

For the 2-D Poisson matrix $\kappa \approx 4/(\pi h)^2$ grows like $h^{-2}$, so $\sqrt\kappa$ and the CG iteration count double when the grid spacing halves. This is why fluid and climate solvers precondition with multigrid, which makes the count almost independent of the grid.

6. Where it shows up #

• Hessian spectra in deep learning. Lanczos with Hessian-vector products (lesson 15) gives the top eigenvalues of a network Hessian without forming it. The largest one bounds stable learning rates ($\eta < 2/\lambda_{\max}$ for gradient descent on a quadratic).
• Gaussian processes. GPyTorch solves $K^{-1}y$ with preconditioned CG and estimates $\log\det K$ with Lanczos, using only kernel matrix-vector products.
• ML for science. Pressure-Poisson solves in incompressible flow (CG with multigrid); implicit time steps of climate and ocean models (Jacobian-free Newton–Krylov: GMRES with Jacobian-vector products only); inverse problems (LSQR, which is Golub–Kahan bidiagonalization built from Krylov vectors); spectral clustering of a sensor graph and the leading EOFs of a huge ensemble (eigsh, svds).

7. Common mistakes #

• Mixing up the QR factorization and the QR algorithm. The algorithm repeats factorizations; one factorization gives no eigenvalues.
• Unshifted QR on a dense matrix. $O(n^3)$ per step and linear convergence. Reduce to Hessenberg/tridiagonal form first, shift, and deflate.
• CG on a matrix that is not SPD. It can stall or diverge. Use GMRES (any matrix), MINRES (symmetric indefinite), or CG on the normal equations only if squaring $\kappa$ is acceptable.
• The plain Krylov basis $b, Ab, A^2b, \dots$ It is numerically rank deficient after a few steps; always orthonormalize (Arnoldi, Lanczos).
• Expecting exact arithmetic. CG does not stop exactly after $n$ steps, and plain Lanczos shows ghost eigenvalues. Stop on the residual (a small residual can hide a large error when $\kappa$ is large, lesson 13).

8. Check yourself #

1. Show that one shifted QR step $A_k - \mu I = Q_kR_k$, $A_{k+1} = R_kQ_k + \mu I$ keeps the eigenvalues.
   Answer

   $R_k = Q_k^\top(A_k - \mu I)$, so $A_{k+1} = Q_k^\top(A_k - \mu I)Q_k + \mu I = Q_k^\top A_kQ_k$: an orthogonal similarity.
2. A symmetric matrix has eigenvalues 10, 9 and 1. In the unshifted QR algorithm, which subdiagonal entry converges slowly, and about how many steps does it need to shrink by $10^{-8}$?
   Answer

   Entry $(2,1)$ decays like $(9/10)^k$: $\ln 10^{-8}/\ln 0.9 \approx 175$ steps. Entry $(3,2)$ decays like $(1/9)^k$: about 9 steps. A shift near 9 or 10 removes the problem.
3. Why is $A$ reduced to Hessenberg form before the QR iteration, and not to triangular form?
   Answer

   The one-time reduction costs $O(n^3)$, after which each QR step costs $O(n^2)$ ($O(n)$ if tridiagonal) and keeps the form. Triangular form would reveal the eigenvalues, which no finite sequence of operations can do for $n \ge 5$.
4. Write the Arnoldi relation with the shapes of all matrices. What changes when $A$ is symmetric?
   Answer

   $AQ_k = Q_{k+1}\tilde H_k$: $A$ is $n\times n$, $Q_k$ is $n\times k$, $Q_{k+1}$ is $n\times(k+1)$ with orthonormal columns, $\tilde H_k$ is $(k+1)\times k$ upper Hessenberg, and $H_k = Q_k^\top AQ_k$. For symmetric $A$, $H_k$ is tridiagonal (Lanczos, three-term recurrence).
5. $\kappa(A) = 10^6$. Use the bounds to compare CG with steepest descent for a $10^{-8}$ reduction of the $A$-norm error.
   Answer

   CG: $\tfrac12\sqrt\kappa\ln(2\times10^8) \approx 9{,}600$ iterations. Steepest descent: $\tfrac12\kappa\ln 10^8 \approx 9.2\times10^6$. A preconditioner that lowers $\kappa$ helps both, CG through $\sqrt\kappa$.
6. An SPD matrix of size 3000 has only the eigenvalues 1, 2 and 3. How many CG steps are needed in exact arithmetic?
   Answer

   At most 3: $p(\lambda) = (1 - \lambda)(1 - \lambda/2)(1 - \lambda/3)$ has degree 3, $p(0) = 1$, and is zero on the spectrum.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Inverse iteration with shift $\sigma$: which eigenvector does it find, and at what rate?

   Show answer

2. How fast does Rayleigh quotient iteration converge for a symmetric matrix, and why?

   Show answer

3. Why does one QR step $A_{k+1} = R_kQ_k$ keep the eigenvalues?

   Show answer

4. Why is $A$ reduced to Hessenberg form before the QR iteration?

   Show answer

5. What is the Wilkinson shift?

   Show answer

6. Define the Krylov subspace $\mathcal K_k(A, b)$.

   Show answer

7. State the Arnoldi relation.

   Show answer

8. What does the Arnoldi matrix $H_k$ become when $A$ is symmetric, and what is the method called?

   Show answer

9. What does CG minimize, and what does GMRES minimize, over $x_0 + \mathcal K_k$?

   Show answer

10. Give the CG convergence bound with $\kappa = \lambda_{\max}/\lambda_{\min}$.

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How does np.linalg.eig work?" — Householder reduction to Hessenberg form, then the shifted (Francis) QR algorithm with deflation; eigh reduces to tridiagonal form and uses divide and conquer; svd bidiagonalizes (Golub–Kahan) and uses divide and conquer (gesdd).
• "Solve a sparse SPD system with $10^6$ unknowns." — CG with a preconditioner (multigrid, incomplete Cholesky): $O(\sqrt{\kappa})$ iterations of one sparse product each, $O(n)$ memory. For a 2-D problem, a sparse Cholesky with a fill-reducing ordering is also good.
• "Top eigenvector of a huge matrix or of a Hessian?" — Lanczos (or power iteration) with matrix-vector (Hessian-vector) products only.
• "CG versus gradient descent?" — on a quadratic, CG picks the optimal polynomial: $\sqrt\kappa$ instead of $\kappa$. Momentum methods (heavy ball, Nesterov) also reach a $\sqrt\kappa$ rate.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Write inverse iteration with one LU factorization, and check its rate $|\lambda_J - \sigma|/|\lambda_K - \sigma|$.
2. Write Rayleigh quotient iteration and watch the number of correct digits triple.
3. Write the unshifted QR algorithm, record the subdiagonal entries, and compare their decay with $|\lambda_{j+1}/\lambda_j|^k$.
4. Write Lanczos with full reorthogonalization, check $AQ_k = Q_kT_k + \beta_kq_{k+1}e_k^\top$, and compare the Ritz values with the extreme eigenvalues of a 1-D Laplacian.
5. Write conjugate gradients, solve a 1-D Poisson problem, and compare its iteration count with steepest descent.

11. Watch and read #

12. Computing Eigenvalues and Singular ValuesGilbert Strang · MIT 18.065, Matrix Methods in Data Analysis (2018) · 49 min

Harvard AM205 video 5.6 - QR algorithmChris Rycroft · Harvard AM205 (numerical methods) · 24 min

Harvard AM205 video 5.9 - Krylov methods: Arnoldi iteration and Lanczos interationChris Rycroft · Harvard AM205 (numerical methods) · 28 min

Harvard AM205 video 5.10 - Conjugate gradient methodChris Rycroft · Harvard AM205 (numerical methods) · 29 min

• Trefethen and Bau, Numerical Linear Algebra (SIAM 1997): Lectures 25–31 (eigenvalue algorithms, Hessenberg reduction, Rayleigh quotient and inverse iteration, QR without and with shifts, computing the SVD) and 32–38, 40 (Arnoldi, GMRES, Lanczos, CG, preconditioning).
• Driscoll and Braun, Fundamentals of Numerical Computation, free online, chapter 8: §8.2 power iteration, §8.3 inverse iteration, §8.4 Krylov subspaces, §8.5 GMRES, §8.6 MINRES and CG and §8.8 preconditioning.
• Jonathan Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" (1994): §6 (steepest descent convergence), §7–8 (conjugate directions, CG), §9 (CG convergence and Chebyshev polynomials).
• Cleve Moler, "The QR Algorithm Computes Eigenvalues and Singular Values" (2019).