Introduction

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,"¹ 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 \times n$ matrices $A, B$ over some field $\mathbb F$, we want to compute $C_{n\times n}$ s.t.

$C_{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 $C$? The authors of this paper posit that there are tight bounds on the number of operations:

$n^{\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.²

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

$\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 $n^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 $n^3$ approach) was optimal. This discovery was published in a paper aptly titled Gaussian elimination is not optimal.² 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: $(T_{ijk})_{ijk} \in \mathbb F^{n\times m\times p}$
• A bilinear map, taking two vectors to a third vector: $T: \mathbb F^n \times \mathbb F^m \rarr \mathbb F^p$
• A trilinear map, taking three vectors to a scalar: $T: \mathbb F^n \times \mathbb F^m \times \mathbb F^p \rarr c$
• a multilinear polynomial: $T: \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 $T$ as a multilinear polynomial over sets of formal variables $X = \{x_1, ..., x_n\}, Y = ..., Z = ...$ so that

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

where $\alpha_{ijk}$ are elements of some field $\mathbb F$ e.g. $\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 $\mathbf{x}, \mathbf{y}$, and every $k \in [p]$ (the number of polynomials we want to compute), we want to compute the coefficients of $z_k$ which is a bilinear polynomial:

$\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,t \in \mathbb N$, define:

$\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 $A_{m\times s}, B_{s\times t}$.

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

$\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:

$\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, $Rank$ is the complexity of a tensor, and a tensor is just the natural generalization of a matrix. The total $Rank$ of a tensor $T$ is the minimum number of $Rank$ one tensors whose product equals $T$. We define a tensor $Rank$ of one if there exists $a_i, b_j, c_k \in \mathbb F$ such that our tensor can be expressed as the product of linear combinations:

$T = \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 $T$ as $a \otimes b \otimes c$ (read "$a$ tensor $b$ tensor $c$"). Matrix $Rank$ is the same just without the third $z_k, c_k$ dimension.

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

At a high level, the decomposition of $T_n$ into $R$ rank one terms provides an algorithm for multiplying arbitrary $n \times n$ matrices using $R$ scalar multiplications.

Properties of Tensor Rank

• Lemma: if $Rank(\langle m, s, t\rangle) = r$, then $\omega \leq 3 \log_{mst} r$ via the recursive application of $Rank(\langle m, m, m\rangle) = r$ implying $\omega \leq 3 \log_m r$. We can express this as the sum of $r$ terms that are linear combinations of the $x$ entries of our tensor times the linear combination of the $y$ entries ... and so on:

$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. $\langle m, s, r\rangle$ representing two $2\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 $\ell = 1, ..., r$ recursively compute

$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)$

2. For coefficients of $z_{ki}$, compute

$\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 $R$ controlling the asymptotic complexity of the algorithm: a matrix $M_{n\times n}$ can be multiplied with asymptotic complexity of $O(n^{\log_N R})$.

Thus, $P_\ell$ is a "block" and a comparatively cheap quantity to compute since it's just a summation with few multiplications compared to our naive $n^3$ approach. Additionally, $\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 $r$, hence why we are concerned with the asymptotic rank of the MMT, and also how Strassen got

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

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

$Rank(\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.⁵ 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 $Rank$ 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\times 2$ MM using only 7 multiplications	Simple, recursive approach
$\vdots$	$\vdots$	$\vdots$
Coppersmith, Winograd 1987	The $CW^{\otimes 2}$ tensor and its powers	Laser Method, introduced by Coppersmith, Winograd, and Strassen
$\vdots$	$\vdots$	$\vdots$
Stothers, 2010	$CW^{\otimes 4}$	Laser Method
Vassilevska Williams, 2011	$CW^{\otimes 8}$	Laser Method
Le Gall, 2014	$CW^{\otimes 32}$	Laser Method
Vassilevska Williams, 2020	$CW^{\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,⁶ Josh Alman proved the correctness of the Laser Method, showing that for tensors $T$ for which it can be applied, if any other method could achieve the optimum of $\omega = 2$ using $T$, 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 $P$ that we have an efficient algorithm for
2. Show how to reduce the problem of Matrix Multiplication to $P$

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 \times 2$ matrices:

$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:

$\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:

$\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 $P$, given terms

$\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:

$\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 $x_{33}, y_{33}$ terms) problem $P$. Our measure of success is removal of some inputs and outputs of $P$ to transform it into the Matrix Multiplication Problem:

$\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 $P$. But how can we cancel those inputs and outputs?

If we phrase $P$ in terms of a tensor that we have a known, good $Rank$ bound for such as the $2 \times 2$ MMT:

$(\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 $P$ where:

$\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, Z$ variable of our instance of $P$ to zero. Setting $x_{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):

$\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 $P$ is reducible to $\langle 2, 2, 2 \rangle$ via zeroing as "$\langle 2, 2, 2 \rangle \leq_{zo} P$" which does have a well-defined $Rank$ bound.

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

if $\langle m, s, t \rangle \leq_{zo} P$, then $\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:

$\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}$

where

$\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}$

