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.

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:

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:

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

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:

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

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

The Matrix Multiplication Problem

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

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:

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:

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:

The algorithm is:

1. For all $\ell = 1, ..., r$ recursively compute

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

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

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

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

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:

evaluate the four polynomials:

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

and suppose we can solve some other problem $P$, given terms

we want to compute the 5 polynomials:

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:

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:

we must now show how to reduce it to $P$ where:

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

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:

where

where

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:

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:

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:

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:

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

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:

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

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:

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:

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:

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

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:

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

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:

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:

The tensor then becomes

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

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:

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

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

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

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

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

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$ that obey those marginals. All remaining variables have exactly the correct frequency for every block $T_i$. This effectively converts the sum to a direct sum. The dimension of this direct sum –the number of $X$ variable sets matching $\pi_i$– is the multinomial coefficient:

and is bounded by the number of variable sets we have by definition of a disjoint sum. So the number of disjoint subtensors is $L \leq M$, and we want to maximize the number of disjoint subtensors. To do this, we should try to find nontrivial partitions s.t. $L \rarr M$. The Laser Method achieves $L \approx M$ if and only if the linear system defined by the marginals has a unique solution:

Iff this linear system has a unique solution, the Laser Method guarantees that we'll get the maximum number of disjoint tensors: $L \approx M$. And, even if no solution exists –that is, the marginals of $\pi$ don't determine $\pi$– the Laser Method still always gives a disjoint sum of roughly $M$ "good" disjoint subtensors of $T^{\otimes N}$ which are consistent with $\pi$.

Additionally, if there are other distributions with the same marginals as $\pi$, then there may be "extra" subtensors also consistent with our distribution and which share variables. We have to remove them, which adds a lot of noise to the zeroing process.

E.g., suppose $T^{ijk}$ uses $X_i, Y_j, Z_k$; and $T^{iii}$ are "good" subtensors. The Laser Method gives us

Removing these noisy terms costs a sometimes/usually significant penalty since we have to remove $s\cdot M$ noisy/bad tensors via zeroing, which also removes some of our $M$ good subtensors. Previous approaches using the laser method yielded the following bound on the remaining viable tensors:

The latest development introduced by Vassilevska Williams shows that it's possible to attain the following improvement:

which impacts the overall $\omega$ bound as follows:

It's Hip to be Square (rooted)

So, how did the authors make a categorical improvement on $s$ bad tensors?

It must be super sophisticated and Peter I can only take so many more confounded summations!

Nay, dear reader.

Whereas the prior method:

1. greedily selected a good subtensor $T^{iii}$ which shared variables with $O(s)$ extra subtensors requiring removal
2. Zeroed out $O(s)$ variables not in the good tensors to remove $O(s)$ bad terms,
3. Each zeroed-out variable would zero out at most one good subtensor

thus, it would remove $O(s)$ good subtensors for each one it would keep, resulting in the $\Theta(M/s) = L$ good remaining subtensors.

The Refined Laser Method approach instead employs a probabilistic selection of candidate terms for removal:

1. Pick $\Theta(M/\sqrt s)$ good subtensors at random, zeroing out every variable they don't use

This would remove all the extra subtensors with non-zero probability.

B-b-but non-zero could be really bad

The analysis is optimal for large a family of tensors. This refined method improves upon the previous greedy approach whenever the marginals of $\pi$ do not uniquely determine $\pi$ which is most of the time for sufficiently large partitioned tensors, including $CW^{\otimes 8}$ and its powers. Since the previous best Matrix Multiplication algorithms use $CW^{\otimes 32}$, this immediately improves $\omega \leftarrow 2.37286$.

The Throne Beckons

Notably, the authors concede that they have not applied their refined Laser Method on larger powers of the Coppersmith Winograd tensor e.g. $CW^{\otimes 64}$... Presumably, this is a computationally intensive task, as searching for a good $\pi$ with that many terms could become intractable,,, maybe? Maybe not. Maybe they just wanted to publish their findings via the minimum necessary candidate tensor. Maybe the academic heart striving for optimality is tempered only by the need to publish.

Even cooler still is that it seems like only ... a few hundred? people are aware of this advancement – judging based on the number of citations the paper currently has. I also suspect that many of these citations might be similar to those referencing the CW paper: handwavey black box appropriation of optimal matrix multiplication performance to pad the literature survey of other adjacent publications.

Sitting next to me, generating heat but not doing too much unless a new Smite god gets dropped anytime soon, is a desktop harboring a GPU that's just begging to mine $\pi$ distributions so that my name could be listed among the ranks of Strassen, LeGall, and Vassilevska Williams.

This probably betrays some fundamental misunderstanding of the problem at a computational level, but I like the idea...

The authors do note that this method is yielding diminishing returns on the order of hundredths of thousands of arithmetic operations being shaved off of $\omega$. In 2014,⁷ Le Gall proved that 2.3725 is the lower bound using $CW^{\otimes 32}_5$ which has been the tensor d'etre for advancements made in the past decade:

• 2.37288 by Vassilevska Williams in '12
• 2.37287 by Le Gall in '14
• 2.37286 by Vassilevska Williams in '20

but just imagine:

• 2.37285 by [your name here] in '23

Additionally, we know that universal generalizations of the Laser Method can't improve past $\omega = 2.18$, so past this lofty point, we'd need a completely different tensor structure than the $CW$, and a different approach, since the Laser Method can find the optimal bound for any family of tensor for which it applies.

AlphaTensor: Those fuckers at DeepMind did it again

AlphaTensor discovers, from scratch, many provably correct Matrix Multiplication Algorithms that improve upon existing algorithms in terms of the number of scalar multiplications required.⁸

AlphaTensor is a Reinforcement Learning agent derived from the famed Go agent AlphaZero for solving planning procedure problems. AlphaTensor is trained to play a single player game whose objective is to find the tensor decomposition for arbitrary matrices in an instance of the Matrix Multiplication Problem within a finite factor space

The results are pretty neat, showing that AlphaTensor succeeded in finding better algorithm embeddings than several State of the Art techniques, including the 4x4 problem, for which it was previously thought that Strassen's algorithm from 1969 was optimal.

The task of finding Matrix Multiplication algorithm embeddings lends itself to DRL since it boils down to finding low rank decomposition of a specific three-dimensional tensor (which is our MMT). While not explicitly mentioned in Vassilevska Williams' paper, it may have been obvious that the reduction strategy for finding sub-cubic Matrix Multiplication Algorithms is $NP$-Hard. The search space is so large ( $|A| \approx 10^{12}$ for most interesting instances of the Matrix Multiplication Problem) that the optimal algorithm for even the 3x3x3 problem is unknown. All previous attempts have relied on human search, continuous optimzation, and combinatorial searches – aided by human-designed heuristics which are probably not optimal.

Algorithm 1

(for reference)