projectsblogpscoaboutreading list

49 | 2.3728596

10 April, 2023 - 33 min read


I ask the reader to call to mind some of the great rivalries over the centuries: Achilles vs. Hector, Tesla vs. Edison, Ford vs. Ferrari, ... and of course Vassilevska Williams vs. Le Gall.

In the past decade, a handful of researchers have been duking it out over matrix multiplication algorithm performance and it seems like no one is paying attention?!?! Today we're gonna poast about this paper "A Refined Laser Method and Faster Matrix Multiplication,"1 by Josh Alman and Virginia Vassilevska Williams which shatters, okay well ... it doesn't really shatter the previous best method – but it makes a notable improvement which can likely be generalized to other classes of tensors to push the upper bound even lower.

The Problem of Matrix Multiplication

Given n×nn \times n matrices A,BA, B over some field F\mathbb F, we want to compute Cn×nC_{n\times n} s.t.

Cij=kAikBkjC_{ij} = \sum_k A_{ik} \cdot B_{kj}

but we wanna go fast. How many operations over the field (addition, multiplication, division) are needed to compute CC? The authors of this paper posit that there are tight bounds on the number of operations:

nω+O(1);2ω3n^{\omega + O(1)}; \quad 2 \leq \omega \leq 3

but if you've ever taken a linear algebra exam, we know that we can certainly do worse than cubic time to apply the standard inner-product algorithm, making several arithmetic mistakes along the way to compute the result matrix.

The Inner Product Algorithm

The first (and usually last method) most of us pedestrian matrix multipliers learn just takes the dot product of all the rows of the first matrix with the columns of the second matrix. This is fine, but slow, as it requires many multiplication operations which are more expensive than additions. I won't get into the weeds (perhaps my next muse...) since there's plenty of weeds at the current layer of abstraction, but suffice it to say that the additions become asymptotically irrelevant and multiplication at the ALU level typically requires several additions, whereas the inverse is not true.2

If we're just limited to pen and paper, and compelled to multiply two 2×22\times 2 matrices XX and YY, the approach we'd probably use would look something like this:

X×Y=[x11x12x21x22][y11y12y21y22]=[x11y11+x12y21x11y12+x12y22x21y11+x22y21x21y12+x22y22]\begin{aligned} X \times Y &= \begin{bmatrix} \color{blue} x_{11} & \color{blue}x_{12} \\ \color{blue}x_{21} & \color{blue}x_{22} \end{bmatrix} \begin{bmatrix} \color{orange} y_{11} & \color{orange}y_{12} \\ \color{orange}y_{21} & \color{orange}y_{22} \end{bmatrix} \\ &= \begin{bmatrix} \color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21} & \color{blue}x_{11}\color{orange}y_{12} \color{black}+ \color{blue}x_{12}\color{orange}y_{22} \\ \color{blue}x_{21}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{21} & \color{blue}x_{21}\color{orange}y_{12} \color{black}+ \color{blue}x_{22}\color{orange}y_{22} \end{bmatrix} \end{aligned}

which is n3=8n^3 = 8 multiplications. Meh.

History of the Problem

Coppersmith and Winograd are heavily cited as a "black box" for other computational performance benchmarks or building blocks in the field of computational theory since matrix multiplication is the foundation of pretty much all applied linear algebraic questions. Operators including matrix inversion, rank, determinant, lower-upper factorization/decomposition, eigen anything-ing, etc. Matrix multiplication also has applications in other fields such as graph theory, complexity theory, etc. so the motivation for wanting to go fast is well established.

Up until 1969, it was thought that it wasn't possible to compute the product of two matrices in faster-than-cubic time. In fact, the discovery of a faster method was made by Volker Strassen when trying to prove that standard Gaussian elimination (an n3n^3 approach) was optimal. This discovery was published in a paper aptly titled Gaussian elimination is not optimal.2 And this was a watershed moment for investigation into faster means of matrix multiplication using various methods:


Till the past decade or so, when only minor improvments on the order of hundredths of thousandths of operations are being shaved off. This might not seem like much, but for massive matrices like the ones you print out to help you on your final exam, these add up. Additionally, the academic heart strives to know the theoretic optimum. Ultimately, the question of Matrix Multiplication optimization boils down to one of asymptotic rank of a Matrix Multiplication Tensor.

The Matrix Multiplication Tensor

There are many representations or interpretations of what a tensor is:

  • a multidimensional array: (Tijk)ijkFn×m×p(T_{ijk})_{ijk} \in \mathbb F^{n\times m\times p}
  • A bilinear map, taking two vectors to a third vector: T:Fn×FmFpT: \mathbb F^n \times \mathbb F^m \rarr \mathbb F^p
  • A trilinear map, taking three vectors to a scalar: T:Fn×Fm×FpcT: \mathbb F^n \times \mathbb F^m \times \mathbb F^p \rarr c
  • a multilinear polynomial: T:Fn×Fm×FpFT: \mathbb F^n \times \mathbb F^m \times \mathbb F^p \rarr \mathbb F

