24.13  Fast Matrix Multiplication

Every algorithm in this chapter is built out of matrix products, and every cost we have counted has been counted against them. Section 12 broke one such count: the Fourier matrix is applied in about $n\log n$ operations, although the definition of a matrix–vector product asks for $n^2$. The definition of a matrix product asks for $n^3$ multiplications. This section asks whether $n^3$ is necessary, answers no, and then asks what the answer costs in accuracy.

24.13.1. What a fast algorithm would be

Count the operations that Definition 0.9.9 prescribes. The product $\mathbf{A}\mathbf{B}$ of two matrices in $M_n(F)$ has $n^2$ entries, each of them a sum of $n$ products, so the recipe performs $n^3$ multiplications and $n^2(n-1)$ additions: $$n^3+n^2(n-1)=2n^3-n^2$$ arithmetic operations in all. Chapter 3 §06 made the same kind of count for an LU factorization and got $\frac{2}{3}{n}^3$, so a product costs about three times a factorization. Nothing so far has suggested that either number could be smaller.

Before asking whether it can be, we need to say what “it” is. A claim that some procedure multiplies matrices in fewer than $n^3$ multiplications is empty until the permitted procedures are fixed, and the fixing has to be done carefully: a procedure allowed to compare two numbers, or to branch on whether an entry is zero, is a different object from one that only computes.

A program is a list of arithmetic instructions, each using the inputs or earlier results, and its cost is the length of the list.

In words: input steps are free, because they only name the data; the three arithmetic operations each cost $1$; and the program must produce all $n^2$ entries as polynomial identities, not merely for the particular numbers at hand. There are no comparisons, no branches, and no loops — a straight-line program does the same thing whatever the data are. Division is excluded, which costs nothing here: none of the algorithms below uses it.

Two examples fix the definition. For $n=1$ the single instruction $g_3=g_1{g}_2$ computes $a_{11}{b}_{11}$, so $M(1)=1$. For $n=2$ the recipe of Definition 0.9.9 is the program that forms the eight products $a_{ik}{b}_{kj}$ and then four sums, at a cost of $12$; so $M(2)\le 12$. A degenerate case worth naming: a program that lists only input steps has cost $0$ and computes nothing, since no entry of $\mathbf{A}\mathbf{B}$ is an indeterminate.

The exponent of matrix multiplication is the smallest growth rate a family of such programs can have.

The infimum is over all real $\tau$ for which some constant $c$ and some $n_0$ give $M(n)\le c{n}^{\tau }$ for every $n\ge n_0$, so $\omega$ is a property of the whole family of problems, not of any one algorithm. The set is nonempty, since $\tau =3$ belongs to it by the count above. It is bounded below, and the reason is the only thing about $\omega$ that this book proves.

So $2\le \omega \le 3$, and the question of the section is where in that interval $\omega$ lies.

24.13.2. Seven products instead of eight

The classical program for $n=2$ forms the eight products $a_{ik}{b}_{kj}$. Its cost is what makes the general recipe cost $n^3$, because a $2\times 2$ block form turns one product of size $n$ into eight of size $n/2$ — and $8=2^3$ is exactly the recursion that reproduces $n^3$. Reducing eight to seven would change the exponent.

Two features of that proof matter more than the result. First, it used only addition, subtraction, multiplication, and their laws in $M_m(F)$ — never the commutativity of a product, never an inverse, never a division. So the identity would hold with the blocks taken in any ring whatever. Second, the case $m=1$ is the ordinary $2\times 2$ product with seven scalar multiplications instead of eight, and the case $m>1$ is the same statement one level up. That self-similarity is what the next theorem exploits.

24.13.3. The recursion, and its exact cost

Apply Theorem 24.13.2 with $m\times m$ blocks, and each of the seven products is a product of two $m\times m$ matrices, to be computed by the same rule. The recursion bottoms out at $1\times 1$ blocks, where a product is one multiplication.

