Exploiting the Structure in Tensor Decompositions for Matrix Multiplication

Abstract: We present a new algorithm for fast matrix multiplication using tensor decompositions which have special features. Thanks to these features we obtain exponents lower than what the rank of the tensor decomposition suggests. In particular for $6\times 6$ matrix multiplication we reduce the exponent of the recent algorithm by Moosbauer and Poole from $2.8075$ to $2.8016$, while retaining a reasonable leading coefficient.

• The paper develops a recursive method that groups rank-one terms forming smaller matrix multiplication tensors, reducing exponents beyond rank-based Strassen recursion while preserving exact computation.
• The method achieves exponents of 2.8016 for 6×6 multiplication and 2.836 for 3×3 multiplication, with the 6×6 scheme outperforming standard algorithms at sufficiently large input sizes.
• The authors combine flip-graph search, sparsification, constrained Hensel lifting, and common-subexpression elimination to find structured decompositions, while identifying memory costs, crossover sizes, and asymptotic overhead as open challenges.

Overview

This paper by Kauers, Moosbauer, and Wood (2602.11041) develops a recursive matrix multiplication algorithm that exploits non-rank-one structure inside tensor decompositions of matrix multiplication tensors. The central observation is that a decomposition of $\langle n,m,p\rangle$ may contain, among its rank-one summands, groups of terms that jointly form smaller matrix multiplication tensors such as $\langle 1,1,2\rangle$ or $\langle 2,2,2\rangle$. Rather than computing these sub-products independently with the standard algorithm, the proposed method treats each group as a single larger multiplication and recurses on it. This yields exponents strictly below what the rank of the decomposition alone would suggest via Strassen's recursion, while keeping the number of base-field multiplications unchanged and retaining a moderate leading coefficient.

The headline results are an exponent of 2.8016 for $6\times 6$ matrix multiplication, improving on the recent Moosbauer–Poole algorithm's $2.8075$, and an exponent of 2.836 for $3\times 3$ multiplication, improving Laderman's classical $2.854$. The authors report that their $6\times 6$-based algorithm outperforms the standard algorithm in total operation count for $n \geq 1000$, and simulations indicate consistent outperformance over Winograd–Strassen from around $10^{10}$.

Background and motivation

The paper works in the standard tensor framework: a restriction $\langle 1,1,2\rangle$0 corresponds to a rank-$\langle 1,1,2\rangle$1 decomposition and, combined with Strassen's block recursion, gives exponent $\langle 1,1,2\rangle$2 (2602.11041). The leading coefficient depends on the number of additions/scalar multiplications $\langle 1,1,2\rangle$3 in the restriction; the authors restate the well-known bound for this coefficient, which recovers $\langle 1,1,2\rangle$4 for Strassen's original scheme ($\langle 1,1,2\rangle$5).

All improvements to the upper bound on $\langle 1,1,2\rangle$6 since 1981 rely on Schönhage's Asymptotic Sum Inequality (ASI) or its generalizations, but these are existence results whose approximate algorithms are impractical due to precision requirements and enormous leading coefficients. The present work instead targets exact algorithms in the practical regime, following the line of research based on flip graph search, numerical optimization, and reinforcement learning. The key precedent is Schönhage's own observation that after three iterations of Pan's algorithm, eight recursive calls form a $\langle 1,1,2\rangle$7 matrix multiplication that can be computed with Strassen's scheme; Schwartz and Zwecher recently applied a similar idea to Pan's trilinear aggregation algorithms. This paper turns that idea into a fully explicit, implementable recursive algorithm achieving the complexity promised by the generalized ASI (the Cohn–Umans variant), which had previously only been used for asymptotic bounds.

The algorithm

Given a restriction of the form

$\langle 1,1,2\rangle$8

where the right-hand side contains copies of non-rank-one tensors, Algorithm 1 computes the linear combinations prescribed by the decomposition, but whenever several recursive calls share an input factor or combine into a common output position, it merges them into one larger matrix product and recurses. Rectangular subproblems eventually bottom out at shapes like $\langle 1,1,2\rangle$9, where the standard algorithm is used; alternating between cyclic permutations of the restriction across recursive levels keeps most recursion paths alive longer.

Two theorems establish the complexity. For square schemes, Theorem 3 proves that Algorithm 1 runs in $\langle 2,2,2\rangle$0 operations where $\langle 2,2,2\rangle$1, with an explicit induction giving a leading-coefficient bound analogous to the classical analysis. Theorem 4 then handles rectangular schemes: taking tensor powers of a restriction together with its cyclic permutations yields a decreasing sequence of exponents $\langle 2,2,2\rangle$2 converging to

