How far are Strassen-like algorithms from being practical?

Mathematically, multiplying two dense $n\times n$ matrices $A$ and $B$ to obtain $C=AB$ using the standard algorithm requires $\Theta(n^3)$ arithmetic operations. More specifically, it performs $n^3$ scalar multiplications and $n^3-n^2$ scalar additions. In practice, implementations of BLAS, particularly the GEMM (general matrix-matrix multiplication) routines, are the de facto standard for dense matrix multiplication. They are heavily optimized for specific hardware architectures and are widely used by higher-level numerical libraries. For example, NumPy can be built against optimized BLAS implementations such as OpenBLAS, BLIS, MKL, or Apple’s Accelerate framework, depending on the platform and build configuration.

I tested the NumPy expression A @ B on my M4 MacBook Pro using power-of-two matrix sizes and observed an empirical scaling between approximately $\Theta(n^{2.6})$ and $\Theta(n^{2.7})$ for $n$ between 512 and 16,384. This is, of course, a measurement of a particular implementation on a particular machine—not the asymptotic complexity of matrix multiplication.

Strassen’s algorithm. In theoretical computer science, there has been a long-standing effort to reduce the asymptotic complexity of matrix multiplication. Strassen’s algorithm is the classic example: it reduces the exponent from $3$ to $\log_2 7 \approx 2.8074$. It is an elegant divide-and-conquer algorithm. Assume for now that $n$ is even. We partition $A$ and $B$ into $2\times2$ block matrices with equally sized blocks:

\[\begin{equation}\label{eqn:partition} A = \begin{bmatrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{bmatrix}, \qquad B = \begin{bmatrix} B_{11} & B_{12}\\ B_{21} & B_{22} \end{bmatrix}. \end{equation}\]

The product $AB$ can then be computed as

\[AB = \begin{bmatrix} M_1 + M_4 - M_5 + M_7 & M_3 + M_5\\ M_2 + M_4 & M_1 - M_2 + M_3 + M_6 \end{bmatrix},\]

where

\[\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}\]

Thus, Strassen replaces the eight block multiplications of the conventional algorithm with seven block multiplications, at the cost of additional additions and subtractions. The original formulation uses 18 block additions/subtractions. If $n$ is a power of 2 and $T(n)$ denotes the number of arithmetic operations, we obtain

\[T(n)=7 \, T\left(\frac n2\right)+18\left(\frac n2\right)^2, \quad\text{with}\quad T(1)=1.\]

Solving this recurrence gives $T(n)=7n^{\log_2 7}-6n^2$, which is where the exponent $\log_2 7$ comes from. Winograd later developed a variant that reduces the number of additions/subtractions from 18 to 15 while retaining seven multiplications. This changes the constant factor but not the asymptotic exponent.

Strassen-like algorithms. Because reducing the number of multiplications is central to improving asymptotic complexity, there has been a long history of searching for algorithms that use fewer multiplications for small matrix sizes. These algorithms need not themselves be recursive. Instead, they provide scalar-level bilinear formulas that can subsequently be used as building blocks for recursive algorithms.

For an $a\times b$ matrix $A$ and a $b\times c$ matrix $B$, the goal is to compute $AB$ using fewer than $abc$ scalar multiplications. A particularly interesting recent development is AlphaTensor, a deep-reinforcement-learning method for discovering matrix multiplication algorithms through tensor decomposition (see the appendix at the end). AlphaTensor has found improved algorithms for a number of small matrix shapes. Importantly, the arithmetic domain matters: some of its most striking results are over the finite field $\mathbb Z_2$, while others apply to standard arithmetic over $\mathbb R$.

For example, AlphaTensor found a 47-multiplication algorithm for $4\times4$ matrix multiplication over $\mathbb Z_2$, improving on the 49 multiplications obtained by applying Strassen twice. Over standard arithmetic, AlphaTensor also found a 76-multiplication algorithm for multiplying a $4\times5$ matrix by a $5\times5$ matrix, improving on the previous bound of 80.

---

Why isn’t Strassen standard in BLAS libraries?

While Strassen’s algorithm is theoretically advantageous over the conventional $\Theta(n^3)$ algorithm, there are several reasons why it does not automatically outperform highly optimized GEMM implementations.

The first issue is data movement. Strassen reduces one matrix multiplication at the cost of additional matrix additions and subtractions. Matrix multiplication has relatively high arithmetic intensity: a large amount of computation is performed on the data loaded from memory. By contrast, matrix addition performs only a small amount of computation per element and is therefore much more sensitive to memory bandwidth. Consequently, the additional additions and subtractions can introduce substantial memory traffic. In a naive implementation that explicitly materializes the intermediate sums, this can require moving considerably more data between levels of the memory hierarchy.

Numerical accuracy is another consideration. Strassen’s method changes the pattern of floating-point operations and can therefore lead to greater rounding-error growth than conventional matrix multiplication. In particular, it forms linear combinations of matrix elements before multiplication. When the elements of $A$ or $B$ have widely differing magnitudes, these additions and subtractions can introduce cancellation and allow rounding errors in the intermediate quantities to affect the final result. This does not make Strassen inherently unusable or categorically unstable, but it means that numerical accuracy must be considered when selecting the algorithm. This consideration is particularly important for a general-purpose numerical library, where users may expect well-characterized and predictable numerical behavior.

