New Bounds for Matrix Multiplication: from Alpha to Omega

Virginia Vassilevska Williams; Yinzhan Xu; Zixuan Xu; Renfei Zhou

Source-linked AI summary

New Bounds for Matrix Multiplication: from Alpha to Omega

Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, Renfei Zhou

arXiv:2307.07970v2Jul 16, 2023cs.DScs.CC

TL;DR

Matrix multiplication seeks faster algorithms approaching the ideal exponent ω = 2, but existing methods face limitations. This paper introduces an improved laser-method variant and obtains stronger bounds, including ω ≤ 2.371552 and α ≥ 0.321334.

Problem
Matrix multiplication seeks ω = 2 for near-linear algorithms, but existing laser-method variants cannot prove ω < 2.3078.
Method
The paper develops an improved laser-method variant to improve bounds for both square and rectangular matrix multiplication.
Results
ω ≤ 2.371552, α ≥ 0.321334, and µ ≤ 0.527661 are obtained, improving bounds for square and rectangular matrix multiplication.
Takeaways & Limitations
The paper establishes improved bounds across square and rectangular matrix multiplication, including the exponents ω, α, and µ.
Takeaways & Limitations
Known laser-method techniques and variants cannot show that ω < 2.3078.

Abstract

from arXiv · show

The main contribution of this paper is a new improved variant of the laser method for designing matrix multiplication algorithms. Building upon the recent techniques of [Duan, Wu, Zhou, FOCS 2023], the new method introduces several new ingredients that not only yield an improved bound on the matrix multiplication exponent $ω$, but also improve the known bounds on rectangular matrix multiplication by [Le Gall and Urrutia, SODA 2018]. In particular, the new bound on $ω$ is $ω\le 2.371552$ (improved from $ω\le 2.371866$). For the dual matrix multiplication exponent $α$ defined as the largest $α$ for which $ω(1,α,1)=2$, we obtain the improvement $α\ge 0.321334$ (improved from $α\ge 0.31389$). Similar improvements are obtained for various other exponents for multiplying rectangular matrices.

1 Introduction

The paper introduces an improved laser-method variant that advances bounds for square and rectangular matrix multiplication. It obtains improved values for ω, α, µ, and related rectangular exponents, with consequences for APSP and clique detection.

Motivation and limitations: ω = 2 is unattainable by known matrix-multiplication techniques, while laser-method variants cannot establish ω < 2.3078.These limitations motivate seeking improvements within the laser-method framework.
Rectangular matrix multiplication: Rectangular multiplication matters for APSP, minimum witnesses, bottleneck paths, and k-clique detection through exponents such as µ and ω(1, 2, 1).µ satisfies ω(1, µ, 1) = 1 + 2µ, while 4-clique detection runs in O(n^ω(1,2,1)+ε) time.
Main contribution: The paper presents a new improved laser-method variant yielding better bounds for both square and rectangular matrix multiplication.Its goal is to improve ω, α, µ, and rectangular matrix-multiplication bounds in general.
Main results: α > 0.321334, µ < 0.527661, and ω(1, 2, 1) < 3.250385, improving previous bounds 0.31389, 0.5286, and 3.25164 respectively.The previous rectangular bounds were given by Le Gall and Urrutia.
Algorithmic consequences: O(n^2.527661) time follows for Zwick’s APSP algorithm, and 4-cliques can be found in O(n^3.250385) time.The same APSP consequence applies to several related algorithms, including minimum witnesses and All-Pairs Bottleneck Paths.

2 Technical Overview

The paper strengthens the laser method by moving zeroing-out to finer level-1 blocks, enforcing split distributions in all three dimensions, and fixing the resulting holes simultaneously. These changes provide more degrees of freedom while requiring the hole fraction to be sufficiently small.