where

$\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:

$\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 $T$.

Direct Summation

For some tensor i.e. $T = 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' = (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 $y_2$ variable. Conversely, a direct sum of two copies of $T$ are completely disjoint: they have no variables in common despite summing over the same set of variables:

$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:⁷

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

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

$\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^{\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, T$ over variable sets $\{X, Y, Z\}, \{X', Y', Z'\}$, respectively, where

$S = \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:

$S \otimes T$ over $X \times X'$, $Y \times Y'$, and $Z \times Z'$ is given by:

$\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:

$T^{\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: $\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 $T^{\otimes N}$ will contain $N$ 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(S \otimes T) \leq Rank(S)\cdot Rank(T)$, so $Rank(T^{\otimes N}) \leq Rank(T)^N$

and sometimes, $Rank(T^{\otimes N}) \lt \lt Rank(T)^N$.

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

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

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

$... \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 $N$^th 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, k$, then we can get a really good bound on $\omega$.

Laser Method

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

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

Special Sauce: Tensor Partitioning

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

$T = \sum^n_a \sum^m_b \sum^p_c \alpha_{abc} x_a y_b z_c$

and partition its variable sets $X = 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 \in [\ell]$, let the subtensor $T_{ijk}$ be

$T_{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 $\alpha_{abc} \cdot x_a y_b z_c$ with $x_a \in X_i, y_b \in Y_j, z_c \in Z_jk$.

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

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

And if all the subtensors $T_{ijk}$ are MMTs, $T$ 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,t$ distinct parts such that we have a direct sum of $1 \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 $q$:

$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+2$ which is optimal since that's the number of $x,y$ or $z$ variables that appear in it e.g. there are always precisely $q+2 \; x_0$ variables in it:

$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 $n^2$, we'd be getting an asymptotic rank of $n^2$, or $\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:

$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

$\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 $N$^th 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

$T$ is constant sized, and partitioned s.t. the $T_{ijk}$ subtensors are matrix products. We take the $N$^th power of $T$:

$T^{\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 $N$^th power, we get the $N$^th power of a sum of subtensors, which in turn gives a sum of $N$-tuples of $T_1$ to $T_N$, where each $T_i$ is one of the subtensors $T_{ijk}$. From the assumption that our subtensors are matrix products, we can conclude that $T^{\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 $T_i$. We do this by associating some distribution $\pi$ to each of our $T_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 $T_{ijk}$ in the tensor power:

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

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

$\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 $T^{\otimes N}$ to get a direct sum of subtensors $T_1 \otimes ... \otimes T_N$ which correspond to $\pi$. This is the Laser Method.

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

$\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 $L$ and $\omega$: $L\uparrow \omega \downarrow$. We want to maximize $L$ 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$:

$\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_\pi, n_\pi, t_\pi$ as functions of $\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 $T_{ijk}$ into a distribution over single variables.

Marginalizing $\pi$

A marginal corresponds to the fraction of times that a term of shape $T_{i**}$ appears in $T_1 \otimes ... \otimes T_N$ for a given index $i$. We define the "$x$" marginals of $\pi$ as

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

Recall that $T^{\otimes N}$ is a tensor over the $N$-tuples $X^N, Y^N, Z^N$, so its $X$ variables are partitioned by $X_{c_1} \times ... \times X_{c_N}$ for $c_1, ..., c_N \in [\ell]$ and each subtensor $T_1 \otimes ... \otimes T_N$ corresponding to $\pi$ uses the $X$ variables from a set

$X_{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. $T_1 = T_{c_1jk}$ for some $j,k$. So, every term in our product(s) has some marginals it obeys. The Laser Method zeroes out variables in $T^{\otimes N}$ that do not have marginals in $\pi_i$ that we desire, yielding a sum of terms $T_1 \otimes ... \otimes T_N$