Existence of Non-Commutative $4 \times 4$ Matrix Multiplication Algorithm

Question

This paper by a Russian gentleman gives an optimal (?) algorithm for $3$ $\times$ $3$ matrix multiplication. It beats a previously known method by reducing the total number of discrete operations from $23$ to $22$.

For technical, architectural and programming-level details these methods are faster in some cases than Strassen's algorithm. Even the naive $O(n^2)$ methods for sufficiently small $n$ matrices may still perform better (see for example this SO question).

So, naturally I am curious to know if there is literature on $4$ $\times$ $4$ versions of the algorithms cited above. Is there?

I'm reminded of this question of mine. I wasn't asking the same thing, but the problems of inverting and multiplying are strongly related. — Jacopo Notarstefano, Dec 15 '15 at 02:01
I think it's the number of multiplications that is $22$, not the number of operations. The reason that is important is that such a method of $3 \times 3$ mulltiplication can be used (as in Strassen's algorithm) as the building block to get an $O(n^p)$ method of multiplication for large matrices that is better than the naive $O(n^3)$. It's not meant to provide a fast way to multiply small matrices. — Robert Israel, Dec 15 '15 at 02:18

score 1 · Answer 1 · answered Nov 10 '17 at 21:47

Late answer, but hopefully useful:

The cited algorithm of Makarov is commutative:

P1 := (b13+a13-a23)*(a11+b31-b32+b33);
P2 := (b12+a12+a22)*(a21-b21+b22-b23);
P3 := (b12+a12+a32)*(a31-b21+b22-b23);
P4 := (b13-a23-a33)*(a31-b31+b32-b33);
P5 := (b11-a13+a23)*a11;
P6 := (b11+a12+a22)*a21;
P7 := (b11+a12+a32+a23+a33)*a31;
P8 := b12*(a11+b21-b22+b23);
P9 := b13*(a21+b31-b32+b33);
P10 := a12*b21;
P11 := a23*b31;
P12 := (a13-a23)*(a11+b31);
P13 := (a12+a22)*(b21-a21);
P14 := (b12+a12)*(b21-b22+b23);
P15 := a22*b23;
P16 := (b13-a23)*(b31-b32+b33);
P17 := a23*b32;
P18 := (a32-a23-a33)*b23;
P19 := (a13+a33-a12-a32)*b32;
P20 := (a12+a32)*(b21-a31+b23+b32);
P21 := (a23+a33)*(a31+b23-b31+b32);
P22 := (a23+a33-a12-a32)*(b23+b32);

c11 := P5+P10+P11+P12;
c12 := P8+P10-P14+P17-P18+P19-P22;
c13 := P1-P11-P12-P16+P17-P18+P19-P22;
c21 := P6-P10+P11+P13;
c22 := P2-P10+P13+P14+P15+P17;
c23 := P9-P11+P15-P16+P17;
c31 := P7-P10-P11+P20-P21+P22;
c32 := P3-P10+P14-P17+P18+P20+P22;
c33 := P4+P11+P16-P17+P18+P21;

As matrix elements of A and B occur on both sides of the products, the algorithm cannot be used recursively to multiply block matrices. This is only feasible with non-commutative algorithms. Currently, no exact non-commutative algorithms for 3x3 with less than 23 products are known.

Here is an approximate algorithm due to Schönhage:

# 3x3x3_21 according to Schönhage's paper
# "Partial and total matrix multiplication" from 1981 (cf. example 2.3 on page 5).
#
# Error constant e = 1E-12
# Maximum error for a single equation: 1E-12 = 1 x e
#    
P1 := (a11 + 1E-24a12) * (1E-24b11 + b21);
P2 := (a21 + 1E-24a22) * (1E-24b12 + b22);
P3 := (a31 + 1E-24a32) * (1E-24b13 + b23);
P4 := (a11 + 1E-24a13) * (b31);
P5 := (a21 + 1E-24a23) * (b32);
P6 := (a31 + 1E-24a33) * (b33);
P7 := (a11) * (b21 + b31);
P8 := (a21) * (b22 + b32);
P9 := (a31) * (b23 + b33);
P10 := (a11 + 1E-24a22) * ( - 1E-12b12 + b21);
P11 := (a11 + 1E-24a32) * ( - 1E-12b13 + b21);
P12 := (1E-24a12 + a21) * ( - 1E-12b11 + b22);
P13 := (a21 + 1E-24a32) * ( - 1E-12b13 + b22);
P14 := (1E-24a12 + a31) * ( - 1E-12b11 + b23);
P15 := (1E-24a22 + a31) * ( - 1E-12b12 + b23);
P16 := (a11 + 1E-24a23) * (1E-12b12 + b31);
P17 := (a11 + 1E-24a33) * (1E-12b13 + b31);
P18 := (1E-24a13 + a21) * (1E-12b11 + b32);
P19 := (a21 + 1E-24a33) * (1E-12b13 + b32);
P20 := (1E-24a13 + a31) * (1E-12b11 + b33);
P21 := (1E-24a23 + a31) * (1E-12b12 + b33);