The last interpretation is the one widely used in the Matrix Multiplication Problem since it's useful for encoding the notion of a "problem" within another problem.

We define a tensor TT as a multilinear polynomial over sets of formal variables X={x1,...,xn},Y=...,Z=...X = \{x_1, ..., x_n\}, Y = ..., Z = ... so that

T=i=1nj=1mk=1pαijkxiyjzk,T = \sum_{i=1}^n \sum_{j=1}^m \sum_{k=1}^p \alpha_{ijk} \cdot x_i y_j z_k,

where αijk\alpha_{ijk} are elements of some field F\mathbb F e.g. C,R,Q\mathbb {C, R, Q}. As mentioned before, every tensor under this interpretation (which is equivalent to any of the others) defines a computationl problem: for given vectors x,y\mathbf{x}, \mathbf{y}, and every k[p]k \in [p] (the number of polynomials we want to compute), we want to compute the coefficients of zkz_k which is a bilinear polynomial:

injmαijkxiyj\sum^n_i \sum^m_j \alpha_{ijk} \cdot x_i y_j

The Matrix Multiplication Tensor (MMT), then, is the one to solve our problem. For m,s,tNm,s,t \in \mathbb N, define:

m,s,t:=i=1mj=1sk=1txijyjkzki\langle m,s,t\rangle := \sum_{i=1}^m \sum_{j=1}^s \sum_{k=1}^t x_{ij} y_{jk} z_{ki}

to be the MMT for multiplying two matrices Am×s,Bs×tA_{m\times s}, B_{s\times t}.

In general, the coefficient of zkiz_{ki} in m,s,t\langle m,s,t\rangle is the (i,ki,k)th entry of

[x11x1sxm1xms]×[y11y1tys1yst]\begin{bmatrix} x_{11} & \cdots & x_{1s} \\ \vdots & \ddots & \vdots \\ x_{m1} & \cdots & x_{ms} \\ \end{bmatrix} \times \begin{bmatrix} y_{11} & \cdots & y_{1t} \\ \vdots & \ddots & \vdots \\ y_{s1} & \cdots & y_{st} \\ \end{bmatrix}

The Matrix Multiplication Problem

The formal problem that we encode as an MMT is the Matrix Multiplication Problem:

n,n,n:=i=1nj=1nk=1nxijyjkzki\langle n, n, n \rangle := \sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n x_{ij} y_{jk} z_{ki}

about which we care for the asymptotic rank.

Tensor Rank

In general, RankRank is the complexity of a tensor, and a tensor is just the natural generalization of a matrix. The total RankRank of a tensor TT is the minimum number of RankRank one tensors whose product equals TT. We define a tensor RankRank of one if there exists ai,bj,ckFa_i, b_j, c_k \in \mathbb F such that our tensor can be expressed as the product of linear combinations:

T=(naixi)(mbjyj)(pckzk)=injmkpaibjckour coeffxiyjzkT = \Big(\sum^n a_i x_i\Big) \Big(\sum^m b_j y_j\Big) \Big(\sum^p c_k z_k\Big) = \sum_i^n \sum_j^m \sum_k^p \underbrace{a_i b_j c_k}_{\text{our coeff}} \cdot x_i y_j z_k

We often express this tensor TT as abca \otimes b \otimes c (read "aa tensor bb tensor cc"). Matrix RankRank is the same just without the third zk,ckz_k, c_k dimension.

A Rank  r=0Rank \; r = 0 mathematical object is just a scalar, an r=1r = 1 object is a vector; r=2r = 2 an n×nn \times n matrix, and Rank3Rank \geq 3 objects are tensors.

At a high level, the decomposition of TnT_n into RR rank one terms provides an algorithm for multiplying arbitrary n×nn \times n matrices using RR scalar multiplications.

Properties of Tensor Rank

  • Lemma: if Rank(m,s,t)=rRank(\langle m, s, t\rangle) = r, then ω3logmstr\omega \leq 3 \log_{mst} r via the recursive application of Rank(m,m,m)=rRank(\langle m, m, m\rangle) = r implying ω3logmr\omega \leq 3 \log_m r. We can express this as the sum of rr terms that are linear combinations of the xx entries of our tensor times the linear combination of the yy entries ... and so on:
Rank(m,s,t)=r(i,jaijxij)(j,kbjkyjk)(k,ickizki)Rank(\langle m, s, t\rangle) = \sum_{\ell}^r \Big(\sum_{i', j} a_{\ell i'j}x_{i'j} \Big) \Big(\sum_{j', k} b_{\ell j'k}y_{j'k} \Big) \Big(\sum_{k, i} c_{\ell ki}z_{ki} \Big)

and we use this expression to develop an algorithm for constant-sized matrices e.g. m,s,r\langle m, s, r\rangle representing two 2×22\times 2 matrices, and we use it to recursively solve for the product of larger matrices whose entries are themselves matrices:

[[    ][    ][    ][    ]]×[[    ][    ][    ][    ]]\begin{bmatrix} \begin{bmatrix} \; & \; \end{bmatrix} & \begin{bmatrix} \; & \; \end{bmatrix} \\ \\ \begin{bmatrix} \; & \; \end{bmatrix} & \begin{bmatrix} \; & \; \end{bmatrix} \end{bmatrix} \times \begin{bmatrix} \begin{bmatrix} \; & \; \end{bmatrix} & \begin{bmatrix} \; & \; \end{bmatrix} \\ \\ \begin{bmatrix} \; & \; \end{bmatrix} & \begin{bmatrix} \; & \; \end{bmatrix} \end{bmatrix}

The algorithm is:

  1. For all =1,...,r\ell = 1, ..., r recursively compute
P=(i,jaijxij)(j,kbjkyjk)P_\ell = \Big(\sum_{i', j} a_{\ell i'j}x_{i'j} \Big) \Big(\sum_{j', k'} b_{\ell j'k'}y_{j'k'} \Big)
  1. For coefficients of zkiz_{ki}, compute
rckiP\sum_{\ell}^r c_{\ell k i} P_\ell

This algorithm can be used to multiply block matrices. Via recursive application, we can multiply matrices of arbitrary size, with the rank RR controlling the asymptotic complexity of the algorithm: a matrix Mn×nM_{n\times n} can be multiplied with asymptotic complexity of O(nlogNR)O(n^{\log_N R}).

Thus, PP_\ell is a "block" and a comparatively cheap quantity to compute since it's just a summation with few multiplications compared to our naive n3n^3 approach. Additionally, rP\sum_{\ell}^r P_\ell is cheap to compute since it's just another linear combination which can be computed in linear time. So, the whole runtime of this algorithm is dependent on the number of recursive calls to compute our blocks which is rr, hence why we are concerned with the asymptotic rank of the MMT, and also how Strassen got

Rank(2,2,2)7    ω=2.81Rank(\langle 2, 2, 2 \rangle) \leq 7 \implies \omega = 2.81

Similarly, in 1978 Pan achieved ω=2.80\omega = 2.80 via

Rank(70,70,70)143,640    ωlog70143,640=2.80Rank(\langle 70, 70, 70 \rangle) \leq 143,640 \implies \omega \leq \log_{70} 143,640 = 2.80

It was around this time that people stopped analyzing fixed size Matrix Multiplication Problems because it quickly becomes unwieldy to work with such large objects. So, the method for analysis changed. In 1987, Coppersmith and Winograd published a new method for designing Matrix Multiplication Algorithms.5 These methods –which have become more or less standard as far as I can tell– rely on two components:

  1. an Algebraic Identity: corresponding to a RankRank representation of some small tensor
  2. a Method for Analyzing the Identity: which is the method of extracting a Matrix Multiplication Algorithm from other tensors
Publication Algebraic Identity Analysis
Strassen, 1969 2×22\times 2 MM using only 7 multiplications Simple, recursive approach
\vdots \vdots \vdots
Coppersmith, Winograd 1987 The CW2CW^{\otimes 2} tensor and its powers Laser Method, introduced by Coppersmith, Winograd, and Strassen
\vdots \vdots \vdots
Stothers, 2010 CW4CW^{\otimes 4} Laser Method
Vassilevska Williams, 2011 CW8CW^{\otimes 8} Laser Method
Le Gall, 2014 CW32CW^{\otimes 32} Laser Method
Vassilevska Williams, 2020 CW32CW^{\otimes 32} Refined Laser Method

The Laser Method

The Laser Method is an indirect technique for extracting Matrix Multiplication Algorithms out of other objects. It was developed by Strassen, and optimized by Coppersmith and Winograd. Crucially, as recently as 2019,6 Josh Alman proved the correctness of the Laser Method, showing that for tensors TT for which it can be applied, if any other method could achieve the optimum of ω=2\omega = 2 using TT, then so can the Laser Method. Therefore, the Laser Method is optimal for a certain subset of tensors.

How does it work?

The Laser Method is a composition of reduction problems:

  1. Take some algebraic problem PP that we have an efficient algorithm for
  2. Show how to reduce the problem of Matrix Multiplication to PP

This is nearly identical to Strassen's original approach to proving sub-cubic time. E.g. given the algebriac problem instance of Matrix Multiplication for two 2×22 \times 2 matrices:

X=[x11x12x21x22],Y=[y11y12y21y22],X = \begin{bmatrix} \color{blue} x_{11} & \color{blue}x_{12} \\ \color{blue}x_{21} & \color{blue}x_{22} \end{bmatrix}, Y = \begin{bmatrix} \color{orange} y_{11} & \color{orange}y_{12} \\ \color{orange}y_{21} & \color{orange}y_{22} \end{bmatrix},

evaluate the four polynomials:

x11y11+x12y21,  x11y12+x12y22x21y11+x22y21,  x21y12+x22y22 \color{blue}x_{11}\color{orange}y_{11} \color{black} + \color{blue}x_{12}\color{orange}y_{21}, \; \color{blue}x_{11}\color{orange}y_{12} \color{black} + \color{blue}x_{12}\color{orange}y_{22} \\ \color{blue}x_{21}\color{orange}y_{11} \color{black} + \color{blue}x_{22}\color{orange}y_{21}, \; \color{blue}x_{21}\color{orange}y_{12} \color{black} + \color{blue}x_{22}\color{orange}y_{22}

Hell, we may as well throw them into the matrix box that they belong in:

[x11y11+x12y21x11y12+x12y22x21y11+x22y21x21y12+x22y22]\begin{bmatrix} \color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21} & \color{blue}x_{11}\color{orange}y_{12} \color{black}+ \color{blue}x_{12}\color{orange}y_{22} \\ \color{blue}x_{21}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{21} & \color{blue}x_{21}\color{orange}y_{12} \color{black}+ \color{blue}x_{22}\color{orange}y_{22} \end{bmatrix}

and suppose we can solve some other problem PP, given terms

x11,x12,x21,x22,x33,y11,y12,y21,y22,y33\color{blue}x_{11}, x_{12}, x_{21}, x_{22},\color{black} x_{33}, \\ \color{orange} y_{11}, y_{12}, y_{21}, y_{22}, \color{black} y_{33}

we want to compute the 5 polynomials:

x11y11+x12y21,x11y12+x12y22+x33y33,x21y11+x22y21,x21y12+x22y22,x11y11+x22y22\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21}, \\ \color{blue}x_{11}\color{orange}y_{12} \color{black}+ \color{blue}x_{12}\color{orange}y_{22} \color{black}+ \color{black} x_{33}y_{33}, \\ \color{blue}x_{21}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{21}, \\ \color{blue}x_{21}\color{orange}y_{12} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}, \\ \color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}