---

How to implement Strassen’s (or Strassen-like) algorithm? Idea 1

Since recursion is not necessary, a natural experiment is to apply Strassen only once. Partition $A$ and $B$ into $2\times2$ blocks, as in \eqref{eqn:partition}. We then perform seven multiplications of $\frac{n}{2} \times \frac{n}{2}$ matrices and 18 block additions/subtractions. A conventional multiplication of two $\frac{n}{2} \times \frac{n}{2}$ matrices requires $\frac{n^3}{8}$ scalar multiplications and $\frac{n^3}{8}-\frac{n^2}{4}$ scalar additions. Therefore, the one-level Strassen scheme requires $\frac{7}{8} n^3$ scalar multiplications and $7(\frac{n^3}{8} - \frac{n^2}{4}) + 18\frac{n^2}{4} = \frac{7}{8}n^3 + \frac{11}{4}n^2$ scalar additions. If we simply count every multiplication and addition as one arithmetic operation, the conventional algorithm uses $2n^3-n^2$ operations, whereas the one-level Strassen scheme uses $\frac{7}{4} n^3 + \frac{11}{4} n^2$. The latter is smaller when $n>15$.

I tested this one-level scheme with NumPy on my M4 MacBook Pro, with several optimizations to reduce temporary matrices and streamline the additions. In my experiments, it did not outperform A @ B for $n<2^{14}$. At $n=2^{14}$, the two timings were very close, with the difference smaller than the observed standard deviation. For $n=2^{15}$, the trend suggested that the one-level scheme might begin to win, but my computer did not have enough memory to complete the experiment.

One may therefore wonder why A @ B can be so fast even though it performs more arithmetic operations. The answer is that arithmetic-operation count is only part of the story. Highly optimized GEMM implementations are designed around the memory hierarchy, vector instructions, cache blocking, register tiling, packing, and hardware-specific micro-kernels. They can achieve extremely high arithmetic throughput. By contrast, the extra additions introduced by Strassen have much lower arithmetic intensity. They must read and write large blocks of data while performing relatively little computation. The resulting memory traffic can dominate the runtime.

---

How to implement Strassen’s (or Strassen-like) algorithm? Idea 2

If partitioning $A$ and $B$ into four large blocks does not fit the hardware well, another possibility is to partition them into many small blocks—for example, using $a\times b$ blocks for $A$ and $b\times c$ blocks for $B$, where $a$, $b$, and $c$ are small. This organization is much closer to the structure of a high-performance GEMM implementation, in which computation is organized around cache blocks, micro-tiles, and register-level kernels. An illustrative discussion of this organization can be found in UT Austin’s LAFF-On Programming for High Performance course.

The idea is to search for Strassen-like formulas whose operands match the dimensions of the small tiles (green and blue) handled by the innermost GEMM micro-kernel. For example, one could imagine finding a Strassen-like algorithm for a small tile shape that reduces the number of scalar multiplications while retaining efficient use of vector instructions, registers, and the L1 cache.

This direction is closely related to the motivation behind AlphaTensor’s hardware-tailored search. In its original work, AlphaTensor was not limited to minimizing the number of multiplications: the authors also modified the reward to optimize actual runtime on a target GPU or TPU. In one experiment, they searched for a $4\times4$ block algorithm for multiplying large matrices and optimized its runtime on specific hardware.

The hardware constraints are important. The useful tile dimensions are determined by factors such as the register file and vector width. They are therefore highly architecture-dependent.

There is another important limitation, too. If a Strassen-like formula is used only as a fixed-size micro-kernel, it reduces the constant factor but does not improve the asymptotic complexity. The overall algorithm remains $\Theta(n^3)$. To obtain a genuinely sub-cubic algorithm, the smaller multiplication scheme must itself be applied recursively.

---

Concluding remarks

Strassen’s algorithm is one of the most beautiful examples of the gap between asymptotic complexity and practical performance.

On paper, replacing eight multiplications by seven seems like an obvious win. In a modern computer, however, the cost of an operation is determined not only by the arithmetic itself, but also by data movement, cache locality, vectorization, register pressure, memory bandwidth, parallelism, and the highly optimized structure of existing GEMM implementations.

This makes fast matrix multiplication an unusually interesting problem for algorithm discovery. The objective need not be simply to minimize the number of scalar multiplications. A truly practical algorithm may need to optimize a much richer objective, including hardware-specific runtime.

This is precisely where AI-based algorithm discovery becomes particularly intriguing. AlphaTensor demonstrated that reinforcement learning can discover non-obvious matrix multiplication algorithms by searching over tensor decompositions, and it also showed that the search objective can be adapted to practical runtime on specific hardware.

The theoretical frontier has continued to move as well: the best published upper bound on the asymptotic matrix-multiplication exponent has been pushed below $2.3712$ in recent work. Yet the gap between such asymptotic results and a faster replacement for production GEMM remains enormous.

