# New 4x4 algorithm found
**URL:** https://discourse.julialang.org/t/new-4x4-algorithm-found/129012
**Category:** Offtopic
**Created:** [May 14, 2025, 7:53pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012 "2025-05-14T19:53:44Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@ShalokShalom](https://discourse.julialang.org/u/ShalokShalom)
#### Post date: [May 14, 2025, 7:53pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/1 "2025-05-14T19:53:44Z")
Someone found a new algorithm for 4x4 matrix multiplication. .
> **[AlphaEvolve: A Gemini-powered coding agent for designing advanced algorithms](https://deepmind.google/discover/blog/alphaevolve-a-gemini-powered-coding-agent-for-designing-advanced-algorithms/)**
>
> New AI agent evolves algorithms for math and practical applications in computing by combining the creativity of large language models with automated evaluators
[ ](https://www.youtube.com/watch?v=sGCmu7YKgPA)
This will save the industry billions.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 14, 2025, 8:08pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/2 "2025-05-14T20:08:21Z")
> [@ShalokShalom](#):
>
> This will save the industry billions.
Probably an overstatement. See the discussion of the previous AlphaTensor results: [Matrix multiply breakthrough, AlphaTensor (could also do for other algorithms): "AlphaTensor discovered algorithms that are more efficient than the state of the art for many matrix sizes."](https://discourse.julialang.org/t/matrix-multiply-breakthrough-alphatensor-could-also-do-for-other-algorithms-alphatensor-discovered-algorithms-that-are-more-efficient-than-the-state-of-the-art-for-many-matrix-sizes/88455)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 14, 2025, 9:02pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/3 "2025-05-14T21:02:58Z")
Also… are they numerically stable? You not only need to care about performance but also the error growth of the new algorithms, and for Strassen’s it requires modifications to make that work:
> **[HiPC-2018.pdf](https://www.comp.nus.edu.sg/~wongwf/papers/HiPC-2018.pdf)**
>
> 985.12 KB
AlphaEvolve didn’t have any measure of numerical stability, and lots of these “faster matrix multiplication” algorithms have difficulty in this area. So I would presume AlphaEvolve will have poor error growth behavior without extra changes and corrections, and those of course slow it down.
---
### Author:  [@e3c6](https://discourse.julialang.org/u/e3c6)
#### Post date: [May 16, 2025, 11:18pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/4 "2025-05-16T23:18:09Z")
Related:
> **[$XX^{t}$ Can Be Faster](https://arxiv.org/abs/2505.09814)**
>
> We present a new algorithm RXTX that computes product of matrix by its transpose $XX^{t}$. RXTX uses $5\\%$ less multiplications and additions than State-of-the-Art and achieves accelerations even for small sizes of matrix $X$. The algorithm was...
---
### Author:  [@photor](https://discourse.julialang.org/u/photor)
#### Post date: [May 17, 2025, 2:03am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/5 "2025-05-17T02:03:13Z")
Where can I find AlphaEvolve’s original paper on this new algorithm?
---
### Author:  [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)
#### Post date: [May 17, 2025, 4:44am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/6 "2025-05-17T04:44:49Z")
> [@ChrisRackauckas](#):
>
> Also… are they numerically stable? You not only need to care about performance but also the error growth of the new algorithms
But wouldn’t it at least work for something like large language models/etc? In large language models, the precision doesn’t really matter that much. People have even been able to quantize the model to 4 bits with minimal performance loss.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [May 18, 2025, 8:12am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/7 "2025-05-18T08:12:38Z")
> [@photor](#):
>
> Where can I find AlphaEvolve’s original paper on this new algorithm?
> **[AlphaEvolve.pdf](https://storage.googleapis.com/deepmind-media/DeepMind.com/Blog/alphaevolve-a-gemini-powered-coding-agent-for-designing-advanced-algorithms/AlphaEvolve.pdf)**
>
> 3.19 MB
---
### Author:  [@greatpet](https://discourse.julialang.org/u/greatpet)
#### Post date: [May 18, 2025, 8:30pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/8 "2025-05-18T20:30:32Z")
I’m very surprised by the use of a general-purpose LLM, Gemini, for dealing with a narrowly defined mathematical problem of matrix multiplication. Gemini’s training texts presumably include e.g. Shakespear’s works and New York Times editorials, both of which could have subtly affected the discovery of a matrix multiplication algorithm, as the transformer neural networks could potentially couple everything.
---
### Author:  [@mrufsvold](https://discourse.julialang.org/u/mrufsvold)
#### Post date: [May 18, 2025, 11:38pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/9 "2025-05-18T23:38:21Z")
When all you have is a hammer that you spent hundreds of millions of dollars on, everything looks like a golden nail that might have some ROI.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2025, 3:52am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/10 "2025-05-19T03:52:26Z")
If anybody is interested in trying to do “algorithm discovery” in Julia, here’s one way you can do it via evolution. No LLMs, just regular genetic algorithms.
```julia
using SymbolicRegression
using MLJBase
```
Let’s try to discover 2x2 matrix algorithms, for some constraints on the operations.
First, we make some random 2x2 matrices and compute their products. These are the labels we will try to regress
```julia
N = 1_000
X = randn(Float32, N, 2 * (2 * 2))
# Label = 4-tuple representing matrix product
y = [
let
A = reshape(x[1:4], 2, 2)
B = reshape(x[5:8], 2, 2)
C = A * B
tuple(C[:]...)
end
for x in eachrow(X)
]
```
For simplicity, let’s say that we want to try to rediscover the Strassen formula for 2x2 multiplication (just scalars). We can use a “template expression” to try to do this, which lets us evolve multiple expressions simultaneously. Here, we will evolve 14 expressions: 10 intermediates (some redundancy for easier evolution), and 4 outputs, one per element.
```julia
expression_spec = @template_spec(expressions=(
m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, # intermediate expressions
c11, c12, c21, c22 # output expressions
)) do a11, a12, a21, a22, b11, b12, b21, b22 # input features
features = (a11, a12, a21, a22, b11, b12, b21, b22)
# First, we compute the intermediates
M1 = m1(features...)
M2 = m2(features...)
M3 = m3(features...)
M4 = m4(features...)
M5 = m5(features...)
M6 = m6(features...)
M7 = m7(features...)
M8 = m8(features...)
M9 = m9(features...)
M10 = m10(features...)
intermediates = (M1, M2, M3, M4, M5, M6, M7, M8, M9, M10)
# Then, compute the outputs
C11 = c11(intermediates...)
C12 = c12(intermediates...)
C21 = c21(intermediates...)
C22 = c22(intermediates...)
# Output needs to be a `ValidVector` type; these automate vectorization and propagate validity checks
return ValidVector(
map(i -> (C11.x[i], C12.x[i], C21.x[i], C22.x[i]), eachindex(C11.x)),
C11.valid && C12.valid && C21.valid && C22.valid
)
end
```
We also need to define a custom loss function to support our 4-tuple outputs:
```julia
tuple_loss(prediction, target) = sum(abs2, map(-, prediction, target))
```
Finally, we throw it into a regressor object. So that we can try to discover Strassen, we can make the `*` operation more “complex” than `+` and `-`:
```julia
model = SRRegressor(
binary_operators=(+, -, *),
unary_operators=(),
complexity_of_operators=[(+) => 1, (-) => 1, (*) => 2], # Make * more expensive
complexity_of_constants=10, # Prevents constants from being used
should_optimize_constants=false, # Also, no need for BFGS constant tuning
expression_spec=expression_spec,
elementwise_loss=tuple_loss,
loss_type=Float32,
batching=true,
batch_size=32,
maxsize=200, # Total complexity overall all expressions
niterations=100000, # How long the search lasts
)
mach = machine(model, X, y; scitype_check_level=0)
fit!(mach)
```
With this the search was able to ‘rediscover’ regular matrix multiplication in about 10minutes on my laptop. With enough compute to burn, this might get to Strassen. I think doing it will need a lot of compute though. Would help to also try tuning with larger values for `population_size` (and maybe `populations` and `ncycles_per_iteration`) to support the large search space needed.
---
### Author:  [@greatpet](https://discourse.julialang.org/u/greatpet)
#### Post date: [May 19, 2025, 10:01am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/11 "2025-05-19T10:01:00Z")
> [@MilesCranmer](#):
>
> With this I was able to ‘rediscover’ regular matrix multiplication in about 10 minutes.
How many “test programs” were generated and executed in the 10 minutes?
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 19, 2025, 12:01pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/12 "2025-05-19T12:01:43Z")
The search above looks to be about ~2e7 expressions. The vast majority of which are only evaluated on a 32-example mini-batch per the `batch_size` argument.
The equivalent brute force: the expression is, if I remember correctly, about 50 “complexity”. So, if you use Catalan numbers and make a bunch of simplifying assumptions, you should find something like 1e50 (give or take 10 orders of magnitude) expressions are possible when there are 8 leaf types (a11 through b22) and 3 node types (+, -, \*). Traditional evolutionary algorithms can do a fairly good job in reducing the search space despite not having any “intelligence.”
---
### Author:  [@Palli](https://discourse.julialang.org/u/Palli)
#### Post date: [May 19, 2025, 12:04pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/13 "2025-05-19T12:04:04Z")
> Despite being general-purpose, AlphaEvolve goes beyond [25], improving the SOTA for 14 matrix multiplication algorithms; notably, for 4 × 4 matrices, AlphaEvolve improves Strassen (1969)’s algorithm by discovering an algorithm using 48 multiplications to multiply 4 × 4 complex-valued matrices.
I.e. improvement of 48/49 = 2.04% improvement/reduction (i.e. of one multiply).
But for the 3×4 times 4×6 they have a higher 4.5% speedup/improvement (three multiplies skipped). And for 4×5 times 5×6 there’s 3.2% improvement, again 3 fewer multiplies (all the other new algorithms have only 2 or 1 fewer), see table 2 for other improvements.
> Table 2 … For ⟨3, 4, 7⟩, ⟨4, 4, 4⟩, and ⟨4, 4, 8⟩, the algorithms discovered by AlphaEvolve use complex-valued multiplications which can be used for exact multiplication of complex or real-valued matrices.
Should we change the title to e.g. [New matrix algorithms found including a 4x4 (complex-valued) algorithm](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012)?
If you know your matrices are real-valued, do we already have better still, despite we can use the one with 48 or older 49? See their footnote 3
> There exist algorithms using fewer than 49 multiplications, but they do not correspond to decompositions of the matrix multiplication tensor, and they cannot be applied recursively to multiplying larger matrices.
I.e. their new 48 isn’t always lowest for 4 × 4 matrix multiply, only when real- or complex valued. See their older work for Tensors from 2022 with AlphaTensor (fig 3, there in regular (open access) Nature):
[https://www.nature.com/articles/s41586-022-05172-4](https://www.nature.com/articles/s41586-022-05172-4)
> multiplying 4 × 4 matrices using 47 multiplications in Z2, thereby outperforming Strassen’s two-level algorithm[2](https://www.nature.com/articles/s41586-022-05172-4#ref-CR2), which involves 72 = 49 multiplications.
They also have an update on AlphaTensor today (intriguingly in [_Nature Machine Intelligence_](https://www.nature.com/natmachintell) volume **7** , not sure why) for Quantum circuit optimization:
[https://www.nature.com/articles/s42256-025-01001-1](https://www.nature.com/articles/s42256-025-01001-1)
---
### Author:  [@photor](https://discourse.julialang.org/u/photor)
#### Post date: [May 19, 2025, 12:24pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/14 "2025-05-19T12:24:27Z")
If I remember it correctly, the 47 scalar multiplications are for binaries, not floating numbers, so it’s much less important than the one found this time.
---
### Author:  [@apo383](https://discourse.julialang.org/u/apo383)
#### Post date: [May 19, 2025, 8:23pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/15 "2025-05-19T20:23:13Z")
> [@greatpet](#):
>
> I’m very surprised by the use of a general-purpose LLM
The LLM is used for code generation, and doesn’t need to know much about matrix multiplication. IMO the key is to couple the LLM to a feedback loop. They apply automated evaluation metrics (rewards), which are usually tricky to specify. This is applied to four problem domains suitable for automated iteration, one of which is matrix multiply. They also keep a database of candidate programs, and have a prompt sampling approach that can feed these in to the LLM. They are optimizing these prompts, so it’s an extra optimization loop on top of the program optimization.
As for the comment elsewhere above about numerical stability, my take is that it could be built into the evaluation, so you could optimize a weighted combo of speed and stability. It is also true that their particular application is AI, so they care more about speed than accuracy.
I wonder if it would make sense to optimize on energy usage. Perhaps there are common bit patterns in AI, where the matmul could reduce bit flipping, or favor pure adds over multiply-adds or such. I presume they want both speed and energy economy.
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [May 19, 2025, 11:13pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/16 "2025-05-19T23:13:19Z")
> [@MilesCranmer](#):
>
> First, we make some random 2x2 matrices and compute their products. These are the labels we will try to regress
The algorithm discovery papers typically cast the problem of optimizing matrix multiplication as that of finding the lowest-rank decomposition of a particular tensor. It works as follows: You can express matrix multiplication as contraction of the matrices with a tensor L:
C\_{ij} = L\_{ijklmn} A\_{kl} B\_{mn}
(Einstein summation implied), where L\_{ijklmn} = \delta\_{ik} \delta\_{lm} \delta\_{jn} . Now your challenge is to find a decomposition as follows:
L\_{ijklmn} = \sum\_{p = 1}^T \gamma\_{ijp} \alpha\_{klp} \beta\_{mnp} \, ,
where T is called the _rank_ of the decomposition. Given such a decomposition, you can reorganize the matrix product such that p is the last index you contract over:
C\_{ij} = \sum\_{i = 1}^T \gamma\_{ijp} \left[\left(\sum\_{kl} \alpha\_{klp} A\_{kl} \right)\left(\sum\_{mn} \beta\_{mnp} B\_{mn}\right)\right] \, .
Clearly, the total number of multiplications is proportional to T. Further, if every element in the \gamma, \alpha, \beta tensors is \in \{-1, 0, +1\}, the tensor contractions simplify to just indexing, addition, and subtraction, so the only multiplication per p is the one between (\alpha\_{klp} A\_{kl}) and (\beta\_{mnp} B\_{mn}). Thus, the total number of multiplications is T.
Moreover, A and B are typically large block matrices or some other “heavy” type such that multiplying elements of A and B is much more expensive than any other operation. Under this assumption, only the multiplications between (\alpha\_{klp} A\_{kl}) and (\beta\_{mnp} B\_{mn}) count in terms of complexity. In the research literature these are called _active_ multiplications, as opposed to all other operations, including any nontrivial scalar multiplications in the tensor contractions, which are _inactive_. Since the number of active multiplications is always equal to the rank T, this remains the minimization objective even when you allow for rational/real/complex tensor elements. For example, in the much-discussed T = 48-decomposition for 4x4 matrices, the tensor elements are \in \{0, \pm 1/2, \pm i/2, \pm (1 \pm i) / 2\}, so the algorithm contains lots of inactive scalar multiplications in addition to the 48 active multiplications.
(See @stevengj’s comment in the previous thread for more context around the application domain of these algorithms: [https://discourse.julialang.org/t/matrix-multiply-breakthrough-alphatensor-could-also-do-for-other-algorithms-alphatensor-discovered-algorithms-that-are-more-efficient-than-the-state-of-the-art-for-many-matrix-sizes/88455/2.](https://discourse.julialang.org/t/matrix-multiply-breakthrough-alphatensor-could-also-do-for-other-algorithms-alphatensor-discovered-algorithms-that-are-more-efficient-than-the-state-of-the-art-for-many-matrix-sizes/88455/2.))
Finally, in this formulation, it’s straightforward to check the correctness of the algorithm: you just verify that
\sum\_{p = 1}^T \gamma\_{ijp} \alpha\_{klp} \beta\_{mnp} = \delta\_{ik} \delta\_{lm} \delta\_{jn} \, .
for each of the 64 combinations of i,j,k,l,m,n \in \{1, 2\}. If you restrict your search to rational tensor elements and perform this verification in exact rational arithmetic, you can be certain that the algorithm you’ve found is exactly correct (in infinite-precision arithmetic, that is—whether it’s stable when taking into account floating-point roundoff is another question).
@MilesCranmer can SymbolicRegression.jl be applied to a problem like this? The tensor decomposition expression has a fixed form, so we’re not allowed to arbitrarily modify the expression tree; the only degree of freedom is to increase or decrease the number of terms T, and the objective is for it to be as small as possible. As for the coefficients, when rediscovering Strassen’s algorithm, we can restrict them to \{-1, 0, 1\}, in which case whether or not we’re on target is a binary property. However, if optimizing over continuous coefficients is preferred, we can relax the coefficients to, say [-1, 1], and try to minimize the distance from the target in a least squares sense. For 2x2 matrices, this problem has 64 equations in 12T unknowns.
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [May 20, 2025, 12:28am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/17 "2025-05-20T00:28:24Z")
> [@danielwe](#):
>
> For example, in the much-discussed T = 48-decomposition for 4x4 matrices, the tensor elements are \in \{0, \pm 1/2, \pm i/2, \pm (1 \pm i) / 2\}, so the algorithm contains lots of inactive scalar multiplications in addition to the 48 active multiplications.
Btw. this means that the following claim from the abstract, now repeated all over the internet, is at best misleading:
> AlphaEvolve developed a search algorithm that found a procedure to multiply two 4 × 4 complex-valued matrices using 48 scalar multiplications
It seems like the AlphaEvolve paper uses “scalar multiplications” the same way other papers use “active multiplications”, making it sound like the algorithm can multiply two 4x4 matrices of plain complex numbers (i.e., not block matrices) using only 48 complex multiplications.
However, here’s what an implementation looks like: [AlphaEvolve-MatrixMul-Verification/matrix\_multiplication\_algorithms.py at da09d7c476353fe888331922f16d3449abb2b21d · PhialsBasement/AlphaEvolve-MatrixMul-Verification · GitHub](https://github.com/PhialsBasement/AlphaEvolve-MatrixMul-Verification/blob/da09d7c476353fe888331922f16d3449abb2b21d/matrix_multiplication_algorithms.py#L417-L719). You can see the 48 active multiplications starting from line 525 under the comment `# Perform the 48 multiplications efficiently`. But the tensor contractions, implemented in the 100 lines before and 140 lines after, contain countless multiplications, which are just as expensive as the active multiplications if your input matrix has scalar elements.
That said, the implementation could be better optimized than this. You could double each tensor component, eliminating inactive multiplications from the tensor contractions and obtaining 8C. Finally, multiply each element by 1/8 to get C; that’s 16 inactive multiplications. In total that’s 64 complex multiplications—exactly the same as in the naive algorithm (but with a lot more additions and subtractions!).
I should add that the claims remain true and valid for block matrices where the distinction between active and inactive multiplications is meaningful. My gripe is only with the choice of words, which makes it sound like the algorithm can do things it cannot.
---
### Author:  [@MilesCranmer](https://discourse.julialang.org/u/MilesCranmer)
#### Post date: [May 20, 2025, 7:02am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/18 "2025-05-20T07:02:12Z")
> [@danielwe](#):
>
> @MilesCranmer can [SymbolicRegression.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/SymbolicRegression) be applied to a problem like this?
Potentially yes; check out [Custom Types · SymbolicRegression.jl](https://ai.damtp.cam.ac.uk/symbolicregression/dev/examples/custom_types). This demonstrates how to evolve expressions on `::String`.
So you could create a new type that implements the same interface for tensors. Then pass operators which can handle that type.
* * *
Although, for searching on non-scalar but numeric types, where you want to still use BFGS for constant opimtization during the evolution, you likely want to set:
```julia
can_optimize(::Type{<:MyType}, _) = true
```
You also need to declare things like `count_scalar_constants`, `pack_scalar_constants`, and `unpack_scalar_constants`, so that the library knows how to get a single vector of all your constants.
To implement this, you have to conform to `DynamicExpressions.ValueInterface` (which is declared and tested using Interfaces.jl). You can see this file for an example: [DynamicExpressions.jl/test/test\_non\_number\_eval\_tree\_array.jl at 605e1116041a0d36a51b0ee90617a4ca9bea39fd · SymbolicML/DynamicExpressions.jl · GitHub](https://github.com/SymbolicML/DynamicExpressions.jl/blob/605e1116041a0d36a51b0ee90617a4ca9bea39fd/test/test_non_number_eval_tree_array.jl#L26-L128).
---
### Author:  [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)
#### Post date: [May 20, 2025, 12:36pm UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/19 "2025-05-20T12:36:36Z")
Maybe we could start a new thread if we’re interested in trying to discover a new matrix multiplication algorithm with Julia, though I doubt it would go anywhere beyond the current SOTA, but who knows?
---
### Author:  [@photor](https://discourse.julialang.org/u/photor)
#### Post date: [May 21, 2025, 3:36am UTC](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012/20 "2025-05-21T03:36:29Z")
Maybe more appropriate names are “cheap multiplications” (N^2 operations each) and “expensive multiplications” (N^3 operations each with the naive method).
[Next page](https://discourse.julialang.org/t/new-4x4-algorithm-found/129012.md?page=2)