At a high level, our goal is to embed matrix multiplication into this similar, and notably larger (by the introduction of the x33,y33x_{33}, y_{33} terms) problem PP. Our measure of success is removal of some inputs and outputs of PP to transform it into the Matrix Multiplication Problem:

x11y11+x12y21,x11y12+x12y22+x33y33,x21y11+x22y21,x21y12+x22y22,x11y11+x22y22\color{blue}x_{11}\color{orange}y_{11} \color{black} + \color{blue}x_{12}\color{orange}y_{21}, \\ \color{blue}x_{11}\color{orange}y_{12} \color{black} + \color{blue}x_{12}\color{orange}y_{22} + \color{red} \cancel{\color{black}x_{33}y_{33}}, \\ \color{blue}x_{21}\color{orange}y_{11} \color{black} + \color{blue}x_{22}\color{orange}y_{21}, \\ \color{blue}x_{21}\color{orange}y_{12} \color{black} + \color{blue}x_{22}\color{orange}y_{22}, \\ \color{red}\cancel{\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}}

which would leave us with just the polynomials needed for Matrix Multiplication as part of PP. But how can we cancel those inputs and outputs?

If we phrase PP in terms of a tensor that we have a known, good RankRank bound for such as the 2×22 \times 2 MMT:

(x11y11+x12y21)z11+(x11y12+x12y22)z12+(x21y11+x22y21)z21+(x21y12+x22y22)z22(\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21}\color{black})z_{11} \color{black}+ (\color{blue}x_{11}\color{orange}y_{12} \color{black}+ \color{blue}x_{12}\color{orange}y_{22}\color{black})z_{12} + \\ (\color{blue}x_{21}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{21}\color{black})z_{21} \color{black}+ (\color{blue}x_{21}\color{orange}y_{12} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}\color{black})z_{22}

we must now show how to reduce it to PP where:

P=(x11y11+x12y21)z11+(x11y12+x12y22+x33y33)z12+(x21y11+x22y21)z21+(x21y12+x22y22)z22+(x11y11+x22y22)z3\begin{aligned} P = &(\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21}\color{black})z_{11} + \\ &(\color{blue}x_{11}\color{orange}y_{12} \color{black}+ \color{blue}x_{12}\color{orange}y_{22} \color{black}+ \color{black} x_{33}y_{33})z_{12} + \\ &(\color{blue}x_{21}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{21}\color{black})z_{21} + \\ &(\color{blue}x_{21}\color{orange}y_{12} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}\color{black})z_{22} + \\ &(\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}\color{black})z_{3} \end{aligned}

The way we achieve this reduction is be setting some of the X,Y,ZX, Y, Z variable of our instance of PP to zero. Setting x33,y33,z30x_{33}, y_{33}, z_3 \leftarrow 0 yields our desired Matrix Multiplication tensor exactly (which is not always the case as we'll see shortly):

P=(x11y11+x12y21)z11+(x11y12+x12y22+x33y33)z12+(x21y11+x22y21)z21+(x21y12+x22y22)z22+(x11y11+x22y22)z3\begin{aligned} P = &(\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21}\color{black})z_{11} + \\ &(\color{blue}x_{11}\color{orange}y_{12} \color{black}+ \color{blue}x_{12}\color{orange}y_{22} \color{black} + \color{red} \cancel{\color{black} x_{33}y_{33}}\color{black})z_{12} + \\ &(\color{blue}x_{21}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{21}\color{black})z_{21} + \\ &(\color{blue}x_{21}\color{orange}y_{12} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}\color{black})z_{22} + \\ &(\color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{22}\color{orange}y_{22}\color{black})\color{red} \cancel{z_{3}} \end{aligned}