So, how far are Strassen-like algorithms from being practical? The answer is: perhaps not as far as conventional wisdom suggests—but closing the gap requires optimizing for the machine, not just for the algebra. I remain optimistic that AI-driven algorithm discovery, combined with a detailed understanding of modern computer architectures, may provide a path toward bridging that gap.

---

Appendix: Tensor decomposition and the discovery of Strassen-like algorithms

At a high level, Strassen’s algorithm expresses each element of $C$ as a sum of products between linear combinations of entries of $A$ and linear combinations of entries of $B$. The standard formula $C[i,j]=\sum_{p=1}^k A[i,p]B[p,j]$ can already be viewed in this form. It is a sum of $k$ products, where the $p$-th product uses the linear combinations \(\sum_{p'=1}^k \mathbf{1}_{p'=p}A[i,p']\) and \(\sum_{q'=1}^k \mathbf{1}_{q'=p}B[q',j]\). Strassen’s algorithm uses a different set of linear combinations. For example, $C[1,1]$ is computed as $C[1,1]=M_1+M_4-M_5+M_7$, where

\[\begin{aligned} M_1 &= (A[1,1]+0+0+A[2,2]) (B[1,1]+0+0+B[2,2]),\\ M_4 &= (0+0+0+A[2,2]) (-B[1,1]+0+B[2,1]+0),\\ -M_5 &= (A[1,1]+A[1,2]+0+0) (0+0+0-B[2,2]),\\ M_7 &= (0+A[1,2]+0-A[2,2]) (0+0+B[2,1]+B[2,2]). \end{aligned}\]

This representation makes the connection to tensor decomposition explicit.

The matrix multiplication tensor. We can extend this formulation to arbitrary dimensions. Let $A\in\mathbb{R}^{m\times k}$, $B\in\mathbb{R}^{k\times n}$, and $C=AB\in\mathbb{R}^{m\times n}$. Define a three-dimensional tensor $\mathcal T$ of dimensions $(mk)\times(kn)\times(mn)$. Its entries can be indexed as $\mathcal{T}[(i’,p’),(q’,j’),(i,j)]$. The tensor encodes the entire matrix multiplication operation:

\[\begin{equation}\label{eqn:T.def} C[i,j] = \sum_{i'=1}^m \sum_{p'=1}^k \sum_{q'=1}^k \sum_{j'=1}^n \mathcal{T}[(i',p'),(q',j'),(i,j)] A[i',p']B[q',j']. \end{equation}\]

The tensor entry is one exactly when $i’=i$, $p’=q’$, and $j’=j$; and is zero otherwise. Thus, the expression reduces to the familiar formula $C[i,j]=\sum_{p=1}^k A[i,p]B[p,j]$. This tensor is fixed by the dimensions $m,k,n$; it does not depend on the particular numerical values of $A$ and $B$. This is the matrix multiplication tensor used in modern work on fast matrix multiplication and algorithm discovery, such as AlphaTensor.

Tensor decomposition. Now consider three collections of vectors, \(\{ \mathbf{u}^r \mid r = 1,\ldots,R \}\), \(\{ \mathbf{v}^r \mid r = 1,\ldots,R \}\), and \(\{ \mathbf{w}^r \mid r = 1,\ldots,R \}\), where the vectors have lengths $mk$, $kn$, and $mn$, respectively. Suppose that $\mathcal T$ can be decomposed as

\[\begin{equation}\label{eqn:T.decomposition} \mathcal T = \sum_{r=1}^R \mathbf u^r\otimes \mathbf v^r\otimes \mathbf w^r. \end{equation}\]

Element-wise, this means

\[\begin{equation}\label{eqn:T.decomposition2} \mathcal{T}[(i',p'),(q',j'),(i,j)] = \sum_{r=1}^R \mathbf u^r[(i',p')] \mathbf v^r[(q',j')] \mathbf w^r[(i,j)]. \end{equation}\]

Substituting \eqref{eqn:T.decomposition2} to \eqref{eqn:T.def} gives

\[C[i,j] = \sum_{r=1}^R \mathbf{w}^r[(i,j)] \left( \sum_{i'=1}^m\sum_{p'=1}^k \mathbf{u}^r[(i',p')] A[i',p'] \right) \left( \sum_{q'=1}^k\sum_{j'=1}^n \mathbf{v}^r[(q',j')] B[q',j'] \right).\]

The significance of this expression is that the two quantities inside parentheses are independent of $(i,j)$. They can therefore be computed once for each $r$, multiplied together, and then distributed to the output according to $\mathbf w^r$.

Thus, a rank-$R$ tensor decomposition \eqref{eqn:T.decomposition} gives a matrix multiplication algorithm using $R$ scalar multiplications, together with the necessary linear combinations and additions. This is the fundamental connection between tensor rank and fast matrix multiplication. For a fixed set of dimensions $m,k,n$, the matrix multiplication tensor is uniquely determined. The algorithm-discovery problem can therefore be viewed as searching for a low-rank decomposition of this tensor. Strassen’s algorithm corresponds to the case $m=k=n=2$ and $R=7$.