The number ${\log }_{2}7=2.8073\dots$ is below $3$, and that is the whole content: an operation that looked like an $n^3$ operation is not one. The mechanism is worth naming, because Section 12 used it too. There, the Fourier matrix ${\mathbf{F}}_n$ was factored into a permutation, two copies of ${\mathbf{F}}_{n/2}$, and a diagonal factor, and the recurrence $C(n)=2C(n/2)+O(n)$ followed. Here the bilinear map $(\mathbf{A},\mathbf{B})\mapsto \mathbf{A}\mathbf{B}$ is factored into seven products and eighteen sums, and the recurrence $T(n)=7T(n/2)+\frac{9}{2}{n}^2$ follows. In both cases the saving comes from a factorization, and the exponent is decided by the number of recursive calls against the size reduction: $2$ calls at half size against a linear overhead gives $n\log n$; $7$ calls at half size against a quadratic overhead gives $n^{{\log }_{2}7}$.

24.13.4. What is known about the exponent, and what is not

Theorem 24.13.3 gives $\omega \le {\log }_{2}7<2.8074$, and Proposition 24.13.1 gives $\omega \ge 2$. Both are proved above. The rest of this subsection is stated and not proved here, and nothing in this book depends on any of it.

Strassen’s bound was the first below $3$, and it has been improved many times since. The decisive technique is the laser method of Coppersmith and Winograd, which brought the bound below $2.376$; successive refinements of that method by several authors have since pushed the published bound below $2.372$. No lower bound better than the $\omega \ge 2$ proved above is known. Whether $\omega =2$ is an open problem.

Two honest caveats belong with those numbers.

• The bounds below ${\log }_{2}7$ are obtained by analyzing families of recursive constructions, not by exhibiting seven-product identities; the algorithms they yield carry constants so large that at every size that arises in practice they are slower than the $2n^3-n^2$ recipe. This is a known feature of the constructions, and it is recorded here, not proved.
• $\omega$ may depend on the field $F$, and the definition above was made field by field for that reason. The bounds just quoted, and Strassen’s, hold over every field.

24.13.5. Accuracy: a norm bound, but not an entrywise one

Everything so far was exact algebra. In the floating-point model of Definition 24.1.1 the two algorithms are not interchangeable, and the difference is sharper than “one is less accurate”.

Start with the classical recipe, which inherits Section 1’s analysis directly. Write $|\mathbf{A}|$ for the entrywise absolute value $(|a_{ij}|)$, as in Chapter 19, and $\hat{\mathbf{C}}$ for the computed product. Everything below is real, because Definition 24.1.1 is.

The bound is entrywise, and that is its strength: the error in the $(i,j)$ entry is measured against the size of the data that produced that entry. If $(|\mathbf{A}||\mathbf{B}|{)}_{ij}$ is tiny, the error there is tiny too. In particular, if $\mathbf{A}$ and $\mathbf{B}$ have nonnegative entries, then $(|\mathbf{A}||\mathbf{B}|{)}_{ij}=(\mathbf{A}\mathbf{B}{)}_{ij}$ and every entry is computed with small relative error.

Strassen’s algorithm satisfies a weaker bound, in which the whole of $\mathbf{A}$ and the whole of $\mathbf{B}$ appear on the right. Here is the base case, proved.

A bound of this shape survives the recursion, with a constant that grows with $n$: for $n=2^p$ there is a constant $c(n)$, polynomial in $n$, with $$\underset{i,j}{max}|{\hat{c}}_{ij}-(\mathbf{A}\mathbf{B}{)}_{ij}|\le c(n)u\underset{i,j}{max}|a_{ij}|\underset{i,j}{max}|b_{ij}|+O(u^2).$$ This is not proved here. A bound of this form was first established by Brent, and the constant has been improved since; the induction is longer than this section, and nothing below uses it. What is proved below is that no improvement of it can be entrywise.