Complete split distribution: Level-1 independence gives more degrees of freedom than independence over level-(ℓ−1) blocks, enabling more copies of the tensor to be retained.The method moves zeroing out level-1 blocks according to γ earlier in the construction.
Enforcing split distributions in all three dimensions: The method enforces complete split distributions across X-, Y-, and Z-variables, extending Dual et al.’s one-dimensional enforcement.All three dimensions must be controlled for the new construction.
Hole analysis: 2^-Ω(n) is the resulting fraction of holes in the X-variable construction, up to ε additive error.Analogous arguments apply to the Y- and Z-variables.
Hole fixing: The generalized hole-fixing method repairs holes in all three dimensions simultaneously for a broad class of tensors.It requires the hole fraction to be below O(1/log N), where N is the number of tensor variables; the preceding step satisfies this condition.

3 Preliminaries · 3.1 Tensors and Tensor Operations

This section defines tensors as trilinear forms over three variable sets and introduces the basic operations used throughout the paper. It also formalizes independence, tensor powers, and isomorphism under variable relabeling.

3.1 Tensors and Tensor Operations: A tensor T is a trilinear form over variable sets X, Y, and Z with coefficients from a field F.The supports of T are the three variable sets X, Y, and Z.
3.1 Tensors and Tensor Operations: All tensors used in the paper have coefficients in {0, 1}, so they can be considered over any field F.This field-independence follows from the coefficient restriction stated in the preliminaries.
3.1 Tensors and Tensor Operations: The sum T + T′ is defined only when T and T′ have identical supports, and it adds corresponding coefficients.The common supports must satisfy (X, Y, Z) = (X′, Y′, Z′).
3.1 Tensors and Tensor Operations: The direct sum T ⊕ T′ first relabels variables to make the supports disjoint and then sums the tensors.For disjoint supports, the ordinary sum equals the direct sum, and the tensors are called independent.
3.1 Tensors and Tensor Operations: The notation T ⊕^n denotes the direct sum of n independent copies of T.This operation repeatedly combines copies over disjoint variable sets.
3.1 Tensors and Tensor Operations: The tensor product, or Kronecker product, T ⊗ T′ is formed over variable sets X × X′, Y × Y′, and Z × Z′.The notation T ⊗^n denotes the n-th tensor power of T.
3.1 Tensors and Tensor Operations: Tensors T and T′ are isomorphic, written T ≡ T′, when their support sizes match and coefficient arrays agree after permutations of the three variable sets.Thus, isomorphism identifies tensors equivalent up to relabeling variables.

3.2 Tensor Rank · 3.3 Degenerations, Restrictions, Zero-outs

The paper defines tensor rank and asymptotic rank, then introduces degenerations, restrictions, and zero-outs as operations that preserve appropriate rank upper bounds. Zero-outs are the restricted transformations used in the laser method.

3.2 Tensor Rank: Tensor rank R(T) is the minimum integer r ≥ 0 for which T admits a rank decomposition.The decomposition is the sum referenced in the definition.
3.2 Tensor Rank: Rank is subadditive under tensor addition: R(T + T′) ≤ R(T) + R(T′).
3.2 Tensor Rank: Asymptotic rank eR(T) is well-defined and satisfies eR(T) ≤ R(T ⊗m)^(1/m) for every fixed integer m > 0.This follows from the stated third rank property and Fekete’s lemma.
3.3 Degenerations, Restrictions, Zero-outs: A degeneration T ⊵ T′ is defined using F[λ]-linear maps between the tensor spaces over a common field F.The supplied definition introduces F[λ] as the polynomial ring in the formal variable λ.
3.3 Degenerations, Restrictions, Zero-outs: Degeneration cannot increase asymptotic rank: if T ⊵ T′, then eR(T′) ≤ eR(T).
3.3 Degenerations, Restrictions, Zero-outs: A restriction is a degeneration whose maps φX, φY, and φZ are F-linear.
3.3 Degenerations, Restrictions, Zero-outs: Restrictions preserve both rank and asymptotic-rank upper bounds: R(T′) ≤ R(T) and eR(T′) ≤ eR(T).
3.3 Degenerations, Restrictions, Zero-outs: The laser method uses zero-outs, which restrict variables by setting selected variables outside X′, Y′, and Z′ to zero.The resulting tensor T′ = T|X′,Y′,Z′ is called a subtensor over X′, Y′, and Z′.