This is called zeroing-out a tensor which is similar to the notion of reducibility in complexity theory. We denote that PP is reducible to 2,2,2\langle 2, 2, 2 \rangle via zeroing as "2,2,2zoP\langle 2, 2, 2 \rangle \leq_{zo} P" which does have a well-defined RankRank bound.

Lemma: If QzoPQ \leq_{zo} P, then Rank(Q)Rank(P)Rank(Q) \leq Rank(P), so

if m,s,tzoP\langle m, s, t \rangle \leq_{zo} P, then ω3logmstRank(P)\omega \leq 3 \log_{mst} Rank(P).

Brief aside about the Strassen method

Strassen showed that it is possible to compute those polynomials using the following identities:

[x11x12x21x22][y11y12y21y22]=[m11m12m21m22]\begin{bmatrix} \color{blue} x_{11} & \color{blue}x_{12} \\ \color{blue}x_{21} & \color{blue}x_{22} \end{bmatrix} \begin{bmatrix} \color{orange} y_{11} & \color{orange}y_{12} \\ \color{orange}y_{21} & \color{orange}y_{22} \end{bmatrix} = \begin{bmatrix} m_{11} & m_{12} \\ m_{21} & m_{22} \end{bmatrix}


m11=x11y11+x12y21m21=w+v+(x11+x12x21x12)y22m21=w+u+x22(y21+y12y11y22)m22=w+u+v\begin{aligned} m_{11} &= \color{blue}x_{11}\color{orange}y_{11} \color{black}+ \color{blue}x_{12}\color{orange}y_{21} \color{black} \\ m_{21} &= w + v + (x_{11} + x_{12} - x_{21} - x_{12})y_{22} \\ m_{21} &= w + u + x_{22}(y_{21} + y_{12} - y_{11} - y_{22}) \\ m_{22} &= w + u + v \end{aligned}


u=(x21x11)(y12y22)v=(x21+x22)(y12y11)w=x11y11+(x21+x22x11)(y11+y22y12)\begin{aligned} u &= (x_{21} - x_{11})(y_{12} - y_{22}) \\ v &= (x_{21} + x_{22})(y_{12} - y_{11}) \\ w &= x_{11}y_{11} + (x_{21} + x_{22} - x_{11})(y_{11} + y_{22} - y_{12}) \\ \end{aligned}

which reduces the number of matrix additions from 18 to 15, and the multiplications from 8 to 7!

Laser Method

The Laser Method is a bit more complicated than just the above reduction, namely due to the fact that it becomes increasingly challenging to zero out terms without removing terms we still need for the Matrix Multiplication proper for larger input matrices.

Thus, the Laser Method instead reduces direct sums of Matrix Multiplication Tensors to powers of tensors with a special structure via the zero-out reduction illustrated above. The process looks like this:

n,n,n    imi,si,tiMM reducable to direct summation of other MMTs    laser methodTN    Tbounding T via its powers\underbrace{\langle n, n, n \rangle \implies \oplus_i \langle m_i, s_i, t_i \rangle}_{\text{MM reducable to direct summation of other MMTs}} \overset{\text{laser method}}{\implies} \underbrace{T^{\otimes N} \implies T}_{\text{bounding } T \text{ via its powers}}

and composing all of these gives us a reduction from Matrix Multiplication to a tensor TT.

Direct Summation

For some tensor i.e. T=x1y1z1+x1y2z2T = x_1y_1z_1 + x_1y_2z_2, we can add it with a slightly modified copy of itself yielding a "regular" sum:

A=T+T=(x1y1z1+x1y2z2)+(x2y2z2+x2y3z3)A = T + T' = (x_1y_1z_1 + x_1\color{red}y_2\color{black}z_2) + (x_2\color{red}y_2\color{black}z_2 + x_2y_3z_3)

Observe that both copies share a y2y_2 variable. Conversely, a direct sum of two copies of TT are completely disjoint: they have no variables in common despite summing over the same set of variables:

B=T+T=(x1y1z1+x1y2z2)+(x3y3z3+x3y4z4)B = T + T' = (x_1y_1z_1 + x_1y_2z_2) + (x_3y_3z_3 + x_3y_4z_4)

Schönhage showed that, because of this disjoint-ness, upper bounds for MMTs also yield upper bounds for ω\omega via the following theorem:7

If the direct sum of kk copies of m,s,t\langle m, s, t \rangle have RankrRank \leq r m then ω3logmstrk\omega \leq 3 \log_{mst} \frac{r}{k}.

This means that we just have to show that the rank of kk disjoint copies of some MMT is smaller than kk times the rank of the original tensor. This is the first step of the method diagrammed above:

n,n,nω3logmstrk    km,s,tRank=r    ...\underbrace{\langle n, n, n \rangle}_{\omega \leq 3 \log_{mst} \frac{r}{k} } \implies \underbrace{\oplus^k \langle m, s, t \rangle}_{Rank = r} \implies ...

Tensor Exponentiation

Now, for the other half of the method, what the hell does this mean: T    TT^{\otimes} \implies T. We need to define what it means to exponentiate a tensor. First, we define the product of two tensors via the Kronecker Product.

For tensors S,TS, T over variable sets {X,Y,Z},{X,Y,Z}\{X, Y, Z\}, \{X', Y', Z'\}, respectively, where

S=i,j,kαijkxiyjzk,T=i,j,kβijkxiyjzkS = \sum_{i,j,k} \alpha_{ijk}x_iy_jz_k, \\ T = \sum_{i',j',k'} \beta_{i'j'k'}x'_{i'}y'_{j'}z'_{k'}

not only does our notation converge on sanscrit, but we get a definition of a tensor product:

STS \otimes T over X×XX \times X', Y×YY \times Y', and Z×ZZ \times Z' is given by:

:SX×Y×Z×TX×Y×ZF,ST=i,j,k,  i,j,kαijkβijk(xixi)(yjyj)(zkzk)\otimes : S^{X \times Y \times Z} \times T^{X' \times Y' \times Z'} \rarr \mathbb F, \\ \\ S \otimes T = \sum_{\color{blue}i,j,k, \; \color{red}i',j',k'} \color{black}\alpha_{\color{blue}ijk}\beta_{\color{red}i'j'k'}\color{black}(\color{blue}x_i \color{red} x'_{i'}\color{black})(\color{blue}y_j \color{red} y'_{j'}\color{black})(\color{blue}z_k \color{red} z'_{k'}\color{black})

which we immediately generalize into the definition of the power of a tensor:

TN=TN copiesTT^{\otimes N} = T \otimes \underbrace{\cdots}_{N \text{ copies}} \otimes T

Lemma: The product of any two MMTs can be combined into a product of their dimensions: a1,a2,a3b1,b2,b3=a1b1,a2b2,a3b3\langle a_1, a_2, a_3 \rangle \otimes \langle b_1, b_2, b_3 \rangle = \langle a_1b_1, a_2b_2, a_3 b_3\rangle which precisely corresponds to our idea of "blocking."

Observe that TNT^{\otimes N} will contain NN tuples over each variable set of the original tensor, so we get larger and larger variable sets, but the rank is guaranteed to be bounded by the individual ranks of the tensors being multiplied per the following:

Lemma: Rank(ST)Rank(S)Rank(T)Rank(S \otimes T) \leq Rank(S)\cdot Rank(T), so Rank(TN)Rank(T)NRank(T^{\otimes N}) \leq Rank(T)^N

and sometimes, Rank(TN)<<Rank(T)NRank(T^{\otimes N}) \lt \lt Rank(T)^N.

From these lemmas, we can get an expression for asymptotic rank:

R=limNRank(TN)1/N\overrightharpoon{R} = \lim_{N \rarr \infty} Rank(T^{\otimes N})^{1/N}

For example, for the 2×22 \times 2 MMT, Rank(2,2,2)=7Rank(\langle 2, 2, 2\rangle) = 7, but the rank of the NNth power of it is nωn^\omega! Thus, asymptotic rank is the secret workhorse of all modern Matrix Multiplication algorithms. This completes the right half of our methodology diagram:

...    laser methodTNRrN    TRank=r... \overset{\text{laser method}}{\implies} \underbrace{T^{\otimes N}}_{\overrightharpoon{R} \leq r^N} \implies \underbrace{T}_{Rank = r}

We can interpret this to mean: if we have a bound on the asymptotic rank of the tensor, we also have a bound on the asymptotic rank of its NNth power. So now we just have to show the reduction via Laser Method allowing us to compose the left and right halves. If we can compose these steps for sufficiently large m,s,t,km, s, t, k, then we can get a really good bound on ω\omega.

Laser Method

Finally, the middle bit: ...    laser method...... \overset{\text{laser method}}{\implies} ...

The Laser Method takes a tensor TT with a "special structure" and uses that structure to reduce the direct sum m,s,t\oplus \langle m, s, t\rangle to TNT^{\otimes N} for large m,s,t,km,s,t, k functions of NN by zeroing out many variables of TNT^{\otimes N} yielding huge disjoint sums.

Special Sauce: Tensor Partitioning

What is the aforementioned "special structure?" We start with some tensor:

T=anbmcpαabcxaybzcT = \sum^n_a \sum^m_b \sum^p_c \alpha_{abc} x_a y_b z_c

and partition its variable sets X=X1...X,Y=Y1...Y,Z=Z1...ZX = X_1 \cup ... \cup X_\ell, Y = Y_1 \cup ... \cup Y_\ell, Z = Z_1 \cup ... \cup Z_\ell into \ell parts. Then, for i,j,k[]i,j,k \in [\ell], let the subtensor TijkT_{ijk} be

Ti,j,k=xaXiybYjzcZjkαabcxaybzcT_{i,j,k} = \sum_{x_a \in X_i} \sum_{y_b \in Y_j} \sum_{z_c \in Z_jk} \alpha_{abc} \cdot x_a y_b z_c

which contains only the triples αabcxaybzc\alpha_{abc} \cdot x_a y_b z_c with xaXi,ybYj,zcZjkx_a \in X_i, y_b \in Y_j, z_c \in Z_jk.

Our whole tensor TT can then be expressed as a sum of subtensors:

T=ijkTijkT = \sum^\ell_i \sum^\ell_j \sum^\ell_k T_{ijk}

And if all the subtensors TijkT_{ijk} are MMTs, TT is a sum of MMTs, though not yet a direct sum. We want them to be a disjoint, direct sum, which is provably possible for all such partitioned tensors. Trivially, we could achieve this by partitioning each variable set into m,s,tm,s,t distinct parts such that we have a direct sum of 1×1×11 \times 1 \times 1 Matrix Multiplication Tensors. Such a decomposition isn't very useful though.

Let's take a look instead at a useful partitioning over a nontrivial tensor such as the Coppersmith-Winograd tensor parameterizd over some natural number qq:

CWq=x0y0zq+1+x0yq+1+xq+1y0z0+iq(x0yizi+ziy0zi+xiyiz0)CW_q = x_0y_0z_{q+1} + x_0y_{q+1} + x_{q+1}y_0z_0 + \sum^q_i (x_0y_iz_i + z_iy_0z_i + x_iy_iz_0)

This tensor is fascinating because its asymptotic rank is q+2q+2 which is optimal since that's the number of x,yx,y or zz variables that appear in it e.g. there are always precisely q+2  x0q+2 \; x_0 variables in it:

CWq=x0+1y0zq+1+x0+1yq+1+xq+1y0z0+iq(x0+qyizi+ziy0zi+xiyiz0)CW_q = \color{blue}\underbrace{x_0}_{+1}\color{black}y_0z_{q+1} + \color{blue}\underbrace{x_0}_{+1}\color{black}y_{q+1} + x_{q+1}y_0z_0 + \sum^q_i (\color{blue}\underbrace{x_0}_{+q}\color{black}y_iz_i + z_iy_0z_i + x_iy_iz_0)

If we had an optimal asymptotic rank for an MMT, since the number of variables is n2n^2, we'd be getting an asymptotic rank of n2n^2, or ω=2\omega = 2, so – we'd really like to embed the Matrix Multiplication Problem in this tensor.

The optimal partitioning for this tensor resembles the following:

X0={x0},X1={x1,...,xq},X2={xq+1},Y0={y0},Y1={y1,...,yq},Y2={yq+1},Z0={z0},Z1={z1,...,zq},Z2={zq+1},X_0 = \{x_0\}, X_1 = \{x_1, ..., x_q\}, X_2 = \{x_{q+1}\}, \\ Y_0 = \{y_0\}, Y_1 = \{y_1, ..., y_q\}, Y_2 = \{y_{q+1}\}, \\ Z_0 = \{z_0\}, Z_1 = \{z_1, ..., z_q\}, Z_2 = \{z_{q+1}\}, \\

The tensor then becomes

CWq=T002+T020+T200+qT011+T101+T100=1,1,1+1,1,1+1,1,1scalars+1,1,q+q,1,1+1,q,1inner products\begin{aligned} CW_q &= T_{002} + T_{020} + T_{200} + \sum^q T_{011} + T_{101} + T_{100} \\ &= \underbrace{\langle 1,1,1 \rangle + \langle 1,1,1 \rangle + \langle 1,1,1 \rangle}_{\text{scalars}} + \underbrace{\langle 1,1,q \rangle + \langle q,1,1 \rangle + \langle 1,q,1 \rangle}_{\text{inner products}} \end{aligned}

The inner products are nontrivial MMTs, and therefore challenging to get any fruitful zeroations out of. Instead, the Laser Method doesn't zero out the tensor itself, but takes the NNth power of it to give us more freedom for zeroing out.

The image of this process in my head is: taking two squares (matrices) we want to multiply, transforming their representation into a cube (tensor), exponentiating that cube into a rectangular prism, and then scratching out a bunch of the numbers we don't need until we're left with something "close" to the number of elements of just two slices of the rectangular prism which would be our initial matrices.

Taking a Large Power

TT is constant sized, and partitioned s.t. the TijkT_{ijk} subtensors are matrix products. We take the NNth power of TT:

TN=(i,j,kTijk)N=(T1,...,TN){Tijki,j,k[]}NT1...TNT^{\otimes N} = \Big (\sum_{i,j,k} T_{ijk}\Big )^{\otimes N} = \sum_{(T_1, ..., T_N) \in \{T_{ijk} \vert i,j,k \in [\ell]\}^N} T_1 \otimes ... \otimes T_N

When we take the NNth power, we get the NNth power of a sum of subtensors, which in turn gives a sum of NN-tuples of T1T_1 to TNT_N, where each TiT_i is one of the subtensors TijkT_{ijk}. From the assumption that our subtensors are matrix products, we can conclude that TNT^{\otimes N} is a tensor product of matrix products, which we've established is also a matrix product. So, it's a sum of matrix products, just really really big ones. Recall that what we really want is a disjoint sum of matrix products – all with the same dimensions. Ours currently do not have the same dimension unless our subtensors all had the same dimension to begin with which was not a constraint of the partition.

So, we zero out some of our TiT_i. We do this by associating some distribution π\pi to each of our T1,...,TN{Tijk}NT_1, ..., T_N \in \{ T_{ijk}\}^N where π\pi is the frequency of each term, or the fraction of times we used a given subtensor TijkT_{ijk} in the tensor power:

πijk=1N{p[N]Tp=Tijk}\pi_{ijk} = \frac{1}{N} \Big \vert \{ p \in [N] \vert T_p = T_{ijk}\} \Big \vert

Now, all of our TiTNT_i \in T^{\otimes N} whose (T1,...,TN)(T_1, ..., T_N) NN-tuple terms correspond to the same π\pi are matrix products with the same dimension: mπ,nπ,tπ\langle m_\pi, n_\pi, t_\pi \rangle. This is because

a,b,ctp,q,ruisoatpu,btqu,ctru\langle a, b, c \rangle^{\otimes t} \otimes \langle p, q, r \rangle^{\otimes u} \overset{iso}{\equiv} \langle a^tp^u, b^tq^u, c^tr^u \rangle

the location of these terms in the sum does not matter, just the number of occurences – hence the ability to proxy dimension with our distribution ascription. So, we just need to pick a π\pi which allows us to zero out the most terms or matrix products of disparate dimension.

Goal: pick π\pi and zero our variables in TNT^{\otimes N} to get a direct sum of subtensors T1...TNT_1 \otimes ... \otimes T_N which correspond to π\pi. This is the Laser Method.

Suppose lasering thusly gives us LL such tensors in the direct summation – we can use Schönhage's theorem and get

ω3logRank(TN)Llogmπnπtπ;Lω\omega \leq 3 \frac{\log \frac{Rank(T^{\otimes N})}{L}}{\log m_\pi n_\pi t_\pi }; \quad L \propto \omega

yielding an inversely proportionate relationship between LL and ω\omega: LωL\uparrow \omega \downarrow. We want to maximize LL for a fixed distribution π\pi and express it in terms of our distribution. Then, we can minimize the upper bound on ω\omega over the choices for π\pi:

n,n,nωlimN3logmstrk    Scho¨nhagekm,s,tRankr    laser methodTNRankrN    recursionT=i,j,kTijkpartitioned tensor withRank=r\underbrace{\langle n, n, n \rangle}_{\omega \lim_{N \rightarrow \infty}\leq 3 \log_{mst} \frac{r}{k} } \overset{\text{Schönhage}}{\implies} \underbrace{\oplus^k \langle m, s, t \rangle}_{Rank \leq r} \overset{\text{laser method}}{\implies} \underbrace{T^{\otimes N}}_{Rank \leq r^N} \overset{\text{recursion}}{\implies} \underbrace{T=\sum_{i,j,k} T_{ijk}}_{\text{partitioned tensor with} Rank = r}

We're left to determmine L,mπ,nπ,tπL, m_\pi, n_\pi, t_\pi as functions of π,N,T\pi, N, T. And already we've got a problem: we're scheming to zero out variables but π\pi is currently a distribution over triples not individual variables. So, we have to do something to convert π\pi from triples of TijkT_{ijk} into a distribution over single variables.

Marginalizing π\pi

A marginal corresponds to the fraction of times that a term of shape TiT_{i**} appears in T1...TNT_1 \otimes ... \otimes T_N for a given index ii. We define the "xx" marginals of π\pi as

πi:=j,k[]πijk\pi_i :=\sum_{j,k \in [\ell]} \pi_{ijk}

Recall that TNT^{\otimes N} is a tensor over the NN-tuples XN,YN,ZNX^N, Y^N, Z^N, so its XX variables are partitioned by Xc1×...×XcNX_{c_1} \times ... \times X_{c_N} for c1,...,cN[]c_1, ..., c_N \in [\ell] and each subtensor T1...TNT_1 \otimes ... \otimes T_N corresponding to π\pi uses the XX variables from a set

Xc1×...×XcN where i[],{p[N]cp=i}=πiNX_{c_1} \times ... \times X_{c_N} \text{ where } \\ \forall i \in [\ell], \Big| \{ p \in [N] \vert c_p = i\}\Big| = \pi_i N

e.g. T1=Tc1jkT_1 = T_{c_1jk} for some j,kj,k. So, every term in our product(s) has some marginals it obeys. The Laser Method zeroes out variables in TNT^{\otimes N} that do not have marginals in πi\pi_i that we desire, yielding a sum of terms T1...TNT_1 \otimes ... \otimes T_N