Computation and Applications

Chapter 21 closed Part V by naming a debt: “every proof here is an existence or a comparison, and not one of them is an algorithm. That is the subject of Chapter 24.” It is the book’s largest single debt — twenty-six promises across eighteen files — and it has been accumulating since Chapter 3, where partial pivoting was recommended for a reason the book could not yet give. Twenty chapters have built objects and proved things about them. Not one has said how to produce those objects from data of finite accuracy, in finitely many arithmetic operations, and not one has been able to say what an answer computed that way even means.

This chapter answers three questions in that order. What does a computed answer mean? It is the exact answer to a nearby question, and the whole of numerical linear algebra lies in making “nearby” precise: a model of arithmetic, backward error, and the rule that ties the two to Chapter 16’s condition number, with Chapter 20’s residual bounds supplying the error bars. How are the objects of this book actually computed? Elimination with pivoting, least squares by reflections, eigenvalues by iteration — because no finite eigenvalue algorithm can exist — Krylov subspaces, conjugate gradients, the fast Fourier transform and fast multiplication. What is linear algebra for? Vibrations, principal components and graphs, each an application in which a theorem already proved does the work.

The chapter is about arithmetic, algorithms and cost. No product, library or language is named anywhere in it: an algorithm here is a rule for producing numbers from numbers, and its cost is a count of arithmetic operations.

Conventions

Seven conventions hold throughout, and each is stated once where it is introduced.

$u$ is the unit roundoff and $\operatorname{fl}(\cdot)$ the rounding map of §01’s model, with $\mathcal{R}$ the set of representable numbers and $\gamma_k = ku/(1-ku)$ the accumulated-rounding constant, defined in §01 and used in every section that mentions arithmetic.
A flop is one arithmetic operation, and costs are leading-order counts, written $\tfrac23 n^3 + O(n^2)$ . The constant matters and the exponent matters more; §§02 and 13 are where that sentence is made quantitative.
Perturbations of data are $\Delta\mathbf{A}$ and $\Delta\mathbf{b}$ , and a computed object carries a hat: $\widehat{\mathbf{x}}$ , $\widehat{\mathbf{L}}$ . Chapters 16, 17 and 20 wrote a perturbation $\mathbf{E}$ ; $\mathbf{E}$ is kept here for a perturbation whose size is given, and $\Delta$ for one produced by an algorithm.
$\kappa(\mathbf{A})$ always means Chapter 16 §08’s condition number (Definition 16.8.1). Chapter 20 §11’s condition number of a simple eigenvalue does not appear in this chapter, so the two never have to be told apart.
Iteration indices are superscripts in parentheses, $\mathbf{x}^{(k)}$ , $\mathbf{A}^{(k)}$ , leaving subscripts free for coordinates.
$\lvert\mathbf{A}\rvert$ is the entrywise absolute value $(\lvert a_{ij}\rvert)$ , as in Chapter 19, and $\mathbf{A} \le \mathbf{B}$ the entrywise order. This is a local override of Chapter 13’s $\lvert\mathbf{A}\rvert = (\mathbf{A}^{*}\mathbf{A})^{1/2}$ , which does not appear in this chapter.
New symbols, each defined at its first use: the growth factor $g_n$ (§02), the Krylov space $\mathcal{K}_k(\mathbf{A}, \mathbf{b})$ and the Ritz values (§06), the energy norm $\left\| \cdot \right\|_{\mathbf{A}}$ and the Chebyshev polynomials $T_k$ (§07, the same polynomials as Chapter 11 §10), the mass and stiffness matrices $\mathbf{M}$ and $\mathbf{K}$ (§09), the data matrix $\mathbf{X}$ with its centered form $\mathbf{X}_c$ and sample covariance $\mathbf{S}$ (§10), the graph Laplacian and its reduced form (§11, following Chapter 3 §07), the unnormalized Fourier matrix $\mathbf{W}_n = \sqrt n\,\mathbf{F}_n$ (§12), and the arithmetic cost $M(n)$ with the exponent $\omega$ (§13). Two letters are reused, each with a warning where it happens: $\mathbf{L}$ is the flowed matrix of §08 and the Laplacian of §11, and $\mathbf{K}$ is §09’s stiffness matrix, not Chapter 19 §09’s second-difference matrix $\mathbf{K}_n$ .

What the chapter assumes

Every earlier chapter assumed only field axioms, and Part V added a short list of facts from analysis. This chapter adds assumptions of a different kind, and both are stated as definitions so that the reader can see exactly where mathematics stops.

