Strassen's Algorithm
Fast Matrix Multiplication via Divide-and-Conquer with $O(n^{\log_2 7})$ Complexity
Advanced
1. Introduction
The naive product $C = AB$ of two $n \times n$ matrices computes each entry $c_{ij} = \displaystyle\sum_{k=1}^n a_{ik} b_{kj}$ with $n$ multiplications and $n-1$ additions, giving $O(n^3)$ in total.
In 1969, Volker Strassen showed that the product of two $2 \times 2$ matrices can be computed with only seven multiplications instead of eight, and that recursive application reduces the overall complexity to
This was the first proof that $O(n^3)$ is not necessary for matrix multiplication, and it launched the field of fast matrix multiplication.
2. Complexity of naive matrix multiplication
Consider the $2 \times 2$ matrix product:
By the standard definition
this requires 8 multiplications and 4 additions. Splitting an $n \times n$ matrix into $n/2 \times n/2$ blocks and applying this recursively gives the recurrence
and by the Master theorem (the theorem that gives the asymptotic solution of divide-and-conquer recurrences $T(n) = a\,T(n/b) + f(n)$ by comparing $f(n)$ with $n^{\log_b a}$) we obtain $T(n) = O(n^{\log_2 8}) = O(n^3)$.
3. Strassen's idea — seven multiplications
Strassen introduced the following seven auxiliary products $M_1, \ldots, M_7$:
Strassen's seven products
$$\begin{aligned} M_1 &= (a_{11} + a_{22})(b_{11} + b_{22}) \\ M_2 &= (a_{21} + a_{22})\, b_{11} \\ M_3 &= a_{11}\, (b_{12} - b_{22}) \\ M_4 &= a_{22}\, (b_{21} - b_{11}) \\ M_5 &= (a_{11} + a_{12})\, b_{22} \\ M_6 &= (a_{21} - a_{11})(b_{11} + b_{12}) \\ M_7 &= (a_{12} - a_{22})(b_{21} + b_{22}) \end{aligned}$$From these seven $M_i$, every entry of $C$ can be recovered using only additions and subtractions:
Recovery of the entries of $C$
$$\begin{aligned} c_{11} &= M_1 + M_4 - M_5 + M_7 \\ c_{12} &= M_3 + M_5 \\ c_{21} &= M_2 + M_4 \\ c_{22} &= M_1 - M_2 + M_3 + M_6 \end{aligned}$$3.1 Verification
For example, expanding $c_{11}$:
which is the correct result. The other entries can be verified analogously.
3.2 Complexity
Splitting an $n \times n$ matrix into $n/2 \times n/2$ blocks and applying Strassen recursively to compute the $M_i$, the multiplication count drops from 8 to 7 and the recurrence becomes
By the Master theorem,
The number of additions/subtractions grows from 4 to 18 per split — these are additions of block matrices, not scalars, and 18 of them are needed at each recursion level. Since they cost only $O(n^2)$, for large $n$ the saving in multiplications dominates.
3.3 Worked example: Strassen on $2 \times 2$ matrices
Take $A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}$ and $B = \begin{pmatrix} 5 & 6 \\ 7 & 8 \end{pmatrix}$.
Step 1. Compute the seven auxiliary products.
Step 2. Recover the result.
Check: $AB = \begin{pmatrix} 1\cdot5+2\cdot7 & 1\cdot6+2\cdot8 \\ 3\cdot5+4\cdot7 & 3\cdot6+4\cdot8 \end{pmatrix} = \begin{pmatrix} 19 & 22 \\ 43 & 50 \end{pmatrix}$ ✓
4. Recursive structure (diagram)
5. Practical considerations
5.1 Crossover (threshold)
Because Strassen has 18 additions/subtractions, the overhead dominates for small matrices. In practice, when the matrix size at the bottom of the recursion drops below a threshold $n_0$, the algorithm switches to the naive product.
Typically $n_0$ is in the tens to hundreds, depending on hardware (cache size, SIMD width) and memory access patterns.
5.2 Numerical stability
Strassen uses many additions/subtractions, so cancellation errors accumulate. In floating-point arithmetic (following Higham, Accuracy and Stability of Numerical Algorithms, Ch. 23, using the $\infty$-norm and the standard recursion that switches to naive multiplication at a threshold), the forward error of the naive product is bounded by
whereas Strassen satisfies a weaker bound
(where $u$ is the unit round-off). For high-precision applications (matrix inversion, eigenvalue decompositions, etc.) this requires care. Multi-precision integer multiplication (no rounding error) does not have this stability issue.
5.3 Odd sizes
If $n$ is odd, the matrix cannot be split evenly. The standard fix is padding (add one row and one column to make the dimension even); the padded entries are removed at the end.
Example: For $n=3$, zero-pad the $3 \times 3$ matrix to $4 \times 4$.
Run Strassen on $C' = A' B'$ and extract the top-left $3 \times 3$ block to recover $AB$. When recursing, repad whenever the current dimension becomes odd.
5.4 Memory usage
Each level of recursion needs temporary storage for the $M_i$, increasing memory consumption over the naive product. The depth of the recursion is $O(\log n)$ and each level needs $O(n^2)$ temporary memory, so the additional memory is $O(n^2 \log n)$ in total.
5.5 Cache efficiency
The recursive partitioning of Strassen breaks the cache locality that the naive triple loop (ijk or ikj order) achieves so naturally. In particular, the creation and destruction of temporaries for each $M_i$ multiplies cache misses, and memory bandwidth often becomes the bottleneck.
For this reason, high-performance BLAS libraries such as BLIS, Intel MKL, and OpenBLAS do not adopt Strassen; instead they keep the asymptotic $O(n^3)$ complexity but extract peak performance with cache blocking (hierarchical block-by-block organisation) combined with SIMD. Strassen wins on multiplication count, but a straightforward recursive implementation incurs extra memory traffic and temporary storage, so it is not always faster in wall-clock time even for $n$ in the few thousands.
5.6 Implementation tips
- Reuse workspace: instead of allocating new matrices at each recursion level, pre-allocate the total memory once and reuse the same buffers across the seven auxiliary products. This eliminates allocator pressure and GC overhead.
- Avoid copies (strided views): represent sub-matrices as views with a leading dimension (a "stride") rather than copying the data. This matches the BLAS
ldaconvention and lets sub-block operations proceed in place. - SIMD-accelerated naive base case: switch off Strassen below a recursion threshold $n_0$ (values in the range 64–256 are common, but $n_0$ depends strongly on the hardware, the data type and how well the base-case multiply is optimised, so it has to be measured) and fall through to a SIMD-vectorised naive block multiply (for example with $4 \times 4$ register blocks). This base case dominates the runtime in practice.
- Parallelism: once the required block sums and differences have been formed, the seven auxiliary products $M_i$ are independent, so a thread/task pool can dispatch them in parallel at each recursion level. The speed-up actually obtained depends on thread-creation cost, memory bandwidth and workspace pressure, so it is well short of a clean seven-fold gain.
6. Winograd variant
S. Winograd showed a variant that uses the same 7 multiplications as Strassen but reduces additions/subtractions from 18 to 15. The asymptotic complexity is the same $O(n^{2.807})$, but the constant factor is improved, so it is somewhat faster in practice.
The Winograd variant is the standard refinement that keeps the seven multiplications while reducing the additions, and it is the form usually considered when a fast matrix multiplication is implemented. As noted in 5.5, however, the mainstream high-performance BLAS implementations do not use Strassen-style algorithms at all.
7. History of the matrix multiplication exponent $\omega$
Writing the number of operations needed for $n \times n$ matrix multiplication as $O(n^\omega)$, $\omega$ is called the matrix multiplication exponent. The theoretical lower bound is $\omega \geq 2$ (since the output has $n^2$ entries).
| Year | Authors | Upper bound on $\omega$ |
|---|---|---|
| — | Naive method | $3.000$ |
| 1969 | Strassen | $2.807$ |
| 1978 | Pan | $2.796$ |
| 1981 | Schönhage | $2.548$ |
| 1986 | Strassen | $2.479$ |
| 1990 | Coppersmith–Winograd | $2.376$ |
| 2012 | Williams | $2.3727$ |
| 2024 | Duan–Wu–Zhou | $2.371339$ |
Whether $\omega = 2$ is achievable is one of the major open problems in computational complexity theory.
7.1 Beyond Strassen and "galactic algorithms"
Post-Strassen improvements rely on more elaborate tensor-decomposition constructions and the laser method. From Coppersmith–Winograd onward, these algorithms have driven the exponent down to $\omega \approx 2.37$, but the hidden constants and lower-order terms are astronomically large. At any realistic matrix size these algorithms lose not only to Strassen but to naive multiplication as well.
Such algorithms—asymptotically fast but practically slower or outright unimplementable—are known as galactic algorithms. In this taxonomy Strassen is one of the few fast matrix-multiplication algorithms that is not galactic, which is why it is practical enough to have inspired real implementations.
8. Summary
- Strassen showed how to compute the product of two $2 \times 2$ matrices using only 7 multiplications.
- Applied recursively, the product of $n \times n$ matrices runs in $O(n^{2.807})$.
- In practice, a hybrid strategy that switches to the naive product below a threshold size is used.
- Floating-point implementations need care with rounding error. In exact arbitrary-precision integer arithmetic the computation is exact, so this rounding-error issue does not arise.
- The Winograd variant reduces additions/subtractions from 18 to 15.
- Minimizing the matrix multiplication exponent $\omega$ is an active area of research in computational complexity.
References
- V. Strassen, "Gaussian elimination is not optimal," Numerische Mathematik, vol. 13, pp. 354–356, 1969. — DOI: 10.1007/BF02165411
- T. H. Cormen, C. E. Leiserson, R. L. Rivest, C. Stein, Introduction to Algorithms, 4th ed., MIT Press, 2022. — Chapter 4.2: Strassen's algorithm for matrix multiplication.
- N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002. — Chapter 23: Fast matrix multiplication.
- Strassen algorithm — Wikipedia
- Matrix multiplication algorithm — Wikipedia
Related articles
- Matrix product — definition and basic properties of matrix multiplication
- LU decomposition
- Fast Fourier transform — another application of divide-and-conquer
Frequently Asked Questions
Q1. What is Strassen's algorithm?
It is a divide-and-conquer algorithm that speeds up the multiplication of $n \times n$ matrices from the standard $O(n^3)$ to $O(n^{2.807})$. It computes the product of two $2 \times 2$ block matrices with 7 multiplications and 18 block additions/subtractions instead of the usual 8 multiplications. It was published by Volker Strassen in 1969.
Q2. What is the practical crossover threshold?
The threshold depends strongly on the implementation, the hardware, the data type and how well the base-case multiply is optimised, so it has to be measured. Values in the range of a few tens to a few hundred are common. Note that the mainstream high-performance BLAS implementations (BLIS, Intel MKL, OpenBLAS) do not use Strassen; they speed up the naive product with cache blocking and SIMD instead (see 5.5).
Q3. What are the drawbacks of Strassen's algorithm?
The main drawbacks are: (1) larger rounding error in floating-point arithmetic, (2) overhead from the extra additions, which dominates for small matrices, (3) extra memory traffic and temporary storage in a straightforward recursive implementation, and (4) the need to pad odd-sized matrices.