The standard algorithm for multiplying two $n\times n$ matrices requires $n^3$ multiplications, but, remarkably, in 1969, Strassen showed it could be done in $\approx n^{2.8}$ multiplications. This is fascinating to me, but the algorithm itself is arcane and seemingly comes out of the blue. In this little post I hope to demystify it by showing how you could have discovered Strassen’s algorithm using a computer-algebra system.

The product of $2\times 2$ matricies

is typically written

but this requires 8 multiplications and multiplication is computationally expensive. Can we do better?

Remarkably, yes. In 1969 Volker Strassen showed it can be done with 7:

While ostensibly about $2\times 2$ matrices, this technique can be used to multiply arbitrary $n\times n$ matricies as well. The basic idea is to break the $n\times n$ matricies into four $n/2\times n/2$ quadrant matrices and use Strassen's algorithm to compute the $n\times n$ product using seven matrix-multiplies over the $n/2\times n/2$ quadrants. Letting $t(n)$ denote the number of multiplies required to multiply $n\times n$ matricies using Strassen's algorithm, this gives the following recurrence relation for the value of $t(n)$:

When $n$ is not a power of two, we can pad the matrices with zeros to make the recurrence terminating. Doing so gives $t(n)=7^{\log_2n}=n^{\log_27}\approx n^{2.8}$.

I find this remarkable: It wouldn't have even occurred to me the typical $n^3$ algorithm for matrix multiplication is sub-optimal. Moreover, the expressions for $P_i$ seemingly comes from thin air. How would one go about "discovering" Strassen's algorithm?

In 1970 Richard Brent gave an algebraic method (p. 33) which we will work with for the remainder of this essay. At a high level, we will set up an algebra problem, the solutions to which correspond to ways to multiply $2\times 2$ matrices using 7 multiplications. This will have 64 equations with 84 variables, so to avoid solving this by hand we will use a SMT solver.

First, let $C=AB$ and to simplify indexing let $\bar{a}=(a_{11},a_{12},a_{21},a_{22})$ be $A$ in row-major form and let $\bar{b}$ and $\bar{c}$ be likewise. Because $a_{ij} a_{kl}$ nor $b_{ij} b_{kl}$ appear in the expression for $C$, the $k$th product will be of the form

where $\alpha^{(k)}_i,\beta^{(k)}_i\in\{-1,0,1\}$ are variables indicating the leading coefficent of $\bar{a}_i$/$\bar{b}_j$ in the $k$th product. For example, in Strassen's algorithm $P_3 = a_{11} (b_{12}-b_{22})$. As $\bar{a}_1$ appears in $P_3$, $\alpha^{(3)}_1=1$. As $-\bar{b}_4$ appears in $P_3$, $\beta^{(3)}_4=-1$.

$\bar{c}_l$ will be a sum of $P_k$s, so letting $\gamma_l^{(k)}\in\{-1,0,1\}$ denote $P_k$'s coefficient in the expression for $\bar{c}_l$,

To make this an algebra problem, we need another formula for $\bar{c}_l$. We can use the standard expression for $C=AB$

and let $T_{ijl}$ be $1$ if $\bar{a}_i\bar{b}_j$ appears in the above expression for $\bar{c}_l$ and $0$ otherwise, giving the desired formula for $\bar{c}_l$ below.

To equate these two expressions, note

so we have

an algebra problem, the solutions to which correspond to ways to compute $C=AB$ using $k$ multiplies. Though, as hinted before, because our matrices have 4 entries each, the $i$, $j$, and $l$ indicies range over 4 values, making $4\cdot 4\cdot 4=64$ equations. Adding to our complications, to recover Strassen's algorithm we'd set $k=7$, meaning there are $3\cdot 4\cdot 7=84$ $w_l^{(k)}$, $\alpha^{(k)}_i$, and $\beta^{(k)}_j$ variables. We will need a computer to solve this algebra problem.

I used Z3 for this, an SMT solver written at Microsoft. The code is straightforward, you can read it if you like, but trying to solve the system of equations above will fail at first blush. To narrow the search space, I added two constraints to the equations:

1. Swapping $P_i$ with $P_j$ yields a new, but uninteresting, solution. To stop the solver from finding these duplicates, I require the $P_k$s to be lexographically ordered, viewing them as strings of $\alpha^{(k)}_i$ and $\beta^{(k)}_j$ variables.
2. Negating $P_k$ and all its corresponding $w^{(k)}_l$ values, or negating both of $P_k$'s summations yields a new, but uninteresting, solution. To stop the solver from finding these duplicates, I require the first non-zero $\alpha^{(k)}_i$ and $\beta^{(k)}_j$ values in $P_k$'s summations to be positive. This picks a sign for $P_k$, as exactly one of $(-1)(\Sigma\dots)$ and $(\Sigma\dots)$ has a positive first non-zero term.

With these symmetries removed, the solver finds several Strassen-like schemes in a minute. Nice - we've rediscovered Strassen's algorithm from first principles.

Thanks to Pete, Jacob, Groceries, and Elam for interesting conversations, and Entrepreneurs First for for hosting the weekend hackathon where I got to weave this together.