Floating-point arithmetic is not a theorem of linear algebra. It is a model of a machine. A machine stores finitely many numbers and rounds every result to one of them, and the statement that it does so with a small relative error is a description of a device, not a fact to be proved. So §01 states it as a definition, Definition 24.1.1, in the register this book uses for definitions: a finite set of representable numbers, a rounding map $\operatorname{fl}$ , and the requirement $\operatorname{fl}(x \odot y) = (x \odot y)(1 + \delta), \qquad \lvert\delta\rvert \le u,$ for $\odot \in \{+, -, \times, \div\}$ and for the square root, with overflow and underflow excluded by assumption. Everything after §01 is a theorem about that model. The model is an idealization: the arithmetic in common use satisfies it for those five operations, with $u$ about $10^{-16}$ , and every rounding claim in this chapter is a consequence of it and of nothing else.

Section 9 does the same for a piece of physics. Newton’s second law and Hooke’s law are not theorems either, so Definition 24.9.1 names them as clauses (P1) and (P2), says that Hooke’s law holds only for small displacements and that no real system is free of friction, and every statement about vibrations is then a theorem about that idealization. Nothing outside §09 uses it.

No new analysis is imported. The list remains Chapter 16’s facts (A1)–(A6) — completeness, monotone convergence, compactness of closed bounded sets, the extreme value theorem, the intermediate value theorem and the mean value theorem — together with the algebra of limits. Each is named again where it is used. Chapter 21’s quoted fact (A7), the Nevanlinna–Pick characterization behind Loewner’s theorem, is not used anywhere in this chapter, and neither is any consequence of it that Chapter 21 marked as resting on it.

What is quoted and not proved

Several classical results in this subject have proofs longer than the sections that need them. The book’s practice is to prove everything or to name the result at its point of use and say plainly that it is not proved. Here is the complete list, section by section; each item is flagged again in the section that quotes it.

§02 — that Chapter 3’s ordering of elimination, which updates the whole trailing block at each step, satisfies the same backward error bound as the ordering §02 displays. The section also declines to quantify the seam between that bound, which is about the factors the model computes, and the growth bound, which is about exact elimination. Separately, and not a quoted theorem at all: partial pivoting is not shown to be backward stable, and §02 says so rather than tidying it.
§03 — Wilkinson’s accumulated backward error for a complete Householder factorization, of which §03 proves the one-reflection case; and the rates at which classical and modified Gram–Schmidt lose orthogonality, $u\kappa_2(\mathbf{A})^2$ and $u\kappa_2(\mathbf{A})$ .
§04 — that Rayleigh quotient iteration converges at a quadratic rate, not a cubic one, for a general diagonalizable matrix. The cubic Hermitian rate is proved.
§05 — the theorem of Abel and Ruffini, quoted from Chapter 11 §08, which is the reason every eigenvalue method in this chapter iterates; and the general convergence theorem of the QR algorithm, of which §05 proves the first column and hence the whole $2 \times 2$ case.
§06 — the Kaniel–Paige–Saad estimate for the extreme Ritz values, of which §06 proves a clean special case; and the analysis of the loss of orthogonality in Lanczos, which §06 describes and names without proving.
§08 — the differentiability of $t \mapsto \mathbf{Q}(t)$ , which the continuous form of the Toda statement rests on. The discrete form, that the QR iterates sit on the flow at integer times, is proved.
§10 — not a quoted theorem but a missing subject. Chapter 13 §10 promised an argument that a sample says something about a population, and that is statistics, which this book does not develop. Section 10 proves everything deterministic, supplies the finite-population reading, and says plainly where it stops.
§13 — the recursive form of Strassen’s error bound; Winograd’s lower bound of seven multiplications for the $2 \times 2$ product; the equality of the two definitions of $\omega$ , through all arithmetic operations and through multiplications alone; every bound on $\omega$ below $\log_2 7$ ; and Laderman’s twenty-three-product identity for $3 \times 3$ matrices, which an exercise hands to the reader.

Nothing in this book is proved from any of them. Every theorem in the chapter either proves its own case or carries what it needs as an explicit hypothesis, and since all of these enter the book here, no earlier chapter can rest on them either. The one place a quoted result is used at all is exercise C1 (b) of §13, where Laderman’s identity is given to the reader as a hypothesis and the deduction from it is the reader’s own.

What you need