c11 := 1E+24 * (P1 + P4 - P7);
c12 := 1E+24 * ( - 1E-12P4 - P8 + P12 + 1E-12P16 + P18);
c13 := 1E+24 * ( - 1E-12P4 - P9 + P14 + 1E-12P17 + P20);
c21 := 1E+24 * ( - 1E-12P5 - P7 + P10 + P16 + 1E-12P18);
c22 := 1E+24 * (P2 + P5 - P8);
c23 := 1E+24 * ( - 1E-12P5 - P9 + P15 + 1E-12P19 + P21);
c31 := 1E+24 * ( - 1E-12P6 - P7 + P11 + P17 + 1E-12P20);
c32 := 1E+24 * ( - 1E-12P6 - P8 + P13 + P19 + 1E-12P21);
c33 := 1E+24 * (P3 + P6 - P9);

The following non-commutative algorithm for 4x4 with 49 products is based on recursive application of Strassen's 2x2 algorithm with 7 products:

P01 := (a11+a22+a33+a44) * (b11+b22+b33+b44);
P02 := (a21+a22+a43+a44) * (b11+b33);
P03 := (a11+a33) * (b12-b22+b34-b44);
P04 := (a22+a44) * (-b11+b21-b33+b43);
P05 := (a11+a12+a33+a34) * (b22+b44);
P06 := (-a11+a21-a33+a43) * (b11+b12+b33+b34);
P07 := (a12-a22+a34-a44) * (b21+b22+b43+b44);
P08 := (a31+a33+a42+a44) * (b11+b22);
P09 := (a41+a42+a43+a44) * (b11);
P10 := (a31+a33) * (b12-b22);
P11 := (a42+a44) * (-b11+b21);
P12 := (a31+a32+a33+a34) * (b22);
P13 := (-a31-a33+a41+a43) * (b11+b12);
P14 := (a32+a34-a42-a44) * (b21+b22);
P15 := (a11+a22) * (b13+b24-b33-b44);
P16 := (a21+a22) * (b13-b33);
P17 := (a11) * (b14-b24-b34+b44);
P18 := (a22) * (-b13+b23+b33-b43);
P19 := (a11+a12) * (b24-b44);
P20 := (-a11+a21) * (b13+b14-b33-b34);
P21 := (a12-a22) * (b23+b24-b43-b44);
P22 := (a33+a44) * (-b11-b22+b31+b42);
P23 := (a43+a44) * (-b11+b31);
P24 := (a33) * (-b12+b22+b32-b42);
P25 := (a44) * (b11-b21-b31+b41);
P26 := (a33+a34) * (-b22+b42);
P27 := (-a33+a43) * (-b11-b12+b31+b32);
P28 := (a34-a44) * (-b21-b22+b41+b42);
P29 := (a11+a13+a22+a24) * (b33+b44);
P30 := (a21+a22+a23+a24) * (b33);
P31 := (a11+a13) * (b34-b44);
P32 := (a22+a24) * (-b33+b43);
P33 := (a11+a12+a13+a14) * (b44);
P34 := (-a11-a13+a21+a23) * (b33+b34);
P35 := (a12+a14-a22-a24) * (b43+b44);
P36 := (-a11-a22+a31+a42) * (b11+b13+b22+b24);
P37 := (-a21-a22+a41+a42) * (b11+b13);
P38 := (-a11+a31) * (b12+b14-b22-b24);
P39 := (-a22+a42) * (-b11-b13+b21+b23);
P40 := (-a11-a12+a31+a32) * (b22+b24);
P41 := (a11-a21-a31+a41) * (b11+b12+b13+b14);
P42 := (-a12+a22+a32-a42) * (b21+b22+b23+b24);
P43 := (a13+a24-a33-a44) * (b31+b33+b42+b44);
P44 := (a23+a24-a43-a44) * (b31+b33);
P45 := (a13-a33) * (b32+b34-b42-b44);
P46 := (a24-a44) * (-b31-b33+b41+b43);
P47 := (a13+a14-a33-a34) * (b42+b44);
P48 := (-a13+a23+a33-a43) * (b31+b32+b33+b34);
P49 := (a14-a24-a34+a44) * (b41+b42+b43+b44);

c11 := P01+P04-P05+P07+P22+P25-P26+P28-P29-P32+P33-P35+P43+P46-P47+P49;
c12 := P03+P05+P24+P26-P31-P33+P45+P47;
c13 := P15+P18-P19+P21+P29+P32-P33+P35;
c14 := P17+P19+P31+P33;
c21 := P02+P04+P23+P25-P30-P32+P44+P46;
c22 := P01-P02+P03+P06+P22-P23+P24+P27-P29+P30-P31-P34+P43-P44+P45+P48;
c23 := P16+P18+P30+P32;
c24 := P15-P16+P17+P20+P29-P30+P31+P34;
c31 := P08+P11-P12+P14+P22+P25-P26+P28;
c32 := P10+P12+P24+P26;
c33 := P01+P04-P05+P07-P08-P11+P12-P14+P15+P18-P19+P21+P36+P39-P40+P42;
c34 := P03+P05-P10-P12+P17+P19+P38+P40;
c41 := P09+P11+P23+P25;
c42 := P08-P09+P10+P13+P22-P23+P24+P27;
c43 := P02+P04-P09-P11+P16+P18+P37+P39;
c44 := P01-P02+P03+P06-P08+P09-P10-P13+P15-P16+P17+P20+P36-P37+P38+P41;

Existence of Non-Commutative $4 \times 4$ Matrix Multiplication Algorithm

1 Answers1