# Fast Matrix Multiplicator Evaluator

A triplet of real tensors $$\(A,B,C)\$$ of size $$\\gamma \times \gamma \times \sigma\$$ represents a matrix multiplication algorithm with $$\\sigma\$$ elementary multiplications iff the following function evaluates to zero:

$$\text{error}(A,B,C) \equiv \sum_{i,j,k,l,m,n=1}^\gamma\left(\sum_{r=1}^{\sigma}{A_{ij}^r B_{kl}^r C_{nm}^r} - \delta_{ni} \delta_{jk} \delta_{lm}\right)^2,$$ where $$\\delta_{ab} = 1\$$ if $$\a=b\$$ and 0 otherwise.

## Reference implementation

A reference implementation in C for the above formula is given below:

double error(size_t gamma, size_t sigma, double * A, double * B, double * C) {
double E = 0.0;
for (int i = 0; i < gamma; ++i) {
for (int j = 0; j < gamma; ++j) {
for (int k = 0; k < gamma; ++k) {
for (int l = 0; l < gamma; ++l) {
for (int m = 0; m < gamma; ++m) {
for (int n = 0; n < gamma; ++n) {
double s = 0.0, e = 0.0;
for (int r = 0; r < sigma; ++r) {
int rij = r*gamma*gamma + i*gamma + j;
int rkl = r*gamma*gamma + k*gamma + l;
int rnm = r*gamma*gamma + n*gamma + m;
s += A[rij] * B[rkl] * C[rnm];
}
e = s - (n==i) * (j==k) * (l==m);
E += e*e;
}}}}}}
return E;
}


## Sample input & output for $$\(\gamma=2,\sigma=7)\$$

double A[] = {1,-1,-1,0,1,-1,-1,1,0,-1,0,0,1,0,-1,0,1,-1,0,0,0,0,0,1,0,0,1,0};
double B[] = {-1,-1,-1,0,0,0,-1,0,1,1,1,1,1,0,1,0,-1,-1,0,0,0,0,0,1,0,-1,0,0};
double C[] = {1,-1,1,0,0,0,-1,0,0,-1,0,0,1,-1,0,0,-1,0,-1,0,0,0,0,1,1,-1,1,-1};
assert (0.0 == error(2, 7, A, B, C));


## Challenge

Write a faster version of the reference implementation. While measuring the speed of your code the elements of $$\(A,B,C)\$$ should be doubles sampled uniformly from the $$\[-1,1]\$$ interval.

## Competition categories

Your implementation can compete in one or more of the following categories:

1. $$\(\gamma=3,\sigma=22)\$$
2. $$\(\gamma=4,\sigma=48)\$$
3. $$\(\gamma=5,\sigma=97)\$$
4. $$\(\gamma=6,\sigma=159)\$$

## Rules

• The submission must state the $$\(\gamma,\sigma)\$$ category/categories in which your code is competing.

• The submission must contain easy-to-follow build instructions.

• Your code must display the speed of your code relative to the reference implementation (e.g. being 3.2x faster than the reference implementation).

• While measuring speed at least 10000 random evaluations must be made.

• For any given input $$\(A,B,C)\$$ the difference between the output of the reference implementation and your implementation must be smaller than $$\10^{-4}\$$.

• Any programming language and open source library (e.g. OpenBLAS) can be used.

• Reference platform: Intel i7-8559U / Linux 5.9.0, GCC 10.3.0. The reference implementation is compiled with the -O3 flag.

## Scoring

• Submissions will be ordered based on their performance on my Intel i7-8559U.
• Welcome to Code Golf! This looks like a good challenge (especially for fastest-code!). You might want to consider having the scoring rely on the timing on your computer, as others' could vary a lot depending on how good their hardware is. Nov 3 '21 at 16:21
• I'd suggest that, in order for the scoring criteria to be truly objective, you run the submissions on your computer and post timings for each answer. That way, programs' scores aren't affected by the poster's computer, only by your own and all the same Nov 3 '21 at 16:23
• @Noodle9: fixed Nov 3 '21 at 17:12
• @RedwolfProgrammed: fixed Nov 3 '21 at 17:25
• The requirement “Your code must display the speed of your code relative to the reference implementation” is hard to satisfy in languages other than C or maybe C++. Benchmarking should be up to the user, not a responsibility of the code itself. Nov 3 '21 at 19:53

# C (gcc), score ≈ 6.7×, 8.5×, 13.6x, 19.5×

For the four categories (3, 22), (4, 48), (5, 97), (6, 159).

This algorithm runs in $$\O(γ^3σ + γ^2σ^2)\$$ rather than $$\O(γ^6σ)\$$, using the rearranged formula:

$$\begin{multline*} \textrm{error}(A, B, C) = γ^3 - 2\sum_{r=1}^σ \sum_{i,j,l=1}^γ A_{ij}^rB_{jl}^rC_{il}^r \\ + \sum_{r,s=1}^σ\left(\sum_{i,j=1}^γ A_{ij}^rA_{ij}^s\right)\left(\sum_{k,l=1}^γ B_{kl}^rB_{kl}^s\right)\left(\sum_{n,m=1}^γ C_{nm}^rC_{nm}^s\right). \end{multline*}$$

double error(size_t gamma, size_t sigma, double *A, double *B, double *C) {
double E = 0.0;
for (int r = 0; r < sigma; ++r) {
for (int i = 0; i < gamma; ++i) {
for (int j = 0; j < gamma; ++j) {
double x = 0.0;
for (int l = 0; l < gamma; ++l) {
int rkl = r * gamma * gamma + j * gamma + l;
int rnm = r * gamma * gamma + i * gamma + l;
x += B[rkl] * C[rnm];
}
int rij = r * gamma * gamma + i * gamma + j;
E += A[rij] * x;
}
}
}
E = gamma * gamma * gamma - 2 * E;
for (int r = 0; r < sigma; ++r) {
for (int s = r; s < sigma; ++s) {
double xa = 0.0, xb = 0.0, xc = 0.0;
for (int ij = 0; ij < gamma * gamma; ++ij) {
int rij = r * gamma * gamma + ij;
int sij = s * gamma * gamma + ij;
xa += A[rij] * A[sij];
xb += B[rij] * B[sij];
xc += C[rij] * C[sij];
}
E += (1 + (r != s)) * xa * xb * xc;
}
}
return E;
}


Try it online!