This chapter draws on almost the whole book. Chapters 5, 14, 15, 22 and 23 are the exceptions: duality, bilinear forms, tensor products, affine and projective geometry, and algebras and representations are not used here. Chapter 21 supplies the sentence the chapter answers, but none of its theorems.

Chapter 0 §06 and §09 — the complex numbers §07 and §12 work in, and the definition of the matrix product with its three views and the trace, which between them supply every operation count in the chapter: Theorem 0.6.2, Theorem 0.6.4, Example (Roots of Unity), Definition 0.9.9, Example (A product of square matrices), Theorem 0.9.2, Theorem 0.9.6.
Chapter 1 §06 — basis extension, which §06 uses to finish an orthonormal list, and “a subspace of the full dimension is the whole space”, which §11 uses to identify a column space: Theorem 1.6.5, Theorem 1.6.8.
Chapter 3 — the invertibility criteria, rank, elimination and the LU factorization, and the incidence matrix of a network that §11 turns into the Laplacian: Theorem 3.4.4, Theorem 3.5.7, Theorem 3.5.5, Lemma 3.6.3, Lemma 3.6.1, Lemma 3.6.2, Theorem 3.6.4, Theorem 3.6.5, Example (An invertible matrix with no LU factorization), Definition 3.7.1, Proposition 3.7.4, Corollary 3.7.5, Exercise (C3: Networks in several pieces).
Chapter 4 §02 — similarity, the relation every eigenvalue iteration in §§04–08 preserves: Definition 4.2.2.
Chapter 6 §01 — that a polynomial over an infinite field is determined by the function it defines, which is how §08 turns an equality of values into an equality of characteristic polynomials: Theorem 6.1.5.
Chapter 7 — determinants, almost all of them for §11: the alternating properties, the definition with the triangular and transpose rules, multiplicativity and the invertibility criterion, Laplace expansion and the adjugate identity, minors and rank, and the Gram determinant that §11 turns into the matrix-tree theorem. Section 6 also uses the tridiagonal recurrence and the Vandermonde determinant, and §§05, 07 and 08 the characteristic polynomial: Theorem 7.2.1, Definition 7.4.1, Theorem 7.4.5, Theorem 7.4.6, Theorem 7.5.4, Theorem 7.5.2, Theorem 7.6.1, Theorem 7.6.2, Theorem 7.7.1, Corollary 7.7.3, Exercise (C1: Counting spanning trees), Theorem 7.8.3, Theorem 7.8.1, Theorem 7.10.2, Theorem 7.10.3.
Chapter 8 — block multiplication, which is what makes both fast algorithms recurse, and Schur complements, which are what elimination is: Theorem 8.1.1, Definition 8.3.1, Theorem 8.3.1, Theorem 8.3.2.
Chapter 9 — diagonalization for §04, the eigenvalues on the diagonal of a triangular form for §05, and Cayley–Hamilton for §08: Theorem 9.4.1, Theorem 9.7.4, Theorem 9.6.1.
Chapter 10 — cyclic subspaces, which are Krylov subspaces under another name, and the matrix exponential, which §08 flows along and §09 solves with: Definition 10.5.1, Definition 10.5.4, Definition 10.5.2, Definition 10.5.3, Proposition 10.5.1, Theorem 10.5.2, Theorem 10.5.6, Definition 10.9.2, Theorem 10.9.4, Corollary 10.9.5, Theorem 10.9.6, Example (A Second-Order Equation with a Repeated Root).
Chapter 11 — the chapter’s largest single source: inner products and Pythagoras, orthonormal bases and Gram–Schmidt, projections and best approximation, least squares, self-adjointness and isometries, reflections and rotations with the QR factorization and Hessenberg form, Parseval’s identity, and the orthogonal polynomials behind the three-term recurrence: Theorem 11.1.5, Corollary 11.1.6, Theorem 11.1.3, Theorem 11.1.4, Definition 11.1.3, Theorem 11.2.2, Theorem 11.2.1, Definition 11.3.2, Definition 11.3.3, Theorem 11.3.2, Theorem 11.3.5, Theorem 11.4.2, Corollary 11.4.4, Proposition 11.6.7, Definition 11.7.2, Proposition 11.7.4, Theorem 11.7.3, Definition 11.8.1, Proposition 11.8.2, Theorem 11.8.3, Theorem 11.8.1, Theorem 11.8.4, Definition 11.8.2, Theorem 11.8.5, Theorem 11.9.2, Definition 11.10.1, Proposition 11.10.1, Theorem 11.10.3.
Chapter 12 — the spectral theorem in both forms, and the circulants and the discrete Fourier transform that §§11 and 12 rest on: Corollary 12.4.3, Corollary 12.5.2, Definition 12.9.1, Theorem 12.9.4, Corollary 12.9.5, Theorem 12.9.3, Theorem 12.9.6, Example (The Fourier matrix of size four), Exercise (B1: Eigenvalues of small circulants).
Chapter 13 — positive definiteness and Cholesky, simultaneous congruence, which is exactly §09’s basis of normal modes, the Loewner order, the singular value decomposition with low-rank approximation for §10, and the pseudoinverse: Definition 13.1.1, Theorem 13.1.5, Theorem 13.1.4, Proposition 13.1.2, Theorem 13.2.3, Example (A Cholesky Factorization, and One That Fails), Definition 13.4.1, Theorem 13.4.2, Theorem 13.4.1, Example (A pencil that does not commute), Definition 13.5.1, Definition 13.8.1, Theorem 13.8.1, Definition 13.10.1, Theorem 13.10.3, Corollary 13.10.4, Definition 13.11.1.
Chapter 16 — norms and operator norms, the derivative of a matrix-valued function, and above all the condition number and the theorem that bounds a relative error by it: Lemma 16.1.1, Theorem 16.3.3, Theorem 16.3.2, Definition 16.6.1, Lemma 16.6.4, Theorem 16.6.5, Theorem 16.6.6, Definition 16.8.1, Proposition 16.8.2, Theorem 16.8.3, Theorem 16.8.4.
Chapter 17 — the Rayleigh quotient and its critical points, the min–max description of the eigenvalues, Weyl’s inequalities, interlacing, and the variational description of the singular values: Definition 17.1.1, Proposition 17.1.1, Proposition 17.1.2, Proposition 17.1.3, Theorem 17.2.3, Corollary 17.2.4, Corollary 17.3.6, Lemma 17.3.5, Theorem 17.4.2, Corollary 17.4.3, Lemma 17.4.5, Theorem 17.9.4, Corollary 17.9.5.
Chapter 18 — two small facts, one from each end of the chapter: the coset description of an affine hull, which §10 minimizes over, and the logarithm of a positive definite matrix, which §08 needs to know which matrices its theorem covers: Proposition 18.1.4, Lemma 18.11.1. Both are cited in the text, in §10’s opening and §08’s warning, and neither is a step in a proof.
Chapter 19 — the cyclic permutation that defeats the power method, the Google matrix and the PageRank rate that §04 explains, and the second-difference matrix that §09’s uniform chain turns out to be: Example (The Cyclic Permutation), Theorem 19.8.2, Proposition 19.8.3, Proposition 19.9.6. These are examples the chapter re-reads rather than results it argues from: they appear in §04’s and §09’s examples, and in no proof.
Chapter 20 — the perturbation bounds that turn a residual into an error bar, which is what makes every iteration in this chapter stoppable: Theorem 20.4.4, Theorem 20.6.1, Theorem 20.6.2. They are quoted where they are needed — the warnings of §§01 and 06 — and applied in the examples and exercises of §§04, 05 and 06, but no proof here is built on one.