$\langle 2,2,2\rangle$3

which is exactly the generalized ASI threshold. Moreover, the convergence is quantified: $\langle 2,2,2\rangle$4, so the gap closes as $\langle 2,2,2\rangle$5. Combining this with Lemma 5's bound $\langle 2,2,2\rangle$6 gives an overall cost of $\langle 2,2,2\rangle$7 for a suitable choice of $\langle 2,2,2\rangle$8. The authors note this bound is likely far from tight — they expect the $\langle 2,2,2\rangle$9 factor could be removed with more careful analysis — and, following prior literature, they adopt the idealized leading coefficient (assuming exact recursion down to $6\times 6$0 blocks) for comparisons, while conceding that unlike standard Strassen-like algorithms, no actual input realizes this assumption.

Applied to Moosbauer–Poole's decomposition $6\times 6$1, the method already gives exponent $6\times 6$2; the authors' improved decomposition with 18 copies of $6\times 6$3 reaches $6\times 6$4 with leading coefficient $6\times 6$5 (versus $6\times 6$6 for Strassen). Notably, the same treatment of the $6\times 6$7 case yields only $6\times 6$8, slightly worse than Strassen's exponent — the improvement is not uniform across formats.

Finding structured decompositions

The decompositions are obtained through a five-step pipeline:

1. Flip graph search on a tensor with a copy of $6\times 6$9 (or $2.8075$0) subtracted, minimizing rank; adding the pattern back yields a decomposition containing it.
2. Structure maximization: a second flip graph search maximizes the number of non-overlapping copies of $2.8075$1 patterns, which are easy to detect as components sharing a common factor.
3. Sparsification via random elements of de Groote's symmetry group, chosen because the structural property is invariant under these symmetries (unlike flip moves).
4. Hensel lifting from $2.8075$2 to $2.8075$3, with additional linear constraints appended to the underdetermined lifting system to preserve both the structure and sparsity. The authors emphasize this preservation is not automatic — e.g., two copies of $2.8075$4 in the Moosbauer–Poole decomposition collapse to a single $2.8075$5 modulo 2 — and report that the constrained lifting often produced coefficients in $2.8075$6.
5. Addition optimization using Mårtensson–Wagner common subexpression elimination.

Step 1 dominates the computation time. The resulting decompositions, summarized across formats from $2.8075$7 up to $2.8075$8, show consistent exponent drops relative to the best known rank-based bounds (e.g., $2.8075$9: from $3\times 3$0 to $3\times 3$1; $3\times 3$2: from $3\times 3$3 to $3\times 3$4), though only the $3\times 3$5 case beats Strassen's $3\times 3$6 outright. All decompositions are released electronically.

Simulation

Because the analytic leading coefficients do not reflect actual operation counts, the authors simulate the recursion for concrete input sizes, switching to Winograd–Strassen below $3\times 3$7 and to the standard algorithm below $3\times 3$8. The simulated counts are lower than the theoretical bounds precisely because of these hybrid switches. Two caveats are stated plainly: zero padding makes the operation count sensitive to input size (producing non-smooth curves), and the simulation ignores memory access costs. Improvements over Strassen appear around $3\times 3$9, with consistent outperformance from $2.854$0 — considerably larger than the $2.854$1 crossover suggested by the pure operation-count comparison against the standard algorithm.

Limitations and open questions

The paper is candid about several limitations. The claimed practical relevance rests on operation counts rather than measured runtimes; memory behavior is unexamined. The rigorous complexity guarantee carries an $2.854$2 overhead term that the authors believe removable but have not removed, and the idealized leading coefficient used for comparison is not attained by any actual input. The improvements depend on finding decompositions with favorable structure, and the search succeeded only sporadically: most formats studied do not beat Strassen's exponent, and the $2.854$3 result is worse. Whether the flip-graph-based structure maximization scales to larger formats, whether the $2.854$4 term can be eliminated, and whether heuristic treatment of rectangular subproblems (beyond periodic splitting into squares) can improve constants all remain open.

Conclusion

The paper converts a theoretical device implicit in Schönhage's and Cohn–Umans' asymptotic arguments into an explicit recursive algorithm whose complexity matches the generalized ASI while remaining exact and implementable. By searching for tensor decompositions containing embedded sub-multiplications and recursing on those structures, it achieves concrete exponent improvements for small base cases — most prominently $2.854$5 for $2.854$6 and $2.854$7 for $2.854$8 — with modest leading coefficients. The contribution is thus a bridge between the asymptotic theory of bilinear complexity and the engineering-oriented search for practical fast multiplication algorithms, with its practical impact still contingent on implementation-level considerations the paper does not resolve.