3.4 Matrix Multiplication Tensors

This section defines matrix multiplication tensors and explains how their rank yields algorithms for square and rectangular matrix multiplication. It introduces the exponents ω and ω(a, b, c), focusing on bounds for ω(1, κ, 1).

Matrix Multiplication Tensors: The tensor ⟨a, b, c⟩ computes the product of an a × b matrix and a b × c matrix, producing an a × c matrix.It is defined over variable sets corresponding to the input and output matrix entries.
Matrix Multiplication Tensors: Tensor products satisfy ⟨a, b, c⟩⊗⟨d, e, f⟩ ≡ ⟨ad, be, cf⟩.
Square Matrix Multiplication: The matrix multiplication exponent ω is defined through the rank of square matrix multiplication tensors and determines circuits of size O(n^(ω+ε)) for every ε > 0.If R(⟨q, q, q⟩) ≤ r, recursive rank decompositions give circuits of size O(n^log_q(r)).
Rectangular Matrix Multiplication: For rectangular multiplication, ω(a, b, c) gives circuits of size O(n^(ω(a,b,c)+ε)) for multiplying matrices of sizes n^a × n^b and n^b × n^c.The paper focuses on values of the form ω(1, κ, 1) for κ > 0, with ω(1, 1, κ) = ω(1, κ, 1) = ω(κ, 1, 1).

3.5 Schönhage’s Asymptotic Sum Inequality · 3.6 The Coppersmith-Winograd Tensor

Section 3.5 explains how Schönhage’s asymptotic sum inequality converts asymptotic-rank bounds for direct sums of matrix multiplication tensors into bounds on ω and rectangular exponents. Section 3.6 defines the Coppersmith-Winograd tensor, identifies its tensor decomposition, and states its known asymptotic-rank bound.

3.5 Schönhage’s Asymptotic Sum Inequality: Asymptotic-rank bounds for matrix multiplication tensors yield upper bounds on ω.Schönhage’s result extends this approach to direct sums of matrix multiplication tensors.
3.5 Schönhage’s Asymptotic Sum Inequality: ω ≤ 3τ, where τ ∈ [2/3, 1] solves the stated equation.This is the quantitative conclusion of Schönhage’s asymptotic sum inequality.
3.5 Schönhage’s Asymptotic Sum Inequality: The asymptotic sum inequality also gives bounds for rectangular matrix multiplication.A separate theorem states the rectangular version for ω(a, b, c).
3.6 The Coppersmith-Winograd Tensor: CW_q is defined for nonnegative integer q ≥ 0 over variables X, Y, and Z, each indexed from 0 through q + 1.The variable sets are X = {x_0, …, x_{q+1}}, Y = {y_0, …, y_{q+1}}, and Z = {z_0, …, z_{q+1}}.
3.6 The Coppersmith-Winograd Tensor: The defining expression includes the terms x_0y_iz_i + x_iy_0z_i + x_iy_iz_0.These terms form part of the Coppersmith-Winograd tensor’s construction.
3.6 The Coppersmith-Winograd Tensor: CW_q is the sum of six matrix multiplication tensors, with the other three being copies of ⟨1, 1, 1⟩.This decomposition is used in the tensor’s analysis.
3.6 The Coppersmith-Winograd Tensor: eR(CW_q) ≤ q + 2 is the known bound due to Coppersmith and Winograd.The bound applies to the Coppersmith-Winograd tensor defined in this section.

3.7 Base Leveled Partition of CWq

The section defines a leveled partition of tensor powers of CW_q by first partitioning CW_q into three variable parts and then coarsening sequence-indexed blocks by their coordinate sums. The resulting level-ℓ constituent structure is characterized by the condition i + j + k = 2ℓ.