Roadmap

Floating point and backward error. The model, the constant $\gamma_k$ , the backward error of an inner product, and the three words the chapter runs on: backward stable, forward error, mixed stability. Then the rule, made a theorem: forward error at most condition number times backward error.
Elimination, pivoting and cost. One step of elimination is a Schur complement of a $1 \times 1$ pivot. The cost is $\tfrac23 n^3$ , and half that for Cholesky. A triangular solve is backward stable, and so is the factorization, up to a growth factor that partial pivoting bounds by $2^{n-1}$ — a bound that is attained, which is why the section stops short of calling partial pivoting stable.
Least squares in practice. The sensitivity of the least-squares solution involves the residual as well as $\kappa_2(\mathbf{A})$ . Forming the normal equations loses twice the digits; the reflection route loses them once, because one reflection is backward stable, and classical Gram–Schmidt loses orthogonality at a rate this section states with credit.
The power method. Repeated multiplication finds the dominant eigenvector at rate $\lvert\lambda_2/\lambda_1\rvert$ , the Rayleigh quotient squares that for Hermitian matrices, a shift and an inverse pick out any eigenvalue, and updating the shift from the Rayleigh quotient makes the convergence cubic. It all fails when two eigenvalues share a modulus.
The QR algorithm. Factor and reverse, repeatedly. The iterates are orthogonal iteration applied to the powers of the matrix, which is where the convergence comes from; Hessenberg form brings the cost of a step down to about $6n^2$ flops; and a shift is what turns slow convergence into fast.
Krylov subspaces. The right home for an iteration that can afford only one matrix–vector product at a time. Gram–Schmidt inside that space is the Arnoldi process and gives a Hessenberg matrix, with a three-term recurrence — the Lanczos process — when the matrix is Hermitian; the eigenvalues of the small matrix interlace the large one’s, and Chapter 20’s bounds apply to them verbatim. In floating point the three-term economy costs orthogonality, and the section says so.
Steepest descent and conjugate gradients. Solving $\mathbf{A}\mathbf{x} = \mathbf{b}$ by minimizing a quadratic. Steepest descent converges like $(\kappa-1)/(\kappa+1)$ ; conjugate gradients minimizes over the whole Krylov space and converges like $(\sqrt\kappa - 1)/(\sqrt\kappa + 1)$ . The square root is the point, and Chebyshev polynomials, small on an interval and large just outside it, are where it comes from.
Lax pairs and the Toda flow. The chapter’s one piece of pure mathematics: a differential equation whose flow preserves the spectrum, and which the QR iterates sit on at integer times. That the flow’s time-one map is a step of the algorithm is derived, but only after granting a differentiability this book has not set up. An algorithm and a differential equation can be the same object seen at two resolutions.
Vibrations. Masses and springs give $\mathbf{M}\ddot{\mathbf{x}} = -\mathbf{K}\mathbf{x}$ , the substitution $\mathbf{x} = e^{i\omega t}\mathbf{v}$ gives a generalized eigenproblem, and Chapter 13 gives the $\mathbf{M}$ -orthonormal basis of normal modes. Stiffening raises every frequency, adding mass lowers them, and pinning a mass interlaces them.
Principal components. The best-fitting $k$ -dimensional affine subspace of a finite set of points is spanned by the leading singular vectors of the centered data matrix, and “variance explained” is a theorem rather than a convention. The covariance matrix is never formed, for §03’s reason, and the inference from a sample to a population is handed back to statistics.
Graph Laplacians. $\mathbf{L} = \mathbf{N}\mathbf{N}^{\top}$ for the incidence matrix of Chapter 3 §07. The multiplicity of the eigenvalue $0$ counts components, the second smallest eigenvalue measures connectivity and the signs of its eigenvector split the graph, and the determinant of a reduced Laplacian counts spanning trees, which completes Chapter 7 §07.
The fast Fourier transform. The Fourier matrix factors into a permutation, two half-size copies of itself and a diagonal factor. Solving the resulting recurrence gives $O(n\log n)$ , and with it fast convolution, fast polynomial multiplication and fast circulant solves.
Fast matrix multiplication. Two $2 \times 2$ matrices multiply with seven products instead of eight, and the recursion gives $O(n^{\log_2 7})$ . The exponent $\omega$ is at least $2$ , which is proved; everything known about it below $\log_2 7$ is stated with credit and not proved. Strassen’s algorithm satisfies a norm error bound but no entrywise one, and a small entry really can come out with every digit wrong. The section closes the chapter.

