Source-linked AI summary
More Asymmetry Yields Faster Matrix Multiplication
Josh Alman, Ran Duan, Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, Renfei Zhou
TL;DR
Fast matrix multiplication algorithms face limitations in the laser method and related approaches. This paper develops a more asymmetric laser-method improvement, yielding improved bounds for square and rectangular matrix multiplication and a better bound for several related problems.
Problem
The laser method and broader approaches to matrix multiplication have known limitations, including limitations preventing them from proving ω = 2.
Method
The paper improves the laser method for powers of CW5 by using more asymmetry, including shared level-k blocks in the Y- and Z-dimensions while preserving independence of level-1 blocks.
Results
The method yields improved bounds on ω and several rectangular matrix multiplication exponents, while improving µ from µ < 0.527661 to µ < 0.5275 for several related problems.
Takeaways & Limitations
The improved µ bound gives O(n2.5275)-time algorithms for APSP in unweighted directed graphs, minimum witnesses of Boolean Matrix Multiplication, and All-Pairs Bottleneck Paths in node-weighted graphs.
Takeaways & Limitations
The approach still requires solving a larger, considerably harder non-convex optimization problem and cannot currently allow level-ℓ X-blocks to be shared.
Abstract
from arXiv · showhide
We present a new improvement on the laser method for designing fast matrix multiplication algorithms. The new method further develops the recent advances by [Duan, Wu, Zhou FOCS 2023] and [Vassilevska Williams, Xu, Xu, Zhou SODA 2024]. Surprisingly the new improvement is achieved by incorporating more asymmetry in the analysis, circumventing a fundamental tool of prior work that requires two of the three dimensions to be treated identically. Together with recent advances in numerical optimization [Dupont et al. arXiv:2608.16884], the method yields a new bound on the square matrix multiplication exponent $ω<2.371177$, improved from the previous bound of $ω<2.371552$. We also improve the bounds of the exponents for multiplying rectangular matrices of various shapes.
1 Introduction
Matrix multiplication remains a central complexity problem, with recent work pushing the exponent ω toward known methodological limits. This paper introduces a laser-method improvement that yields improved square and rectangular matrix multiplication bounds, while requiring substantially harder numerical optimization.
- Results: ω < 2.371552 was the previous best bound, obtained by optimizing the recent approach of Duan, Wu, and Zhou.The paper’s new bound is further improved in later work using larger-scale numerical optimization.
- Motivation: ω has a lower bound of 2, while the laser method on CW5 cannot prove ω < 2.3078.Recent matrix multiplication research has focused on narrowing the gap between the best upper bounds and this limitation.
- Contribution: The paper’s main contribution is a new improvement to the laser method applied to powers of the Coppersmith-Winograd tensor CW5.The method builds on the approaches of Duan, Wu, and Zhou and Vassilevska Williams, Xu, Xu, and Zhou.
- Results: The method also improves rectangular exponents ω(1, k, 1), defined for multiplying n × n^k by n^k × n matrices.Table 1 compares the new bounds with those of Vassilevska Williams, Xu, Xu, and Zhou; the k = 1 entry uses parameters found in later work, while other entries analyze the fourth tensor power.
- Limitations: The resulting optimization problems are larger and more complex, and solving a fixed ω(1, k, 1) instance can take many days.The authors also report that they could not solve the constraint program for Coppersmith’s parameter α.
2 Technical Overview
The paper develops a more asymmetric laser-method analysis for tensor powers of Coppersmith–Winograd tensors, allowing controlled sharing of variable blocks while retaining many matrix multiplication tensors. This approach improves prior bounds but introduces larger optimization problems and leaves structural limitations.
- The Laser Method: The laser method converts Kronecker powers of tensors into large direct sums of matrix multiplication tensors, whose rank can be bounded using Schönhage’s asymptotic sum inequality.The construction starts from CWq, takes a large Kronecker power, applies zeroing out, and then uses the resulting direct sum to obtain an algorithmic bound.
- Recursive Applications: Recursive laser applications become harder at higher levels because each level-ℓ subtensor is a sum of many level-(ℓ−1) subtensors.The increasing complexity creates algorithmic and computational challenges, and prior work also addressed penalties caused by non-unique marginal distributions.
- More Asymmetry: The new approach requires level-ℓ subtensors to avoid sharing level-ℓ X-blocks, while allowing them to share higher-level Y- and Z-blocks when their level-1 blocks remain independent.This relaxes the prior disjointness requirements and allows more level-ℓ subtensors to survive the zeroing process.
- More Asymmetry: Different zeroing-out procedures for X, Y, and Z introduce asymmetry that circumvents a major limitation of more symmetric prior analyses.Although the final goal treats the dimensions symmetrically, the analysis uses dimension-specific procedures to preserve independence at the level-1 blocks.
- Bounds and Optimization: The new analysis improves the square-matrix bound from ω < 2.374399 to ω < 2.37432 for the analyzed CWq power.The final parameter selection depends on a large, non-convex optimization problem, whose expanded formulation is harder to solve because the X- and Y-dimensions are no longer symmetric.
- Limitations: The technique still requires unique level-ℓ X-block ownership and is therefore described as reaching the limit of generalizing the techniques introduced by Duan, Wu, and Zhou.The authors also report that solving the resulting nonlinear programs can take many days and that the relevant Coppersmith exponent could not be computed.
3 Preliminaries
The preliminaries define tensor operations, rank and asymptotic rank, matrix multiplication tensors, and the exponent functions used to analyze square and rectangular multiplication. They then introduce the Coppersmith–Winograd tensor and the partitioning underlying the laser method.
- Tensor preliminaries: Tensor rank is the minimum number of rank-1 tensors whose sum equals a tensor, and asymptotic rank captures its limiting growth under tensor powers.Restrictions and zero-outs preserve the relevant rank upper bounds, supporting the laser method’s transformations.
- Matrix multiplication exponents: The matrix multiplication exponent ω is defined from ranks of square matrix multiplication tensors, while ω(a, b, c) measures arithmetic complexity for rectangular dimensions.The paper focuses on bounds for ω(1, κ, 1).
- Asymptotic sum inequality: Schönhage’s asymptotic sum inequality converts asymptotic-rank bounds for direct sums of matrix multiplication tensors into upper bounds on ω and rectangular exponents.The resulting square bound has the form ω ≤ 3τ for a parameter τ determined by an equation.
- Coppersmith–Winograd tensor: The Coppersmith–Winograd tensor decomposes into a sum of matrix multiplication tensors, but the sum is not direct, motivating zero-outs in tensor powers through the laser method.The method uses leveled partitions of tensor powers into constituent tensors.
4 Algorithm Outline
The algorithm follows a multistage interface-tensor framework, but introduces more asymmetry in hashing and compatibility handling. This requires six dimension permutations rather than the three used in prior work.
- Framework: The framework transforms tensor powers into independent interface tensors, then recursively applies global and constituent stages to obtain matrix multiplication tensors.Interface tensors carry constituent-tensor parameter lists and complete split distributions across levels.
- Asymmetry: The algorithm treats X, Y, and Z asymmetrically and applies its procedure across six permutations of the dimensions.Prior work treated two dimensions symmetrically and therefore used three regions.
- Asymmetric hashing: After asymmetric hashing, each level-ℓ X-block lies in a unique triple, while Y- and Z-block multiplicities are ordered rather than both being forced to uniqueness.This weaker intermediate property is sufficient for later compatibility zero-outs.
- Compatibility and holes: The method zeroes out incompatible and non-useful Y-blocks before applying analogous Z-compatibility steps, then fixes holes in the resulting broken interface tensors.The hole-fixing theorem applies when the hole fraction is sufficiently small and enough broken copies are available.
5 Global Stage
The global stage partitions the tensor power into six permutation-specific regions and applies asymmetric hashing followed by marginal, compatibility, and usefulness zero-outs. It produces many independent interface tensors subject to the method’s distributional constraints.
- Global-stage output: Theorem 5.3 states that the global stage yields 2^o(n) independent copies of a level-ℓ ε-interface tensor under the proposition’s constraints.These outputs are the intermediate objects passed to later stages.
- Region partitioning: The global stage partitions the tensor power into six regions because the method treats the three dimensions asymmetrically.Each region uses a different permutation of X, Y, and Z.
- Initial zero-out: The first zero-out retains blocks consistent with the marginal distributions αX, αY, and αZ of a chosen constituent-tensor distribution α.The remaining block counts are estimated through the entropies of these marginals.
- Asymmetric hashing: Asymmetric hashing guarantees unique triple containment for every level-ℓ X-block, unlike prior procedures that also guaranteed it for Y-blocks.The hashing procedure uses Salem–Spencer sets to eliminate incompatible block combinations.
5.3 Y -Compatibility Zero-Out
The Y-compatibility stage removes Y-blocks that fail split-distribution or uniqueness conditions, after which analogous Z processing and hole analysis prepare broken interface tensors for repair. The resulting copies have controlled holes and can be converted into unbroken copies.
- Y-compatibility: Y-blocks are first filtered for split-distribution consistency, then retained only when compatible with a unique level-ℓ triple.Compatibility checks include boundary cases and aggregate conditions over indices.
- Y-usefulness: Y-usefulness requires every relevant split distribution to match βY,i,j,k for the unique triple containing the block.Blocks failing any equality are zeroed out.
- Z processing: The analogous Z steps enforce compatibility and usefulness, but Z-blocks may initially belong to multiple triples and are zeroed out when this occurs.Afterward, each remaining Z-block is associated with one compatible triple.
- Hole fixing: The remaining subtensors are broken copies of the target interface tensor, with holes caused by compatibility zero-outs.The expected number of such copies is large enough for the hole-fixing theorem when the hole fraction is bounded.
5.6 Summary
The algorithm degenerates the tensor into independent copies of a level-ℓ interface tensor, using six regions associated with permutations of the X, Y, and Z dimensions.
- The algorithm degenerates the tensor into independent copies of a level-ℓ interface tensor.
- It applies the same procedure across six regions, permuting the roles of the X-, Y-, and Z-blocks in lexicographic order.The outputs from all six regions are combined by tensor product.
6 Constituent Stage
The constituent stage transforms level-ℓ interface tensors into a matrix multiplication tensor together with independent copies of lower-level interface tensors. Its formal construction uses split distributions, region weights, and recursively compatible parameters.
- The constituent stage combines a matrix multiplication tensor with independent copies of level-(ℓ−1) interface tensors.Terms with it = 0, jt = 0, or kt = 0 yield the matrix multiplication tensor; the remaining terms produce lower-level interface tensors.
- The construction defines distributions γX,t, γY,t, and γZ,t that split each level-ℓ index sequence into two level-(ℓ−1) sequences.
- The joint distribution αt has marginals consistent with the three γ distributions and is paired with level-(ℓ−1) complete split distributions βX,t, βY,t, and βZ,t.
- Proposition 6.3 degenerates an s′-term level-ℓ interface tensor into independent copies of a level-(ℓ−1) interface tensor under six-region parameter constraints.The regions are weighted by At,r, with nonnegative weights summing to one for each t.
- Theorem 6.4 extends the constituent-stage degeneration to 2^o(n) independent copies and preserves the lower-level interface-tensor structure.Its proof is stated to follow the corresponding theorem of VXXZ24.
6.1 Dividing into Regions
The algorithm divides each term into six regions, assigns region proportions and complete split distributions, and zeroes out blocks inconsistent with those distributions. The remaining tensor has the claimed structured form.
- Each term is divided into six regions with nonnegative proportions At,r summing to one.
- Each region receives complete split distributions for the X-, Y-, and Z-dimensions, subject to the conditions in Remark 6.1.
- The algorithm retains level-1 blocks only when their regional portions are ε-approximately consistent with the selected split distributions.The same consistency-based zeroing is applied to X-, Y-, and Z-blocks.
- After zeroing out, the remaining tensor is isomorphic to the structured tensor stated in Claim 6.5.
6.2 More Asymmetric Hashing
More asymmetric hashing first enforces marginal consistency, then uses independent hash functions and a progression-free bucket set to ensure that surviving block triples have controlled sharing. The resulting procedure satisfies the stated hashing properties.
- The algorithm first zeroes out blocks inconsistent with the marginal distributions γX,t, γY,t, and γZ,t.It also records the numbers NBX, NBY, and NBZ of surviving X-, Y-, and Z-blocks.
- The standard hashing step uses independent hash functions for X, Y, and Z blocks and retains only blocks whose hash values lie in a Salem-Spencer set B.The set B has no three-term arithmetic progressions and size M^1−o(1).
- The more asymmetric procedure zeroes out an X-block when multiple surviving triples share it and have the same hash value, then retains only triples consistent with {αt}t∈[s′].Every surviving X-block is contained in a unique consistent triple after this step.
- Lemma 6.7 states that same-bucket triples sharing an X-block cannot differ in both their Y- and Z-blocks, with analogous properties for shared Y- or Z-blocks.
- The surviving consistent triples satisfy the additional structural condition stated in the second property of Lemma 6.7.
6.3 Y -Compatibility Zero-Out
The Y-compatibility zero-out removes level-1 Y-blocks that fail prescribed split conditions, then retains only blocks associated uniquely and usefully with level-(ℓ−1) block triples.
- Compatibility filtering: The Y-compatibility step targets level-1 Y-blocks whose split distributions do not match the required βY values.The zero-out checks the prescribed split conditions across relevant indices and parts.
- Compatibility filtering: A remaining Y-block is compatible when it satisfies both the boundary-index and aggregate j′ split requirements.The two conditions require agreement for k′ = 0 cases and for every permitted j′ aggregate.
- Compatibility filtering: Claim 6.9 establishes that every remaining level-1 Y-block is compatible with the level-(ℓ−1) block triple containing it.This provides the correctness guarantee after the initial compatibility zero-out.
- Uniqueness and usefulness: Blocks compatible with multiple level-(ℓ−1) triples are zeroed out, leaving each surviving Y-block associated with a unique compatible triple.The subsequent zero-out resolves ambiguity before usefulness is checked.
- Uniqueness and usefulness: A Y-block is useful for its unique triple exactly when all its split values equal the corresponding βY values.The final usefulness zero-out removes blocks that fail any required split equality.
6.4 Z-Compatibility Zero-Out
The Z-compatibility zero-out applies analogous filtering to level-1 Z-blocks, using prescribed split conditions, uniqueness of triple membership, and a final usefulness test.
- Compatibility filtering: The Z-compatibility stage begins by removing level-1 Z-blocks whose split distributions disagree with required βZ values.The initial test checks the relevant aggregate k′ splits for every t.
- Compatibility filtering: A Z-block is compatible when it satisfies both the boundary-index conditions and the prescribed aggregate k′ split conditions.The boundary condition covers indices with i′ = 0 or j′ = 0, while the second covers each permitted k′.
- Compatibility filtering: Claim 6.12 states that every remaining level-1 Z-block is compatible with its containing level-(ℓ−1) block triple.The claim supplies the guarantee after the compatibility filtering step.
- Uniqueness and usefulness: Z-blocks compatible with multiple level-(ℓ−1) triples are zeroed out, so each survivor belongs to one compatible triple.This uniqueness step precedes the usefulness check.
- Uniqueness and usefulness: A Z-block is useful when every split value matches the corresponding βZ value for its unique triple.The final zero-out removes Z-blocks that are not useful for that triple.
6.5 Fixing Holes
The hole-fixing analysis classifies losses from input restrictions and compatibility zero-outs, then bounds compatibility ambiguity through probabilities for Y- and Z-blocks.
- Sources of holes: The analysis focuses on holes created when Y- and Z-blocks are compatible with multiple level-(ℓ−1) triples.These compatibility holes are analyzed alongside holes caused by enforcing approximate complete split distributions.
- Typical blocks: Typical Y- and Z-blocks are defined by requiring their complete split distributions to lie within ε in L∞ distance of the corresponding β distributions.The typicalness conditions are imposed for every t ∈ [s′].
- Compatibility probabilities: The quantities pY,comp and pZ,comp measure the probability that a uniformly random consistent level-(ℓ−1) triple is compatible with a fixed Y- or Z-block.Symmetry makes the Y compatibility probability independent of the chosen block and containing block.
- Compatibility probabilities: For Y-blocks, compatibility counting compares all consistent blocks with those additionally compatible with the triple, yielding pY,comp = Q/P.The proof bounds Q by relaxing the complete-split-distribution condition and uses the L∞ distance bound ε.
- Compatibility probabilities: The Y compatibility characterization partitions the relevant interface indices into k′ = 0 subsets and positive-k′ aggregate subsets, enabling the count of compatible blocks.The resulting possibilities are combined with the number of consistent triples to bound pY,comp.
6.6 Summary
The algorithm degenerates the input into independent copies of a lower-level interface tensor, and applies the same construction across all six regions with permuted dimension roles.
- Algorithmic conclusion: The algorithm degenerates the input into Ns′ independent copies of a level-(ℓ−1) interface tensor.The output tensor has the parameter list specified immediately after this conclusion.
- Algorithmic conclusion: The construction described for the first region extends to the other regions by permuting the roles of the X, Y, and Z dimensions.The final output is the tensor product of the outputs from all six regions.
7 Numerical Result
The paper combines a multilevel degeneration procedure, asymptotic sum inequalities, and nonlinear optimization to obtain bounds for rectangular matrix multiplication exponents. Its analysis gives ω ≤2.371339, while later optimization and scaling to the eighth CW-tensor power yield ω <2.371177.
- Degeneration: The procedure repeatedly degenerates interface tensors into Kronecker products with matrix multiplication tensors, eventually producing a direct sum of matrix multiplication tensors.The construction starts from independent level-ℓ* interface tensors and proceeds through lower levels until level 1.
- Exponent bounds: The resulting degeneration is combined with Schönhage’s asymptotic sum inequality to bound ω(1, κ, 1) after taking n →∞ and ε →0.The bound depends first on ε and then yields the limiting exponent bound.
- Optimization: The best parameters are found by formulating the construction as a nonlinear optimization problem subject to the constraints from Theorems 5.3, 6.2, and 6.4.SNOPT is used to seek feasible parameter settings; optimality is not guaranteed, but any sufficiently good feasible solution gives a useful exponent bound.
- Numerical results: The numerical program evaluates bounds for ω(1, κ, 1) across values of κ using the fourth power of the CW tensor with q = 5.The reported parameter sets and code are made available by the authors.
- Numerical results: ω <2.371177 results when modern optimization and machine-learning techniques find better parameters and scale the same analysis to the eighth CW-tensor power.The improved parameters are substituted into the existing analysis rather than requiring a different algorithmic framework.