Tensor-power construction: T^(ℓ) is defined as the (2ℓ−1)-th tensor power of CW_q for ℓ ≥ 1.
Level-1 Partition: CW_q is partitioned at level 1 into three parts for each variable set X^(1), Y^(1), and Z^(1).The corresponding subtensors are called level-1 constituent tensors.
Level-1 Partition: At level 1, a constituent tensor is nonzero if and only if i + j + k = 2.
Higher-level partitions: The level-1 partition induces partitions of higher tensor powers indexed by {0, 1, 2}-sequences of length 2ℓ−1.The same construction applies to the X-, Y-, and Z-variable sets.
Higher-level partitions: The level-ℓ partition merges sequence-indexed parts with the same sum, yielding level-ℓ variable blocks whose nonzero constituents satisfy i + j + k = 2ℓ.This coarsening can also be viewed recursively from the level-(ℓ−1) partition.

3.8 Leveled Partition for Large Tensor Powers of CWq

The section defines leveled partitions for large tensor powers of CWq at levels 1 and ℓ. It characterizes nonzero subtensors through index-sum conditions and relates fine level-1 blocks to coarser level-ℓ blocks.

Level-1 partition: At level 1, variable blocks are indexed by length-N sequences in {0, 1, 2}^N, and compatible triples satisfy Î_t + Ĵ_t + K̂_t = 2 for every t.Such blocks define level-1 triples that partition (T (1))^⊗N.
Level-ℓ partition: At level ℓ, (CWq)^⊗N is viewed as (T (ℓ))^⊗n with n = N/2^(ℓ−1), and blocks are indexed by length-n sequences in {0, 1, . . . , 2^ℓ}.The level-ℓ partition is induced from the partition of T (ℓ).
Level-ℓ partition: Level-ℓ blocks form nonzero triples exactly when I_t + J_t + K_t = 2^ℓ for every t, denoted X_IY_JZ_K.These subtensors are called level-ℓ triples.
Relation between levels: Because the level-ℓ partition coarsens the level-1 partition, a level-1 block X_Î is contained in X_I when summing consecutive length-2^(ℓ−1) subsequences of Î yields I.This containment is written Î ∈ I and X_Î ∈ X_I.

3.9 Distributions and Entropy

This section defines finite-support distributions and their base-2 entropy, then introduces a combinatorial lemma and joint distributions. For integer-sequence supports, the joint distribution is formed through concatenation.

A distribution has finite support, assigns nonnegative probabilities to its support, and the probabilities sum to 1.
Entropy H(α) is defined using logarithms with base 2.
Lemma 3.3 states a well-known combinatorial fact for a distribution over [s] and a positive integer N.
For distributions α and β on S and S′, α × β is defined over pairs in S × S′ according to the product construction.
When supports contain integer sequences, α × β instead ranges over concatenations of sequences, with s ◦ s′ denoting concatenation.

3.10 Complete Split Distributions

This section defines complete split distributions to describe how level-1 variable blocks occur within level-ℓ constituent tensors. It also formalizes consistency, approximation, distance, subset restriction, and tensor notation based on these distributions.

Definition and motivation: A complete split distribution for T_i,j,k with i + j + k = 2^ℓ assigns probabilities to all length-2^ℓ−1 sequences over {0, 1, 2}.These distributions characterize the level-1 variable blocks contained in level-ℓ variable blocks.
Consistency: A level-1 index sequence is consistent with β when each length-2^ℓ−1 index type occurs with the proportion specified by β.The definition checks every sequence in {0, 1, 2}2^ℓ−1.
Induced distributions: Every level-1 index sequence induces a complete split distribution by measuring the proportions of its length-2^ℓ−1 consecutive chunks.The same construction can be applied after restricting the sequence to a subset S ⊆ [n].
Approximate consistency and tensor notation: Approximate consistency allows L∞ error at most ε, where the distance between β1 and β2 is the maximum coordinate-wise difference over all index sequences.For a level-ℓ constituent tensor, the notation aggregates level-1 X-, Y-, and Z-variable blocks approximately consistent with βX, βY, and βZ.