The constant $1/u$ is about ${10}^{16}$ in the arithmetic in common use, so the conclusion is that no entrywise bound worth stating exists: Proposition 24.13.4 holds with $c={\gamma }_2/u$, which is close to $2$, and for Strassen’s algorithm the same inequality at $n=2$ needs a constant larger than the reciprocal of the unit roundoff. The hypotheses of the theorem are not exotic.

They are not vacuous either, and the fourth one fixes the scale of $\epsilon$. Since $1+\epsilon$ is representable, Definition 24.1.1 writes $\mathrm{fl}((1+\epsilon {)}^2)=(1+\epsilon {)}^2(1+\delta )$ with $|\delta |\le u$, and the fourth hypothesis says the left side is $1+2\epsilon$, which forces $$\delta =-\frac{{\epsilon }^2}{(1+\epsilon {)}^2},\text{hence}{\epsilon }^2\le u(1+\epsilon {)}^2.$$ So $\epsilon$ is at most about $\sqrt{u}$: the theorem asks for a quantity small enough that its square falls off the end of the significand, and for nothing else. The next example produces one.

24.13.6. Whether to use it

The operation count and the error analysis pull in opposite directions, and the honest summary is a list of exact numbers rather than a verdict. Suppose the recursion is stopped at a base size $n_0$, below which the classical recipe is used, and let $n=n_0{2}^p$. The same induction as in Theorem 24.13.3, started from $T(n_0)=2n_0^3-n_0^2$, gives $$T(n_0{2}^p)=7^p(2n_0^3-n_0^2)+6n_0^2(7^p-4^p).$$ Three readings of this formula. With $n_0=1$ the recursion runs to scalars, and it is slower than the classical recipe until $n$ exceeds about $654$; at $n=512$ it costs $280902385$ operations against $268173312$, and at $n=1024$ it costs $1971035287$ against $2146435072$. With a classical base case the picture changes completely: a single level of the recursion, $p=1$, costs $14n_0^3+11n_0^2$ against the classical $16n_0^3-4n_0^2$ at size $2n_0$, and the first wins as soon as $2n_0^3>15n_0^2$, that is from $n_0=8$, so from $n=16$. And at $n=1024$ with $n_0=64$ the count is $1301696512$, about $0.61$ of the classical one.

So on operation count alone the trade is favorable at sizes that are not large. Three things keep the classical algorithm the one in ordinary use at small and moderate sizes, and none of them is visible in an operation count.

• Additions are not free in the way the count assumes. The seven products must be stored, and the eighteen block sums read and write whole blocks that the classical recipe never touches. An operation count charges $1$ for an addition and $1$ for a multiplication and charges nothing for moving a number; on a real machine the movement is often the expensive part, and the recursion moves much more data than it saves in arithmetic.
• The error bound is weaker, and it is weaker in a way that matters exactly when it is not expected — see Example (A small entry lost entirely).
• The sizes have to cooperate. The recursion wants even dimensions at every level; an odd size must be padded or handled specially, and the bookkeeping erodes the saving at precisely the moderate sizes where it is thinnest.

The conclusion is not that Theorem 24.13.2 is a curiosity. It is used, at large sizes, with a classical base case and with the accuracy caveat understood. What it is not is a replacement for a definition: $\mathbf{A}\mathbf{B}$ still means what Definition 0.9.9 says it means, and every theorem in this book about products is untouched. A fast algorithm changes the cost of an operation, never the operation.

24.13.7. Summary and transfer

Chapter 21 closed Part V with a complaint and a promise: “every proof here is an existence or a comparison, and not one of them is an algorithm. That is the subject of Chapter 24.” Thirteen sections later, the promise is paid, and what “compute” turned out to mean can be said in one sentence: a computed answer is the exact answer to a nearby question, and how near depends on the algorithm, while how much that nearness costs depends on the problem.

The chapter’s one new assumption was the floating-point model of Definition 24.1.1, a description of a machine stated as a definition. Everything else was proved from it and from the linear algebra of Chapters 0 to 21.