Named moves

State the machine as a definition, then prove theorems about it. When a subject rests on something that is not mathematics, isolate it in one definition, say what kind of statement it is, and derive the rest. Everything in this chapter that mentions rounding is a consequence of Definition 24.1.1 alone, and §09 treats its two physical laws the same way.
Measure an answer by a change in the data that makes it exact. The forward error is what one wants and cannot see; the backward error is what one can bound. Prove the backward statement and let the condition number convert it.
Count the operations, and solve the recurrence exactly. An asymptotic order hides the constant that decides whether an algorithm is worth using. Every cost in this chapter is a closed form or an explicit leading term.
A residual is an error bar. Chapter 20’s perturbation bounds turn a quantity the algorithm already has into a certificate, which is what allows an iteration to know when to stop.
Iterate in a subspace, not along a direction. The Krylov space is generated by the one operation an iteration can afford, and it is large enough to contain the answer to a polynomial approximation problem whose solution is the algorithm.
A conserved spectrum comes from a commutator. If $\dot{\mathbf{L}} = [\mathbf{B}, \mathbf{L}]$ , then $\mathbf{L}(t)$ is a similarity orbit and its eigenvalues do not move; differentiating a conjugation is the whole proof.
A fast algorithm is a factorization. Factor the linear map (the Fourier matrix) or the bilinear map (the matrix product) into cheap pieces and recursive calls, then count the calls against the size reduction.