3.11 Salem-Spencer Sets

The hashing step uses a large dense subset of ZM that avoids nontrivial 3-term arithmetic progressions. Such a set B exists for every positive integer M.

Hashing step: The laser method’s hashing step uses a large dense subset of ZM with no nontrivial 3-term arithmetic progressions.This ingredient is recalled from prior work.
Salem-Spencer property: For every positive integer M > 0, there exists a subset B ⊆ ZM containing no nontrivial 3-term arithmetic progressions.For a, b, c ∈ B, a + b ≡ 2c (mod M) holds if and only if a = b = c.

4 Algorithm Outline

The algorithm degenerates a large tensor power into independent matrix multiplication tensors by passing through a hierarchy of interface tensors. Its stages enforce compatibility and independence through asymmetric zero-outs, then repair the resulting holes.

Algorithm framework: The algorithm transforms a large tensor power into independent matrix multiplication tensors of dimensions ⟨m, mκ, m⟩, yielding an upper bound on ω(1, κ, 1).The first level accepts a large tensor power, while subsequent levels reduce interface tensors until the target matrix multiplication tensors are obtained.
Interface tensors: Interface tensors connect algorithmic levels, and their tensor products remain interface tensors with concatenated parameter lists.A level-ℓ interface tensor is built from constituent tensors and complete split distributions for the X-, Y-, and Z-variables.
Algorithm framework: The framework first applies the global stage, then recursively applies constituent tensor stages from level ℓ* down to level 2, and finally degenerates level-1 interface tensors.Each constituent stage outputs a matrix multiplication tensor together with independent copies of a lower-level interface tensor.
Global-stage procedure: The global-stage procedure partitions the input into three regions, balances variable blocks across dimensions, and performs the subsequent procedure separately up to rotation.The regions satisfy A1 + A2 + A3 = 1, with each allowing level-ℓ block sharing in a different dimension.
Global-stage procedure: The procedure uses distribution-based zero-outs, asymmetric hashing, compatibility tests, usefulness filtering, and hole fixing to obtain level-1-independent interface tensors.Asymmetric cleanup makes X- and Y-blocks belong to unique level-ℓ triples, enabling later compatibility and usefulness checks.
Hole fixing: If each of r broken copies has at most a 1/(8N) fraction of level-1 blocks as holes and r ≥ 2C1N/log N, their direct sum degenerates into an unbroken copy.This hole-fixing guarantee is supplied by Corollary 4.2 for sufficiently large C1 > 0.

5 Global Stage

The global stage converts multiple copies of a tensor power into independent copies of a level-ℓε-interface tensor using asymmetric hashing across three regions. Its analysis combines marginal-consistent block counting, Salem–Spencer filtering, and control of holes in the resulting copies.

Global-stage guarantee: Theorem 5.3 yields 2^o(n) independent copies of a level-ℓε-interface tensor from the prescribed tensor-power input, under the Proposition 5.1 constraints.The output uses a perturbed interface parameterization and incurs a small 2^o_{1/ε}(n) loss in the number of copies retained.
Regional hashing: The three tensor-product regions are analyzed independently using asymmetric hashing that shares Z-blocks, Y-blocks, and X-blocks in the first, second, and third regions, respectively.Each region is degenerated into independent copies of a level-ℓ interface tensor, whose tensor product forms the output.
Block counting: N_{αX,αY,αZ} = 2^(H(α)+P_α)·A_1n±o(n) counts the block triples surviving the marginal-consistency filtering.Here P_α is the entropy penalty from maximizing over distributions with the same three marginals.
Hash filtering and useful tensors: A Salem–Spencer subset of size M^(1−o(1)) filters hashed blocks while excluding nontrivial 3-term arithmetic progressions modulo M.The subsequent union-bound analysis gives copies of T* with hole fraction at most 1/8N, and their overall count is N_α·M^(−1−o(1)).

6 Constituent Stage

