Paper deep dive
Improving the matrix multiplication exponent with modern optimization and AlphaEvolve
Emilien Dupont, Marvin Eisenberger, Borislav Kozlovskii, Abbas Mehrabian, Francisco J. R. Ruiz, Abigail See, Renfei Zhou, Josh Alman, Virginia Vassilevska Williams, Matej Balog
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 89%
Last extracted: 8/23/2026, 3:02:42 AM
Summary
This paper improves the upper bound on the matrix multiplication exponent omega from 2.371339 to 2.371177. The authors achieve this by reformulating the optimization problem inherent in combination loss analysis, leveraging machine learning-based gradient descent to solve it at a higher recursion level (l*=4), and refining the solution using AlphaEvolve.
Entities (12)
Relation Signals (8)
Strassen → demonstrated → sub-cubic time matrix multiplication
confidence 95% · The pioneering work of Strassen (22) showed that two n×n matrices can be multiplied in sub-cubic time
combination loss analysis → refines → laser method
confidence 95% · The current best bounds on the matrix multiplication exponent ω are obtained through a refinement of the laser method called combination loss analysis
JAX → implements → Gradient-based optimization
confidence 90% · our gradient-based optimization, implemented in Jax
AlphaEvolve → improves → optimization algorithm
confidence 90% · we refine the resulting optimization algorithm with AlphaEvolve
Gradient Descent → solves → optimization problem
confidence 90% · address the non-convex optimization problem of combination loss analysis using a gradient descent approach
Duan et al. → contributedto → combination loss analysis
confidence 80% · combination loss analysis (Duan et al., 2022
Williams et al. → contributedto → combination loss analysis
confidence 80% · combination loss analysis (Duan et al., 2022; Williams et al., 2024
Alman et al. → contributedto → combination loss analysis
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The current best bounds on the matrix multiplication exponent $\omega$ are obtained through a refinement of the laser method called combination loss analysis (Duan et al., 2022; Williams et al., 2024; Alman et al., 2025). In this note, we address the optimization problem at the core of this approach and propose several improvements. First, we reformulate the optimization problem allowing us to solve it in a larger setting than was previously possible. Second, we leverage recent advances in machine learning to design a new optimization algorithm for this problem. Finally, we refine the resulting optimization algorithm with AlphaEvolve. Our combined approach yields an upper bound of $\omega$ < 2.371177, improving the previous best bound of 2.371339.
Tags
Links
- Source: https://arxiv.org/abs/2608.16884v1
- Canonical: https://arxiv.org/abs/2608.16884v1
Trouble viewing inline? Open PDF directly →
Full Text
40,941 characters extracted from source content.
Expand or collapse full text
Improving the matrix multiplication exponent with modern optimization and AlphaEvolve Emilien Dupont Affiliation: Equal contribution in alphabetical order Affiliation: Marvin Eisenberger Affiliation: Equal contribution in alphabetical order Affiliation: Borislav Kozlovskii Affiliation: Equal contribution in alphabetical order Affiliation: Abbas Mehrabian Affiliation: Equal contribution in alphabetical order Affiliation: Francisco J. R. Ruiz Affiliation: Equal contribution in alphabetical order Affiliation: Abigail See Affiliation: Equal contribution in alphabetical order Affiliation: Renfei Zhou Affiliation: Equal contribution in alphabetical order Affiliation: Carnegie Mellon University Josh Alman Affiliation: Columbia University Virginia Vassilevska Williams Affiliation: MIT Matej Balog Affiliation: Abstract The current best bounds on the matrix multiplication exponent ω are obtained through a refinement of the laser method called combination loss analysis (10; 23; 1). In this note, we address the optimization problem at the core of this approach and propose several improvements. First, we reformulate the optimization problem allowing us to solve it in a larger setting than was previously possible. Second, we leverage recent advances in machine learning to design a new optimization algorithm for this problem. Finally, we refine the resulting optimization algorithm with AlphaEvolve. Our combined approach yields an upper bound of ω<2.371177ω<2.371177, improving the previous best bound of 2.3713392.371339. 1 Introduction From accelerating machine learning computations to enabling realistic computer graphics, matrix multiplication is a fundamental operation underpinning critical applications in computer science. Despite its prominence, the computational complexity of matrix multiplication—the number of arithmetic operations needed to multiply large matrices—is unknown, and determining it is a major open question in theoretical computer science (4). The pioneering work of Strassen (22) showed that two n×n× n matrices can be multiplied in sub-cubic time—specifically, (nω+o(1))O(n^ω+o(1)) operations for ω<2.81ω<2.81—spurring a line of work attempting to further reduce the complexity exponent ω (16; 3; 18; 17; 6; 21; 7; 20; 24; 14; 2; 10; 23; 1). All improvements in the past 40 years rely on the laser method, a mathematical technique to indirectly design matrix multiplication algorithms. The current best bound is achieved by a refinement of the laser method called combination loss analysis (10), a technique that requires solving a non-convex optimization problem as part of the computer-assisted proof, and yields ω<2.371339ω<2.371339 (1). Here, we improve the bound to ω<2.371ω<2.371177 using a two-step approach. First, we leverage recent advances in machine learning and adjacent areas to address the non-convex optimization problem of combination loss analysis using a gradient descent approach; this alone improves the previous state-of-the-art (SOTA) bound by ≈0.97×10−4≈ 0.97× 10^-4. Second, we use AlphaEvolve (15) to improve our optimization algorithm; this raises the improvement over the SOTA to ≈1.62×10−4≈ 1.62× 10^-4. 10 2.3712.371866 23 2.3712.371552 1 2.3712.371339 This note 2.3712.371177 Table 1: Recent improvements to ω. The optimization problem at the core of combination loss analysis is formulated in 1, where it was also shown that any feasible solution provides an upper bound on ω. To achieve our new bound, we target a slightly different optimization problem. Specifically, combination loss analysis has a parameter—the maximum recursion level, denoted by ℓ∗ ^*—that introduces a trade-off between the complexity of the optimization (which grows doubly exponentially in ℓ∗ ^*) and the best possible bound on ω it can achieve. The previous SOTA bound was found with ℓ∗=3 ^*=3 (1). In contrast, our gradient-based optimization, implemented in Jax (5), allows for hardware parallelization and can handle ℓ∗=4 ^*=4 (the number of optimizable parameters increases from approximately 25k to 7 million when moving from ℓ∗=3 ^*=3 to ℓ∗=4 ^*=4). In this note, we describe the full optimization problem from 1 in detail, explain how we numerically solved it and applied AlphaEvolve to it, and finally show how we rigorously certified the resulting omega bound. 2 Optimization problem The high-level structure of the optimization problem can be captured by a rooted tree, where every node is associated with a collection of optimizable parameters. Integers q≥1q≥ 1 and ℓ∗≥2 ^*≥ 2 are fixed as hyperparameters. Intuitively, the tree structure describes a recursive way to decompose the tensor CWq⊗2ℓ∗CW_q 2 ^* into smaller tensors, where CWqCW_q is the Coppersmith-Winograd tensor (7); see 1 for more details about this correspondence. We start by defining the tree structure, and will then define the optimizable parameters associated to its nodes. 2.1 Tree structure Denote [k]≔1,…,k[k] \1,…,k\. Each non-root node T in the tree is associated with a level, a shape sTs_T, and a region rTr_T, defined as follows. • The level of a non-root node is a positive integer ℓ∈2,3,…,ℓ∗ ∈\2,3,…, ^*\, describing its depth in the tree; higher is closer to the root. Direct children of the root node have level ℓ∗ ^*, which is a hyperparameter fixed in advance. • A level-ℓ shape is a triple s=(sX,sY,sZ)s=(s_X,s_Y,s_Z), where sX,sY,sZs_X,s_Y,s_Z are non-negative integers summing to 2ℓ2 . We use ℓ≔(i,j,k)∈ℤ≥03|i+j+k=2ℓS_ \(i,j,k) _≥ 0^3\; |\;i+j+k=2 \ to denote the set of all level-ℓ shapes. We use dimensions X,Y,ZX,Y,Z to refer to the indices of the three coordinates in a shape. • There are six regions, indexed by an integer r∈[6]r∈[6]. Each region is associated with πr _r, the r-th permutation over symbols X,Y,Z\X,Y,Z\ in lexicographic order. We call a non-root node T a positive-shape node if all coordinates in its shape sTs_T are positive. Otherwise it is a zero-shape node. The root node is denoted by G and does not have the aforementioned attributes—level, shape, and region. Next, we define the tree structure by describing the child nodes for different types of nodes. • Root: For every region r∈[6]r∈[6] and level-ℓ∗ ^* shape s∈ℓ∗s _ ^*, the root G has a child at level ℓ∗ ^* with shape s and region r, denoted as G[s,r]G[s,r]. • Positive-shape node: For every level-ℓ shape s=(sX,sY,sZ)s=(s_X,s_Y,s_Z), we define Split(s)≔u∈ℓ−1| 0≤uX≤sX, 0≤uY≤sY, 0≤uZ≤sZ. (s) \u _ -1\; |\;0≤ u_X≤ s_X,\;0≤ u_Y≤ s_Y,\;0≤ u_Z≤ s_Z \. Fixing a positive-shape node T at level ℓ≥3 ≥ 3, for every region r∈[6]r∈[6] and level-(ℓ−1)( -1) shape u∈Split(sT)u∈ (s_T), the node T has a child at level ℓ−1 -1 with shape u and region r, denoted as T[u,r]T[u,r]. Note that the region index r of the child can be different from that of T itself. • Zero-shape nodes and level-2 positive-shape nodes do not have child nodes. We call them the leaves. 2.2 Optimizable parameters Next we list the free variables of the optimization problem, grouped by the associated node on the tree. For any finite set D, we use Δ(D) (D) to denote the simplex of probability distributions on D: Δ(D)≔p:D→[0,1]|∑x∈Dp(x)=1 (D) \p:D→[0,1]\; |\; _x∈ Dp(x)=1 \ Root node. The optimizable parameters associated with the root node G are: • A distribution AG=(AG(1),…,AG(6))∈Δ([6])A_G=(A_G^(1),…,A_G^(6))∈ ([6]) over the six regions; • For r∈[6]r∈[6], a distribution αG(r)∈Δ(ℓ∗) _G^(r)∈ (S_ ^*) over all level-ℓ∗ ^* shapes. Positive-shape node. Each level-ℓ positive-shape node T with ℓ≥3 ≥ 3 is associated with the following parameters: • A distribution AT=(AT(1),…,AT(6))∈Δ([6])A_T=(A_T^(1),…,A_T^(6))∈ ([6]) over the six regions; • For r∈[6]r∈[6], a distribution αT(r)∈Δ(Split(sT)) _T^(r)∈ ( (s_T)), where sTs_T is the shape of T. Zero-shape node. For a level ℓ and an integer a with 0≤a≤2ℓ0≤ a≤ 2 , we define ℓ,a≔L∈0,1,22ℓ−1|∑p=12ℓ−1Lp=a.C_ ,a \L∈\0,1,2\^2 -1\; |\; _p=1^2 -1L_p=a \. This is the set of length-2ℓ−12 -1 vectors over 0,1,2\0,1,2\ whose entries sum to a. A distribution over ℓ,aC_ ,a is called a level-ℓ complete split distribution of a. For every level-ℓ zero-shape node T, let W∈X,Y,ZW∈\X,Y,Z\ be the first dimension where sT,Ws_T,W is nonzero. T is then associated with a complete split distribution βT,W∈Δ(ℓ,sT,W) _T,W∈ (C_ ,s_T,W). Level-2 node. The only remaining nodes are level-2 nodes with shapes (1,1,2)(1,1,2), (1,2,1)(1,2,1), or (2,1,1)(2,1,1), which are the only valid strictly positive shapes for level 2, since their coordinates must sum to 22=42^2=4. Each such node is associated with a scalar μT∈[0,1/2] _T∈[0,1/2]. 2.3 Derived quantities With the free variables fixed, all remaining quantities can be computed deterministically. In the following, all logarithms and entropies are in base two. For a distribution ρ with finite support, we denote its entropy by H(ρ)≔−∑x∈supp(ρ)ρ(x)logρ(x).H(ρ) - _x∈ (ρ)ρ(x) ρ(x). We also use H(p1,p2,…,pk)H(p_1,p_2,…,p_k) to denote the entropy of a distribution over [k][k] with probability masses p1,…,pkp_1,…,p_k. If D is a set of shapes and ρ∈Δ(D)ρ∈ (D), let ρW _W denote the marginal of ρ in coordinate W, namely, ρW(w)≔∑a∈D:aW=wρ(a) _W(w) _a∈ D:a_W=wρ(a). Define HDmax(ρ)≔supρ′∈Δ(D)ρW′=ρW for W∈X,Y,ZH(ρ′),H_D (ρ) _ subarraycρ ∈ (D)\\ ρ _W= _W for W∈\X,Y,Z\ subarrayH(ρ ), as well as the penalty notion PD(ρ)≔HDmax(ρ)−H(ρ).P_D (ρ ) H_D (ρ )-H (ρ ). Masses. Each non-zero-shape node T of the tree has a real number mT∈[0,1]m_T∈[0,1] associated with it, called its mass. The masses of non-root nodes are computed top-down as follows: • The root G has mG=1m_G=1. • For the children of the root, we set mG[s,r]≔AG(r)⋅αG(r)(s),∀s∈ℓ∗,r∈[6].m_G[s,r] A_G^(r)· _G^(r)(s), ∀\,s _ ^*,\ r∈[6]. • For any positive-shape node T of level ℓ≥3 ≥ 3, we set mT[u,r]=mT⋅AT(r)⋅(αT(r)(u)+αT(r)(sT−u)),∀u∈Split(sT),r∈[6].m_T[u,r]=m_T· A_T^(r)· ( _T^(r)(u)+ _T^(r)(s_T\!-\!u) ), ∀\,u∈ (s_T),\ r∈[6]. Complete split distributions. Every level-ℓ non-root node T carries, for each dimension W∈X,Y,ZW∈\X,Y,Z\, a complete split distribution βT,W∈Δ(ℓ,sT,W) _T,W∈ (C_ ,s_T,W ). These distributions are calculated as follows. Let 0→≔(0,0,…,0) 0 (0,0,…,0) and 2→≔(2,2,…,2) 2 (2,2,…,2) denote vectors of length 2ℓ−12 -1. For a zero-shape node T, let W0∈X,Y,ZW_0∈\X,Y,Z\ be the first zero coordinate of sTs_T, W1W_1 be the first nonzero coordinate of sTs_T, and W2W_2 be the other coordinate. βT,W0 _T,W_0 is the point mass distribution at the length-2ℓ−12 -1 vector 0→ 0; βT,W1 _T,W_1 was defined as optimizable parameters in section 2.2; βT,W2≔βT,W1∨ _T,W_2 β _T,W_1, where for any complete split distribution β, we define β∨(2→−L)≔β(L),∀L∈supp(β).β ( 2\!-\!L) β(L), ∀\,L∈ (β). (Notice, in particular, that if the components of L sum to sT,W1s_T,W_1, then the components of 2→−L 2-L sum to 2⋅2ℓ−1−sT,W1=sT,W22· 2 -1-s_T,W_1=s_T,W_2, keeping the mapping validly within the correct domain.) For a positive-shape node T at level ℓ≥3 ≥ 3, we define βT,W(r)≔∑u∈Split(sT)αT(r)(u)⋅(βT[u,r],W×βT[sT−u,r],W),∀r∈[6], _T,W^(r) _u∈ (s_T) _T^(r)(u)· ( _T[u,r],W× _T[s_T-u,r],W ), ∀\,r∈[6], where × denotes the Cartesian product of complete split distributions; that is, if LleftL_left and LrightL_right are two sequences of length 2ℓ−22 -2 in supports of βleft≔βT[u,r],W _left _T[u,r],W and βright≔βT[sT−u,r],W _right _T[s_T-u,r],W, respectively, then we form a new sequence L, of length 2ℓ−12 -1, by concatenating LleftL_left and LrightL_right, and letting (βleft×βright)(L)≔βleft(Lleft)⋅βright(Lright).( _left× _right)(L) _left(L_left)· _right(L_right). Then, the complete split distributions of T are calculated by βT,W≔∑r=16AT(r)βT,W(r). _T,W _r=1^6A_T^(r) _T,W^(r). The only remaining case is a level-2 positive-shape node, which must have shape (1,1,2)(1,1,2), (2,1,1)(2,1,1), or (1,2,1)(1,2,1). For a node T with shape (1,1,2)(1,1,2), recall that μT∈[0,1/2] _T∈[0,1/2] is the only optimizable parameter associated with T. We let δa,b _a,b denote the point mass at (a,b)∈0,1,22(a,b)∈\0,1,2\^2 and set βT,X=βT,Y=12δ0,1+12δ1,0,βT,Z=μTδ0,2+μTδ2,0+(1−2μT)δ1,1. _T,X= _T,Y= 12 _0,1+ 12 _1,0, _T,Z= _T _0,2+ _T _2,0+(1-2 _T) _1,1. Similarly, for a node T with shape (2,1,1)(2,1,1), we set βT,Y=βT,Z=12δ0,1+12δ1,0,βT,X=μTδ0,2+μTδ2,0+(1−2μT)δ1,1, _T,Y= _T,Z= 12 _0,1+ 12 _1,0, _T,X= _T _0,2+ _T _2,0+(1-2 _T) _1,1, and for a node T with shape (1,2,1)(1,2,1), we set βT,X=βT,Z=12δ0,1+12δ1,0,βT,Y=μTδ0,2+μTδ2,0+(1−2μT)δ1,1. _T,X= _T,Z= 12 _0,1+ 12 _1,0, _T,Y= _T _0,2+ _T _2,0+(1-2 _T) _1,1. Retained exponent at root node. Next, we define how to calculate the retained exponent EGE_G associated with the root G, which will be obtained from intermediate quantities EG(r)E_G^(r) for regions r∈[6]r∈[6]. We start by introducing how to compute EG(1)E^(1)_G for region r=1r=1, where π1 _1 is the identity permutation. In all summations below, s ranges over ℓ∗S_ ^*. We will use the symbols *,+,*,+, and < in our notation below as mnemonic placeholders corresponding to the dimensions, indicating the meanings “unconstrained / any value”, “strictly positive”, and “strictly less than the parent’s coordinate”, respectively. We define ηG,Y(1)≔∑s:sZ=0αG(1)(s)⋅H(βG[s,1],Y)+∑j=02ℓ∗αG(1)(*,j,+)⋅H(β¯G,Y,*,j,+(1)), _G,Y^(1) _s\,:\,s_Z=0 _G^(1)(s)· H ( _G[s,1],Y )+ _j=0^2 ^* _G^(1)(*,j,+)· H ( β_G,Y,*,j,+^(1) ), where αG(1)(*,j,+)≔∑s:sY=j,sZ>0αG(1)(s),β¯G,Y,*,j,+(1)≔1αG(1)(*,j,+)∑s:sY=j,sZ>0αG(1)(s)⋅βG[s,1],Y. _G^(1)(*,j,+) _s\,:\,s_Y=j,\;s_Z>0 _G^(1)(s), β_G,Y,*,j,+^(1) 1 _G^(1)(*,j,+) _s\,:\,s_Y=j,\;s_Z>0 _G^(1)(s)· _G[s,1],Y. We then define ηG,Z(1)≔∑s:sX=0orsY=0αG(1)(s)⋅H(βG[s,1],Z)+∑k=02ℓ∗αG(1)(+,+,k)⋅H(β¯G,Z,+,+,k(1)), _G,Z^(1) _s\,:\,s_X=0\ or\ s_Y=0 _G^(1)(s)· H ( _G[s,1],Z )+ _k=0^2 ^* _G^(1)(+,+,k)· H ( β_G,Z,+,+,k^(1) ), where αG(1)(+,+,k)≔∑s:sX>0,sY>0,sZ=kαG(1)(s),β¯G,Z,+,+,k(1)≔1αG(1)(+,+,k)∑s:sX>0,sY>0,sZ=kαG(1)(s)⋅βG[s,1],Z. _G^(1)(+,+,k) _s\,:\,s_X>0,\;s_Y>0,\;s_Z=k _G^(1)(s), β_G,Z,+,+,k^(1) 1 _G^(1)(+,+,k) _s\,:\,s_X>0,\;s_Y>0,\;s_Z=k _G^(1)(s)· _G[s,1],Z. To avoid the issue of denominators αG(1)(⋅) _G^(1)(·) being zero, we regard 0⋅undefined≔00· undefined 0, so that Section 2.3 and Section 2.3 are always well-defined. The quantity EG(1)E_G^(1) is then given by EG(1)≔minH((αG(1))X)−Pℓ∗(αG(1)),H(β¯G,Y,*,*,*(1))−ηG,Y(1),H(β¯G,Z,*,*,*(1))−ηG,Z(1),E_G^(1) \H ( ( _G^(1) )_X )-P_S_ ^* ( _G^(1) ),\;H ( β_G,Y,*,*,*^(1) )- _G,Y^(1),\;H ( β_G,Z,*,*,*^(1) )- _G,Z^(1) \, where β¯G,W,*,*,*(1)≔∑s∈ℓ∗αG(1)(s)⋅βG[s,1],W,∀W∈X,Y,Z. β_G,W,*,*,*^(1) _s _ ^* _G^(1)(s)· _G[s,1],W, ∀\,W∈\X,Y,Z\. For r=2,…,6r=2,…,6, we define EG(r)E_G^(r) by applying the same formula Section 2.3 after relabeling the coordinates X,Y,ZX,Y,Z as πr(X),πr(Y),πr(Z) _r(X), _r(Y), _r(Z), and replacing the region index 11 with r. Finally, we set EG≔∑r=16AG(r)⋅EG(r).E_G _r=1^6A_G^(r)· E_G^(r). Retained exponents at level ℓ≥3 ≥ 3. For each positive-shape node T at level ℓ≥3 ≥ 3, we will calculate intermediate quantities ET,W(r)E_T,W^(r) for each region r∈[6]r∈[6] and W∈X,Y,ZW∈\X,Y,Z\. These quantities will later be aggregated to form the retained exponent of level ℓ . Note that the region index r in the calculation can be different from the region index rTr_T of node T itself. As in the above, we start by defining the quantities for region r=1r=1 where π1 _1 is the identity permutation. In the summations below, u ranges over Split(sT) (s_T). We first define ηT,Y(1) _T,Y^(1) ≔∑u:uZ=0(αT(1)(u)+αT(1)(sT−u))⋅H(βT[u,1],Y) _u\,:\,u_Z=0 ( _T^(1)(u)+ _T^(1)(s_T\!-\!u) )· H ( _T[u,1],Y ) +∑j=0minsT,Y, 2ℓ−1(αT(1)(*,j,+)+αT(1)(*,sT,Y−j,<))⋅H(β¯T,Y,*,j,+(1)), + _j=0 \s_T,Y,\,2 -1\ ( _T^(1)(*,j,+)+ _T^(1)(*,\,s_T,Y\!-\!j,\,<) )· H ( β_T,Y,*,j,+^(1) ), where αT(1)(*,j,+)≔∑u:uY=j,uZ>0αT(1)(u),αT(1)(*,sT,Y−j,<)≔∑u:uY=sT,Y−j,uZ<sT,ZαT(1)(u), _T^(1)(*,j,+) _u\,:\,u_Y=j,\;u_Z>0 _T^(1)(u), _T^(1)(*,\,s_T,Y\!-\!j,\,<) _ subarraycu\,:\,u_Y=s_T,Y-j,\\ u_Z<s_T,Z subarray _T^(1)(u), β¯T,Y,*,j,+(1)≔1αT(1)(*,j,+)+αT(1)(*,sT,Y−j,<)⋅∑u:uY=j,uZ>0(αT(1)(u)+αT(1)(sT−u))⋅βT[u,1],Y. β_T,Y,*,j,+^(1) 1 _T^(1)(*,j,+)+ _T^(1)(*,\,s_T,Y\!-\!j,\,<)· _u\,:\,u_Y=j,\;u_Z>0 ( _T^(1)(u)+ _T^(1)(s_T\!-\!u) )· _T[u,1],Y. We also define ηT,Z(1) _T,Z^(1) ≔∑u:uX=0oruY=0(αT(1)(u)+αT(1)(sT−u))⋅H(βT[u,1],Z) _u\,:\,u_X=0\ or\ u_Y=0 ( _T^(1)(u)+ _T^(1)(s_T-u) )· H ( _T[u,1],Z ) +∑k=0minsT,Z, 2ℓ−1(αT(1)(+,+,k)+αT(1)(<,<,sT,Z−k))⋅H(β¯T,Z,+,+,k(1)), + _k=0 \s_T,Z,\,2 -1\ ( _T^(1)(+,+,k)+ _T^(1)(<,\,<,\,s_T,Z\!-\!k) )· H ( β_T,Z,+,+,k^(1) ), where αT(1)(+,+,k)≔∑u:uX>0,uY>0,uZ=kαT(1)(u),αT(1)(<,<,sT,Z−k)≔∑u:uX<sT,X,uY<sT,Y,uZ=sT,Z−kαT(1)(u), _T^(1)(+,+,k) _u\,:\,u_X>0,\;u_Y>0,\;u_Z=k _T^(1)(u), _T^(1)(<,\,<,\,s_T,Z\!-\!k) _ subarraycu\,:\,u_X<s_T,X,\;u_Y<s_T,Y,\\ u_Z=s_T,Z-k subarray _T^(1)(u), β¯T,Z,+,+,k(1)≔1αT(1)(+,+,k)+αT(1)(<,<,sT,Z−k)⋅∑u:uX>0,uY>0,uZ=k(αT(1)(u)+αT(1)(sT−u))⋅βT[u,1],Z. β_T,Z,+,+,k^(1) 1 _T^(1)(+,+,k)+ _T^(1)(<,\,<,\,s_T,Z\!-\!k)· _u\,:\,u_X>0,\;u_Y>0,\;u_Z=k ( _T^(1)(u)+ _T^(1)(s_T\!-\!u) )· _T[u,1],Z. Again, as we define 0⋅undefined≔00· undefined 0, Section 2.3 and Section 2.3 are well-defined even when some denominators αT(1) _T^(1) are zeros. Then, ET,W(1)E^(1)_T,W are given by ET,X(1) E_T,X^(1) ≔mTAT(1)⋅(H((αT(1))X)−PSplit(sT)(αT(1))), m_TA_T^(1)· (H ( ( _T^(1) )_X )-P_ (s_T) ( _T^(1) ) ), ET,Y(1) E_T,Y^(1) ≔mTAT(1)⋅(H(βT,Y(1))−ηT,Y(1)), m_TA_T^(1)· (H ( _T,Y^(1) )- _T,Y^(1) ), ET,Z(1) E_T,Z^(1) ≔mTAT(1)⋅(H(βT,Z(1))−ηT,Z(1)). m_TA_T^(1)· (H ( _T,Z^(1) )- _T,Z^(1) ). For regions r=2,…,6r=2,…,6, we define ET,πr(X)(r)E_T, _r(X)^(r), ET,πr(Y)(r)E_T, _r(Y)^(r), and ET,πr(Z)(r)E_T, _r(Z)^(r) by applying Section 2.3, Section 2.3, and Section 2.3 after relabeling coordinates X,Y,ZX,Y,Z as πr(X),πr(Y),πr(Z) _r(X), _r(Y), _r(Z) and replacing region index 11 by r. To calculate the retained exponent for level ℓ∈[3,ℓ∗] ∈[3, ^*], we let ℓ+T_ ^+ denote the set of positive-shape nodes at level ℓ . Then, the retained exponent at level ℓ is given by Eℓ≔∑r=16min∑T∈ℓ+ET,X(r),∑T∈ℓ+ET,Y(r),∑T∈ℓ+ET,Z(r).E_ _r=1^6 \ _T _ ^+E_T,X^(r),\; _T _ ^+E_T,Y^(r),\; _T _ ^+E_T,Z^(r) \. Retained exponent at level 2. For a positive level-22 node T of shape (1,1,2)(1,1,2), we define (ET,X,ET,Y,ET,Z)≔mT⋅(1, 1,H(μT,μT,1−2μT)). (E_T,X,E_T,Y,E_T,Z ) m_T· (1,\;1,\;H ( _T, _T,1-2 _T ) ). Similarly, for a node T of shape (2,1,1)(2,1,1), we define (ET,X,ET,Y,ET,Z)≔mT⋅(H(μT,μT,1−2μT), 1, 1), (E_T,X,E_T,Y,E_T,Z ) m_T· (H ( _T, _T,1-2 _T ),\;1,\;1 ), and for a node T of shape (1,2,1)(1,2,1), we define (ET,X,ET,Y,ET,Z)≔mT⋅(1,H(μT,μT,1−2μT), 1), (E_T,X,E_T,Y,E_T,Z ) m_T· (1,\;H ( _T, _T,1-2 _T ),\;1 ), Then, letting 2+T_2^+ denote the set of positive-shape nodes at level 22, the retained exponent at level 2 is given by E2≔min∑T∈2+ET,X,∑T∈2+ET,Y,∑T∈2+ET,Z.E_2 \ _T _2^+E_T,X,\; _T _2^+E_T,Y,\; _T _2^+E_T,Z \. Local matrix size for zero-shape nodes. For every zero-shape node T, we let W0∈X,Y,ZW_0∈\X,Y,Z\ be the first zero coordinate of sTs_T and W1W_1 be the first nonzero coordinate of sTs_T. Then, we define the local matrix size of T, written (MT,X,MT,Y,MT,Z)(M_T,X,M_T,Y,M_T,Z), as MT,W0≔mT⋅(H(βT,W1)+∑L∈supp(βT,W1)βT,W1(L)⋅|p∈[2ℓ−1]|Lp=1|⋅logq);M_T,W_0 m_T· (H ( _T,W_1 )+ _L∈ ( _T,W_1) _T,W_1(L)· | \p∈[2 -1]\; |\;L_p=1 \ |· q ); the other two coordinates in the local matrix size are zeros. Local matrix size for (1,1,2)(1,1,2)-nodes. For a node T with shape (1,1,2)(1,1,2), we define its local matrix size as (MT,X,MT,Y,MT,Z)≔mT⋅((1−2μT)logq,(1−2μT)logq, 2μTlogq). (M_T,X,M_T,Y,M_T,Z ) m_T· ((1-2 _T) q,\;(1-2 _T) q,\;2 _T q ). Similarly, for a node T with shape (2,1,1)(2,1,1), we define its local matrix size as (MT,X,MT,Y,MT,Z)≔mT⋅(2μTlogq,(1−2μT)logq,(1−2μT)logq), (M_T,X,M_T,Y,M_T,Z ) m_T· (2 _T q,\;(1-2 _T) q,\;(1-2 _T) q ), and for a node T with shape (1,2,1)(1,2,1), we define its local matrix size as (MT,X,MT,Y,MT,Z)≔mT⋅((1−2μT)logq, 2μTlogq,(1−2μT)logq). (M_T,X,M_T,Y,M_T,Z ) m_T· ((1-2 _T) q,\;2 _T q,\;(1-2 _T) q ). 2.4 Final assembly Finally, we define the total retained exponent as Etotal≔EG+E2+∑ℓ=3ℓ∗Eℓ,E_total E_G+E_2+ _ =3 ^*E_ , and define the total matrix size as Mtotal≔min∑T∈ℒMT,X,∑T∈ℒMT,Y,∑T∈ℒMT,Z,M_total \ _T M_T,X,\; _T M_T,Y,\; _T M_T,Z \, where ℒL is the set of leaves. We can now write the full optimization problem. minimizeΩsubject toEtotal+Mtotal⋅Ω≥2ℓ∗−1log(q+2),all free variables lie in the domains stated in section 2.2. array[]lminimize& \\[1.99997pt] subject to&E_total+M_total· ≥ 2 ^*-1 (q + 2),\\ &all free variables lie in the domains stated in @crefcreftype @tilde refnumsubsec:optimizable-parameters. array (11) Theorem 1 (1). Any feasible solution of Equation 11 implies ω≤Ωω≤ , where ω is the asymptotic matrix multiplication exponent. 2.5 Dealing with maximum entropies Among the quantities introduced above, HDmax(ρ)H_D (ρ) is the only type that cannot be computed analytically. However, it is easy to see that replacing each occurrence of HDmaxH_D in the computation of EtotalE_total with its upper bound can only decrease EtotalE_total, so any solution that satisfies the constraints in Equation 11 after this replacement is feasible for Equation 11 itself, and still yields a valid upper bound ω≤Ωω≤ . In this section, we describe how to derive an upper bound of HDmax(ρ)H_D (ρ). For a set D of shapes and a probability distribution ρ∈Δ(D)ρ∈ (D), a valid certificate for HDmax(ρ)H_D (ρ) consists of: 1. a distribution y∈Δ(D)y∈ (D) satisfying y(a)>0y(a)>0 for every a∈Da∈ D and yW=ρWy_W= _W for every W∈X,Y,ZW∈\X,Y,Z\; 2. a real number λ0 _0 and, for each W∈X,Y,ZW∈\X,Y,Z\ and each marginal value w∈aW|a∈Dw∈\a_W\;|\;a∈ D\, a real number λW(w) _W(w) (these values play the role of Lagrange multipliers for the maximization defining HDmax(ρ)H_D (ρ)); and 3. a real number ε≥0 ≥ 0 satisfying, for every a=(aX,aY,aZ)∈Da=(a_X,a_Y,a_Z)∈ D, |logy(a)−(λ0+λX(aX)+λY(aY)+λZ(aZ))|≤ε. y(a)- ( _0+ _X(a_X)+ _Y(a_Y)+ _Z(a_Z) ) ≤ . Lemma 1. For any valid certificate, we have H(y)≤HDmax(ρ)≤H(y)+2εH(y)≤ H_D (ρ)≤ H(y)+2 . Proof. The left inequality follows from the definition: y has the same marginals as ρ, so it is feasible for the maximization defining HDmax(ρ)H_D (ρ) (see Section 2.3). For the right inequality, fix any ρ′∈Δ(D)ρ ∈ (D) with ρW′=ρWρ _W= _W for all W∈X,Y,ZW∈\X,Y,Z\; to prove the right inequality we need only show H(ρ′)≤H(y)+2εH(ρ )≤ H(y)+2 . The entropy function H is concave on Δ(D) (D), and since y is strictly positive, H is differentiable at y with ∇H(y)a=−logy(a)−loge∇ H(y)_a=- y(a)- e (recall that all logarithms are in base 2). Therefore H(ρ′)≤H(y)+⟨∇H(y),ρ′−y⟩.H(ρ )≤ H(y)+ ∇ H(y),\,ρ -y . Write g(a)≔λ0+λX(aX)+λY(aY)+λZ(aZ)g(a) _0+ _X(a_X)+ _Y(a_Y)+ _Z(a_Z), and view g, logy y, and the constant loge e as vectors indexed by a∈Da∈ D. Then, ⟨g+loge,ρ′−y⟩ g+ e,\,ρ -y =(λ0+loge)⋅∑a∈D(ρ′(a)−y(a))+∑W∈X,Y,Z∑a∈DλW(aW)⋅(ρ′(a)−y(a)) = ( _0+ e )· _a∈ D (ρ (a)-y(a) )+ _W∈\X,Y,Z\\, _a∈ D _W(a_W)· (ρ (a)-y(a) ) =(λ0+loge)⋅(1−1)+∑W∈X,Y,Z∑w∈aW∣a∈DλW(w)⋅(ρW′(w)−yW(w))=0, = ( _0+ e )·(1-1)+ _W∈\X,Y,Z\\, _w∈\a_W a∈ D\ _W(w)· (ρ _W(w)-y_W(w) )=0, where the last step uses ρW′=ρW=yWρ _W= _W=y_W for each W∈X,Y,ZW∈\X,Y,Z\. Hence, ⟨∇H(y),ρ′−y⟩=⟨g−logy,ρ′−y⟩≤ε⋅‖ρ′−y‖1≤2ε, ∇ H(y),\,ρ -y = g- y,\,ρ -y ≤ · ρ' - y_1≤ 2 , where the first inequality uses Item 3 and the second one uses the fact that the ℓ1 _1-distance between two probability distributions is at most 22. Thus, H(ρ′)≤H(y)+2εH(ρ )≤ H(y)+2 , proving the lemma. ∎ 3 Solving the optimization problem numerically We wish to numerically minimize Ω from Equation 11 while satisfying all the constraints, with q=5q=5 and ℓ∗=4 ^*=4. In 1, this was done (for ℓ∗=3 ^*=3) via a sequential quadratic programming (SQP) algorithm using the software package SNOPT (12). In contrast, we take a gradient-based approach to tackle the non-convex minimization. We develop an algorithm applying several techniques from machine learning and adjacent areas, such as optimal transport. 3.1 Differentiable objective To apply a gradient-based optimization algorithm, we need a differentiable objective function to be optimized. However, obtaining a differentiable objective from Equation 11 is challenging due to the problem constraints. Since most of the parameters to be optimized are distributions, this creates the constraints that the probabilities must be non-negative and must add up to one. To sidestep that issue and obtain unconstrained optimization, we parameterize distributions in terms of their logits, and obtain the distribution by applying the softmax operation on the free parameters. We initialize the algorithm with random (logit) parameters. The other relevant set of constraints comes from the maximum entropy distributions. Given an input probability distribution, we must find the distribution with maximum entropy that shares the same marginals as the input. This is reminiscent of a problem that typically arises in optimal transport, where the Sinkhorn-Knopp algorithm (19) is used to obtain a stable and differentiable objective (9). Thus, unlike 1, we do not treat the maximum entropy distributions (and the corresponding Lagrange multipliers) as free parameters to be optimized together with the rest of distributions; rather, we obtain the maximum entropy distributions (and the Lagrange multipliers) using the Sinkhorn-Knopp algorithm. We compute the gradients of the objective using automatic differentiation. To improve stability, we employ implicit differentiation for backpropagating through the Sinkhorn-Knopp algorithm (8; 11). We update the parameters using Adam (13). 3.2 Software implementation To obtain a fast algorithm that can leverage hardware platforms (such as GPUs) and scale up to level ℓ∗=4 ^*=4 (with nearly 7 million parameters to optimize) we implement our algorithm in Jax (5), applying some techniques to improve efficiency. The key to an effective speed-up is finding a new representation for the problem, switching from a loop over graph nodes (as done in the SQP implementation from 1) to a naturally parallelizable computation over tensors. In the implementation from 1, the free parameters are stored in the nodes of a graph where updates happen through message-passing (both from child nodes to parent nodes and vice versa). The graph is non-uniform enough (e.g., different nodes have different number of children and different types of them) to make parallelization hard, but we solve this via two techniques. Firstly, we introduce phantom nodes in the graph with masking; as a trade-off, we pay the cost of increasing the number of optimization parameters by up to a factor of 33. Secondly, we cluster the nodes into a limited number of highly specialized groups (“stages”). With these techniques, we are able to represent the full graph with multi-dimensional tensors, enabling parallel processing over a number of axes (up to 1010), challenging the limits of tensor-processing backends. 3.3 Applying AlphaEvolve We use AlphaEvolve (15) to further improve the optimization algorithm. Specifically, we let AlphaEvolve modify the optimization program, which is then executed (taking approximately 5 hours on a single GPU) to output a bound on omega. AlphaEvolve then evolves the code to minimize omega. We found improved results by using AlphaEvolve’s “evolving constructions” feature, where the optimization algorithm at each generation starts at the best solution point found by the parent algorithm. 4 Rigorous verification of the omega upper bound To rigorously certify our omega bound, we run a separate verification step computing all quantities in rational arithmetic to guard against floating point errors. More specifically, we round the floating point solution obtained at the end of the optimization to rational numbers, ensuring the maximum-entropy certificates stay valid. We then evaluate all derived quantities in exact rational arithmetic, and replace each logarithm with a rational bound rounded in the proper direction that ensures each constraint in Equation 11 is satisfied, so the certified bounds are free from numerical errors. We are preparing a repository in which we will release the verification code and our discovered solution. 5 Discussion In this note, we improve the SOTA upper bound on ω to 2.371177. Our improvement is comparable in magnitude to most improvements in the last 40 years since ω<2.376ω<2.376 was attained by 7. We achieve this by leveraging modern optimization techniques and AlphaEvolve to design a better optimization algorithm for the problem described by 1. While further modest improvements may be obtained in this manner, achieving larger improvements to ω likely requires new mathematical ideas and is an exciting area of research. Acknowledgments. The work presented in this note was motivated by initial conversations during the Complexity and Linear Algebra program at Simons Institute for the Theory of Computing in Fall 2025. References Alman et al. (2025) J. Alman, R. Duan, V. V. Williams, Y. Xu, Z. Xu, and R. Zhou More asymmetry yields faster matrix multiplication. In Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), p. 2005–2039. External Links: Document, Link, https://epubs.siam.org/doi/pdf/10.1137/1.9781611978322.63 Cited by: Table 1, §1, §1, §1, §2, §3.1, §3.2, §3, §5, Theorem 1, Abstract. Alman and Williams (2024) J. Alman and V. V. Williams A refined laser method and faster matrix multiplication. TheoretiCS 3 (11261). Cited by: §1. Bini et al. (1979) D. Bini, M. Capovani, F. Romani, and G. Lotti (n2.7799)O(n^2.7799) Complexity for n×n× n approximate matrix multiplication. Information Processing Letters 8 (5), p. 234–235. External Links: ISSN 0020-0190 Cited by: §1. Bläser (2013) M. Bläser Fast matrix multiplication. Theory of Computing, p. 1–60. Cited by: §1. Bradbury et al. (2018) J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang JAX: composable transformations of Python+NumPy programs. External Links: Link Cited by: §1, §3.2. Coppersmith and Winograd (1981) D. Coppersmith and S. Winograd On the asymptotic complexity of matrix multiplication. In 22nd Annual Symposium on Foundations of Computer Science (sfcs 1981), Vol. , p. 82–90. Cited by: §1. Coppersmith and Winograd (1990) D. Coppersmith and S. Winograd Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation 9 (3), p. 251–280. Note: Computational algebraic complexity editorial External Links: ISSN 0747-7171 Cited by: §1, §2, §5. Cuturi et al. (2020) M. Cuturi, O. Teboul, J. Niles-Weed, and J. Vert Supervised quantile normalization for low-rank matrix approximation. In International Conference on Machine Learning, Cited by: §3.1. Cuturi (2013) M. Cuturi Sinkhorn distances: lightspeed computation of optimal transport. In International Conference on Neural Information Processing Systems, p. 2292–2300. Cited by: §3.1. Duan et al. (2023) R. Duan, H. Wu, and R. Zhou Faster matrix multiplication via asymmetric hashing. In 2023 IEEE 64th Annual Symposium on Foundations of Computer Science (FOCS), Vol. , p. 2129–2138. External Links: Document Cited by: Table 1, §1, Abstract. Eisenberger et al. (2022) M. Eisenberger, A. Toker, L. Leal-Taixé, F. Bernard, and D. Cremers A unified framework for implicit Sinkhorn differentiation. In IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 499–508. Cited by: §3.1. Gill et al. (2002) P. E. Gill, W. Murray, and M. A. Saunders SNOPT: an SQP algorithm for large-scale constrained optimization. SIAM Journal on Optimization 12 (4), p. 979–1006. Cited by: §3. Kingma and Ba (2015) D. P. Kingma and J. Ba Adam: a method for stochastic optimization. In International Conference on Learning Representations (ICLR), Cited by: §3.1. Le Gall (2014) F. Le Gall Algebraic complexity theory and matrix multiplication. In Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation, ISSAC ’14, New York, NY, USA, p. 23. External Links: ISBN 9781450325011 Cited by: §1. Novikov et al. (2025) A. Novikov, N. Vũ, M. Eisenberger, E. Dupont, P. Huang, A. Z. Wagner, S. Shirobokov, B. Kozlovskii, F. J. R. Ruiz, A. Mehrabian, M. P. Kumar, A. See, S. Chaudhuri, G. Holland, A. Davies, S. Nowozin, P. Kohli, and M. Balog AlphaEvolve: a coding agent for scientific and algorithmic discovery. arXiv. External Links: 2506.13131, Link Cited by: §1, §3.3. Pan (1978) V. Y. Pan Strassen’s algorithm is not optimal trilinear technique of aggregating, uniting and canceling for constructing fast algorithms for matrix operations. In 19th Annual Symposium on Foundations of Computer Science, Vol. , p. 166–176. Cited by: §1. Romani (1982) F. Romani Some properties of disjoint sums of tensors related to matrix multiplication. SIAM Journal on Computing 11 (2), p. 263–267. Cited by: §1. Schönhage (1981) A. Schönhage Partial and total matrix multiplication. SIAM Journal on Computing 10 (3), p. 434–455. Cited by: §1. Sinkhorn and Knopp (1967) R. Sinkhorn and P. Knopp Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematics 21, p. 343–348. Cited by: §3.1. Stothers (2010) A. J. Stothers On the complexity of matrix multiplication. PhD thesis, University of Edinburgh, Edinburgh, United Kingdom. Cited by: §1. Strassen (1986) V. Strassen The asymptotic spectrum of tensors and the exponent of matrix multiplication. In 27th Annual Symposium on Foundations of Computer Science, Vol. , p. 49–54. Cited by: §1. Strassen (1969) V. Strassen Gaussian elimination is not optimal. Numerische mathematik 13 (4), p. 354–356. Cited by: §1. Williams et al. (2024) V. V. Williams, Y. Xu, Z. Xu, and R. Zhou New bounds for matrix multiplication: from alpha to omega. In Proceedings of the 2024 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), p. 3792–3835. External Links: Document, Link, https://epubs.siam.org/doi/pdf/10.1137/1.9781611977912.134 Cited by: Table 1, §1, Abstract. Williams (2012) V. V. Williams Multiplying matrices faster than Coppersmith-Winograd. In Proceedings of the Forty-Fourth Annual ACM Symposium on Theory of Computing, STOC ’12, New York, NY, USA, p. 887–898. External Links: ISBN 9781450312455 Cited by: §1.