• What a computed answer means. Section 1 stated the model, proved Lemma 24.1.1 and Theorem 24.1.2, and gave the words the rest of the chapter needs: Definition 24.1.3, Definition 24.1.2, and the rule Theorem 24.1.3, that the forward error is at most the condition number times the backward error.
• Solving a system. Section 2 read elimination as a sequence of Schur complements (Theorem 24.2.1), counted it (Proposition 24.2.3), and proved backward error analyses in the cases it could reach: Theorem 24.2.5 outright, and Theorem 24.2.6 for the ordering it displays, the same bound for Chapter 3’s ordering being stated with credit. With it came the growth factor Definition 24.2.1, the bound Theorem 24.2.8 — attained, by Example (A matrix that really grows) — and the pivot-free safety of Cholesky (Theorem 24.2.9).
• Least squares. Section 3 proved that the sensitivity of the least-squares solution involves the residual as well as ${\kappa }_2(\mathbf{A})$ (Theorem 24.3.1), that the normal equations lose twice the digits (Theorem 24.3.2), and that a reflection does not (Theorem 24.3.3) — that last for one reflection applied to a matrix, the accumulation over the whole factorization being stated with credit.
• Eigenvalues. Section 4 built the power method (Theorem 24.4.1), inverse iteration (Theorem 24.4.3), and Rayleigh quotient iteration (Theorem 24.4.4); Section 5 identified the QR algorithm with orthogonal iteration (Theorem 24.5.2), proved its convergence (Theorem 24.5.3), and explained the Hessenberg reduction (Proposition 24.5.4) and the shifts.
• Iterating in a subspace. Section 6 built the Krylov space (Definition 24.6.1), the Arnoldi process (Theorem 24.6.2), and its Hermitian three-term form (Theorem 24.6.3), and showed that the Ritz values (Definition 24.6.2) interlace (Theorem 24.6.4) and carry Chapter 20’s residual bounds. Section 7 minimized over the same space: Theorem 24.7.4, the Chebyshev bound Lemma 24.7.6, and the square root in Theorem 24.7.7.
• One piece of pure mathematics. Section 8 showed that an algorithm and a differential equation can be the same object: the Lax pair Definition 24.8.1, the isospectral theorem Theorem 24.8.1, and the Toda flow’s relation to the QR algorithm (Theorem 24.8.5).
• Applications. Section 9 turned a mass-and-spring system into a generalized eigenproblem (Theorem 24.9.2); Section 10 proved that the best-fitting affine subspace is spanned by the leading singular vectors of the centered data matrix (Theorem 24.10.2) and turned “variance explained” into a theorem (Theorem 24.10.4); Section 11 built the graph Laplacian (Definition 24.11.1), identified its kernel (Theorem 24.11.2), introduced the Fiedler value (Definition 24.11.2) with the minimum it solves (Proposition 24.11.4) and the partition the signs of a Fiedler vector give, and proved the matrix-tree theorem (Theorem 24.11.8).
• Speed. Section 12 factored the Fourier matrix (Theorem 24.12.1) and solved the resulting recurrence (Theorem 24.12.2); this section factored the two-by-two product (Theorem 24.13.2), solved the resulting recurrence (Theorem 24.13.3), proved the one lower bound on $\omega$ this book has (Proposition 24.13.1), and showed what the speed costs (Theorem 24.13.6).

Four moves carried the chapter.