The constituent stage handles boundary terms that are already matrix multiplication tensors, then applies a three-region splitting and hashing procedure to degenerate level-ℓ interface tensors into independent lower-level copies. The analysis guarantees many surviving block triples while bounding holes from split-distribution constraints.

Boundary terms: Boundary terms with it = 0, jt = 0, or kt = 0 are handled as matrix multiplication tensors before the laser method is applied.When kt = 0, the tensor is isomorphic to an inner product tensor ⟨1, M, 1⟩ for some M ≥ 0.
Degeneration guarantee: A level-ℓ interface tensor can be degenerated into independent copies of a level-(ℓ−1) interface tensor under the parameter constraints of Proposition 6.2.Theorem 6.3 extends this to 2^o(n) independent input copies, yielding 2^(E1+E2+E3−o(n)−o_1/ε(n)) independent output copies.
Three-region decomposition: The constituent stage divides each remaining term into three regions with weights At,1, At,2, At,3 satisfying At,1 + At,2 + At,3 = 1.Each region uses its own complete split distributions, and incompatible level-1 blocks are zeroed out.
Hashing analysis: N^α · M^−1−o(1) copies of T* are expected after hashing, with the fraction of holes caused by uniqueness and input split-distribution constraints equal to 1/n^2.This combines the surviving compatible block triples with the two stated hole sources.

7 Fixing Holes

This section presents a hole-fixing method that degenerates sufficiently many broken copies of a partitioned tensor into an unbroken copy, under symmetry and small-hole conditions. Applied to interface tensors, the method yields an unbroken tensor from a direct sum of broken copies with at most a 1/8N fraction of holes in each dimension.

General hole-fixing theorem: The method degenerates sub-polynomially many broken copies with small hole fractions in all three dimensions into an unbroken copy of the tensor.The result generalizes a result by Duan and applies in a more general setting.
General hole-fixing theorem: The proof uses partition-preserving, tensor-structure-preserving permutations whose uniformly random elements send any given part to a uniformly random part.These conditions constitute Property 7.1 and support the permutation lemma used in the recursive construction.
General hole-fixing theorem: Theorem 7.2 requires each broken copy to have at most 1/4 log MX, 1/4 log MY, and 1/4 log MZ fractions of holes in the X-, Y-, and Z-parts, respectively.Here M = max{MX, MY , MZ}.
Proof strategy: The recursive construction first covers most terms with one broken copy, decomposes the missing terms into 7 smaller subtensors, and recursively handles those subproblems.The resulting recursion bounds the number of broken copies needed to produce a complete copy of T.
Application to interface tensors: For a level-ℓ interface tensor, a direct sum of r broken copies degenerates into an unbroken copy when each dimension has at most a 1/8N fraction of holes and r ≥ 2C1N/log N.The corollary applies the general theorem to level-1 blocks of the interface tensor.

8 Numerical Result

The numerical optimization uses MATLAB with SNOPT to evaluate bounds for ω(1, κ, 1) across different κ. Analyzing the fourth power of the CW tensor with q = 5 yields the key bounds ω ≤2.371552, α ≥0.321334, and µ ≤0.527661.

Optimization strategy: The κ = 1 optimization starts from a feasible solution matching Le Gall’s upper bound on ω.Sequential quadratic programming is then applied to obtain an upper bound for ω = ω(1, 1, 1).
Optimization strategy: For κ ≠ 1, each optimization is initialized from a solution for a nearby κ and proceeds sequentially across κ values.For example, the solution for ω(1, 1, 1) initializes the optimization for ω(1, 0.95, 1), followed by ω(1, 0.90, 1).
Numerical results: ω ≤2.371552, α ≥0.321334, and µ ≤0.527661 are the key numerical bounds obtained from the optimization.All bounds analyze the fourth power of the CW tensor with q = 5.
Numerical results: Bounds for ω(1, κ, 1) are computed for different κ using a MATLAB program assisted by SNOPT.The resulting values are reported in Table 1.