• Conditioning belongs to the problem, stability to the algorithm, and the two multiply. Chapter 16 §08 stated the distinction and could not use it; Theorem 24.1.3 makes it a theorem, and every accuracy claim in the chapter is an instance. Transfer: when an answer is wrong, ask first whether any method could have done better. If the condition number says no, do not change the algorithm — change the problem.
• Measure an answer by a change in the data that makes it exact. A forward error is usually unknowable and a backward error is usually computable, so prove the backward statement and let Theorem 24.1.3 supply the rest (§§01, 02, 03). Transfer: when a quantity you want to bound is out of reach, look for a perturbation of the input that the computed output answers exactly; bounding that is often elementary.
• Every eigenvalue algorithm iterates, because no finite one can exist. Chapter 11 §08 gave the reason — a finite procedure using arithmetic and radicals would solve the general quintic — and Sections 4 to 7 are what is left once that door is shut: a sequence of subspaces, each cheap to extend, each carrying a certified error bar from Chapter 20. Transfer: when exact solution is impossible in principle, the design question becomes what to iterate in; the Krylov space is the answer whenever the only affordable operation is one matrix–vector product.
• A fast algorithm is a factorization. The Fourier matrix factors into sparse pieces and the two-by-two product factors into seven multiplications; in both cases the cost follows by solving a recurrence exactly, and in both cases the underlying object is unchanged (§§12, 13). Transfer: to make an operation fast, do not look for a shortcut in the arithmetic. Look for a factorization of the operation itself, and count the recursive calls against the size reduction.

The chapter paid the following promises, in the words in which they were made.

Some things the chapter did not do, and the chapter’s index collects them in one list. Section 2 proved the backward error of a factorization for the ordering it displays and stated with credit that Chapter 3’s ordering satisfies the same bound; it also showed that the growth bound for partial pivoting is far too large to prove stability, and said so rather than tidying it. Section 3 proved the backward error of a single reflection and stated with credit both the accumulation over a whole factorization and the rates at which classical and modified Gram–Schmidt lose orthogonality. Section 4 stated, and did not prove, that Rayleigh quotient iteration converges at a quadratic rate, not a cubic one, for a general diagonalizable matrix. Section 5 quoted the theorem of Abel and Ruffini, which is why every eigenvalue method here iterates, and proved the first column of the convergence theorem of the QR algorithm while stating the general form with credit. Section 6 stated the classical convergence estimate for the extreme Ritz values and proved the clean special case alone, and described the loss of orthogonality in Lanczos without analyzing it. Section 8 derived the continuous Toda statement only after granting that $t\mapsto \mathbf{Q}(t)$ is differentiable, which is analysis this book has not set up; its discrete half, Theorem 24.8.5, rests on nothing unproved. Section 10 proved every deterministic statement about principal components and said plainly that the inference from a sample to a larger population is statistics, which this book does not develop; Chapter 13 §11 promised that inference, and promised as well an account of how a singular value decomposition is actually computed for a large matrix, which no section here gives — §10 says which route to avoid and why, and points at the orthogonal reductions of §§03 and 05, and there it stops. Both over-promises are recorded rather than hidden. In this section, the recursive form of Strassen’s error bound, Winograd’s lower bound of seven multiplications, the equality of the two definitions of $\omega$, everything known about $\omega$ below ${\log }_{2}7$, and Laderman’s twenty-three-product identity were stated with credit and not proved.

Nothing in this book is proved from any of them. Each is quoted exactly where a reader would otherwise expect a proof, and every theorem in the chapter either proves its own case or carries what it needs as an explicit hypothesis; since all of them 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: Would twenty-three products for three-by-three beat Strassen?) (b), where Laderman’s identity is handed to the reader as a hypothesis and the deduction from it is the reader’s own. Chapter 21’s quoted fact (A7) was not used anywhere in this chapter, and no analysis beyond Chapter 16’s (A1)–(A6) was imported.

The book’s account of computation ends here. Twenty chapters built objects and proved things about them; this one asked what it costs to produce those objects from data of finite accuracy, and found that the answer splits cleanly in two. How much accuracy the data deserve is a question about the problem, answered by a condition number. How much of that the algorithm delivers is a question about the algorithm, answered by a backward error. Everything in this chapter — the pivoting, the reflections, the shifts, the Krylov spaces, the seven products — is an attempt to keep the second small when the first is out of our hands.

24.13.8. Exercises

A. Check your understanding

B. Practice