Local Search for Fast Matrix Multiplication
Marijn J. H. Heule Manuel Kauers Martina Seidl


 

Department of Computer Science, The University of Texas at Austin, United StatesInstitute for Algebra, J. Kepler University Linz, AustriaInstitute for Formal Models and Verification, J. Kepler University Linz, Austria

 

 Local Search for Fast Matrix Multiplication
    Marijn J.H. Heule1Supported by NSF grant CCF-1813993 and AFRL Award FA8750-15-2-0096. Manuel Kauers2Supported by the Austrian FWF grants P31571-N32 and F5004. Martina Seidl3Supported by the Austrian FWF grant NFN S11408-N23 
and the LIT AI Lab funded by the State of Upper Austria.
    Received: date / Accepted: date
=============================================================================================================================================================================================================================================================================================



  
Laderman discovered a scheme for computing the product of two 3×3 matrices
using only 23 multiplications in 1976. Since then, some
more such schemes were proposed, but nobody knows how many such schemes
there are and whether there exist schemes with fewer than 23
multiplications. In this paper we present two independent
SAT-based methods for finding new schemes using 23 multiplications. 
Both methods allow
computing a few hundred new schemes individually, and many
thousands when combined. Local search SAT solvers outperform CDCL
solvers consistently in this application.
 



§ INTRODUCTION


Matrix multiplication is a fundamental operation with applications in nearly any
area of science and engineering. However, after more than 50 years of work on
matrix multiplication techniques (see,
e.g., <cit.>),
the complexity of matrix multiplication is still a mystery. Even for small matrices, the problem
is not completely understood, and understanding these cases better can provide valuable hints
towards more efficient algorithms for large matrices.

The naive way for computing the product C of two 2× 2 matrices
A,B requires 8 multiplications:

    [ a_11 a_12; a_21 a_22 ][ b_11 b_12; b_21 b_22 ] =
     [      a_11b_11 + a_12b_21 .5em a_11b_12 + a_21b_22;      a_21b_11 + a_22b_21 .5em a_21b_12 + a_22b_22 ]
     =
     [ c_11 c_12; c_21 c_22 ]

Strassen observed 50 years ago that C can also be computed with only 7
multiplications <cit.>.
His scheme proceeds in two steps.
In the first step he introduces auxiliary variables M_1, …, M_7
which are defined as the product of certain linear combinations of
the entries of A and B. In the second step the entries of C are obtained as
certain linear combinations of the M_i:

 M_1    = (a_11 + a_22) (b_11 + b_22)
             c_11     =  M_1 + M_4 - M_5 + M_7 

 M_2    = (a_21 + a_22) (b_11)
      c_12     =  M_3 + M_5 

 M_3    = (a_11) (b_12 - b_22)
      c_21     =   M_2 + M_4 

 M_4    = (a_22) (b_21 - b_11)
      c_22     =   M_1 - M_2 + M_3 + M_6 

 M_5    = (a_11 + a_12)(b_22)

 M_6    = (a_21 - a_11) (b_11+ b_12)

 M_7    = (a_12 - a_22) (b_21 + b_22)

Recursive application of this scheme gave rise to the first algorithm for multiplying
arbitrary n× n matrices in subcubic complexity.
Winograd <cit.> showed
that Strassen's scheme is optimal in the sense
that there does not exist a similar scheme with fewer than 7 multiplications,
and de Groote <cit.>
showed that Strassen's scheme is essentially unique.

Less is known about 3× 3 matrices.
The naive scheme requires 27 multiplications, and in 1976
Laderman <cit.> found one with 23.
Similar as Strassen, he defines M_1,…,M_23 as
products of certain linear combinations of
the entries of A and B. The entries of C=AB are then obtained
as linear combinations of M_1,…,M_23.
It is not known whether 23 is optimal (the best lower bound is 19 <cit.>).
It is known however that Laderman's scheme is not unique.
A small number of intrinsically different schemes have been found
over the years. Of particular interest are schemes in which all
coefficients in the linear combinations are +1, -1, or 0.
The only four such schemes (up to equivalence) we are aware of
are due to Laderman, Smirnov <cit.>, Oh et al. <cit.>, and Courtois
et al. <cit.>. 

While Smirnov and Oh et al. found their multiplicatoin schemes 
with computer-based search using non-linear numerical optimization 
methods, Courtois found his multiplication scheme using a SAT solver. This is
also what we do here. We present two approaches which allowed us to generate more
than 13,000 mutually inequivalent new matrix multiplication schemes for 3×
3 matrices, using altogether about 35 years of CPU years. We believe that the
new schemes are of interest to the matrix multiplication community. We therefore
make them publicly available in various formats and grouped by invariants at


<http://www.algebra.uni-linz.ac.at/research/matrix-multiplication/>.









§ ENCODING AND WORKFLOW


^ 

To search for multiplication schemes of 3 × 3 matrices having
the above form, we define the M_i as product of linear combination
of all entries of A and B with undetermined coefficients α_ij^(ℓ), β_ij^(ℓ):
 
   M_1   =
   (α_11^(1)a_11+ ⋯+ α_33^(1)a_33)
   (β_11^(1)b_11+ ⋯+ β_33^(1)b_33)

      ⋮

   M_23   =
   (α_11^(23)a_11+ ⋯+ α_33^(23)a_33)
   (β_11^(23)b_11+ ⋯+ β_33^(23)b_33)
 
 Similarly, we define the c_ij as linear combinations of the
 M_i with undetermined coefficients γ_i,j^(ℓ):
 
   c_11   = γ_11^(1)M_1+ ⋯+ γ_11^(23)M_23,
        …,     
   c_33   = γ_33^(1)M_1+ ⋯+ γ_33^(23)M_23
 Comparing the coefficients of all terms a_i_1i_2b_j_1j_2c_k_1k_2 in the equations c_ij=∑_k a_ikb_kj
 leads to the polynomial equations
 
    ∑_ℓ=1^23α_i_1i_2^(ℓ)β_j_1j_2^(ℓ)γ_k_1k_2^(ℓ) =
     δ_i_2j_1δ_i_1k_1δ_j_2k_2

 for i_1,i_2,j_1,j_2,k_1,k_2∈{1,2,3}.  These 729 cubic equations with 621
 variables are also known as Brent equations <cit.>.  The
 δ_uv on the right are Kronecker-deltas, i.e., δ_uv=1 if u=v
 and δ_uv=0 otherwise.  Each solution of the system of these equations
 corresponds to a matrix multiplication scheme.  The equations become slightly
 more symmetric if we flip the indices of the γ_ij, and since this is
 the variant mostly used in the literature, we will also adopt it from now on.

 Another view on the Brent equations is as follows. View the
 α_i_1i_2^(ℓ), β_j_1j_2^(ℓ),
 γ_k_1k_2^(ℓ) as variables, as before, and regard
 a_i_1i_2,b_j_1j_2,c_k_1k_2 as polynomial indeterminants. Then the task
 consists of instantiating the variables in such a way that
 
    ∑_ℓ=1^23(α_11^(ℓ)a_11+⋯)(β_11^(ℓ)b_11+⋯)(γ_11^(ℓ)c_11+⋯)
     =∑_i=1^3∑_j=1^3∑_k=1^3 a_ijb_jkc_ki

 holds as equation of polynomials in the variables
 a_i_1i_2,b_j_1j_2,c_k_1k_2. Expanding the left hand side and equating
 coefficients leads to the Brent equations as stated before (but with indices of
 γ flipped, as agreed). In other words, expanding the left hand side,
 all terms have to cancel except for the terms on the right. We found it
 convenient to say that a term a_i_1i_2b_j_1j_2c_k_1k_2 has “type m”
 if m = δ_i_2j_1+δ_j_2k_1+δ_k_2i_1.  With this
 terminology, all terms of types 0, 1, 2 have to cancel each other, and all
 terms of type 3 have to survive. Note that since all 27 type 3 terms
 must be produced by the 23 summands on the left, some summands must produce
 more than one type 3 term.

 For solving the Brent equations with a SAT solver, we use Z_2 as
 coefficient domain, so that multiplication translates into `and' and
 addition translates into `xor'. When, for example, the variable α_i_1i_2^(ℓ) is true
 in a solution of the corresponding SAT instance, this indicates that
 the term a_i_1i_2 occurs in M_ℓ, and likewise for the b-variables.
 If γ_k_1k_2^(ℓ) is true, this means that M_ℓ appears in the
 linear combination for c_k_1k_2. We call  α_i_1i_2^(ℓ), β_j_1j_2^(ℓ), and γ_k_1k_2^(ℓ)
 the base variables.

 In order to bring the Brent equations into CNF, we use Tseitin transformation, i.e., we introduce definitions for
 subformulas to avoid exponential blow-up. To keep the number of fresh variables
 low, we do not introduce one new variable for every cube but only for pairs of literals,
 i.e., we encode a cube (αβγ) as u ↔ (αβ) and
 v ↔ (u γ). In this way, we can reuse u.
 Furthermore, a sum v_1 ⊕…⊕ v_m with m≥4 is encoded by
 w↔(v_1⊕ v_2⊕ v_3) and v_4 ⊕…⊕ v_m⊕ w,
 with the latter sum being encoded recursively. This encoding seems to require the
 smallest sum of the number of variables and the number of clauses—a commonly used optimality heuristic. The used scripts are available at 


<https://github.com/marijnheule/matrix-challenges/tree/master/src>.


 The generation of new schemes proceeds in several steps, with SAT solving being
 the first and main step. If the SAT solver finds a solution, we next check
 whether it is equivalent to any known or previously found solution modulo de Groote's symmetry group <cit.>.
 If so, we discard it. Otherwise, we next try to simplify the new scheme by searching for an element in
 its orbit which has a smaller number of terms. The scheme can then be used to initiate
 a new search. In the fourth step, we use
 Gröbner bases to lift the scheme from the coefficient domain Z_2 to Z.
 Finally, we cluster large sets of similar solutions into parameterized families.
 

   [xscale=1.7]
     (b) at (1,0) ∙ (1,0) node[left,xshift=-2mm] known-2ptschemes;
     (c) at (2,0) ∙;
     (d) at (3,0) ∙;
     (e) at (4,0) ∙;
     (f) at (5,0) ∙;
     (g) at (6,0) ∙ (6,0) node[right,xshift=2mm] new-2ptschemes;
     [->] (b)–(c) node[midway,above] solve;
     [->] (c)–(d) node[midway,above] filter;
     [->] (d)–(e) node[midway,above] simplify;
     [->] (e)–(f) node[midway,above] lift;
     [->] (f)–(g) node[midway,above] cluster;
     [->] (e)–++(0,-.5)–++(-3,0)–(b);
   
 

 In the present paper, we give a detailed description of the first step in this workflow.
 The subsequent steps use algebraic techniques unrelated to SAT and will be described in <cit.>.

 

§ SOLVING METHODS



 The core of a scheme is the pairing of the type 3 terms.
 
 Our first method focuses on finding schemes with new cores, while our second method searches for schemes that are similar to an
 existing one and generally has the same core.
 For all experiments we used the local search SAT solver yalsat <cit.> as this solver performed best on instances from this application.
 We also tried solving these instances using CDCL solvers, but the performance was disappointing. We observed a possible explanation:
 The runtime of CDCL solvers tends to be exponential in the average backtrack level (ABL) on unsatisfiable instances.
 For most formulas arising from other applications, ABL is small (<50), while on the matrix multiplication instances ABL is large (> 100).

 

 §.§ Random Pairings of Type 3 Terms and Streamlining


 Two of the known schemes, those of Smirnov <cit.> and Courtois et al. <cit.>,
 have the property that each type 3 term occurs exactly once and at most two type 3 terms occur in the same summand. We decided to use this pattern
  to search for new schemes: randomly pair four type 3 terms and assign the remaining type 3 terms to the other 19 summands.
  Only in very rare cases, random pairing could be extended to a valid scheme in reasonable time, say a few minutes.
  In the other cases it is not known whether the pairing cannot be extended to a valid scheme or whether finding such a scheme is very hard.

  Since the number of random pairings that could be extended to a valid scheme was very low, we
  tried 
  adding streamlining constraints <cit.> to formulas. A streamlining constraint is a set of clauses that
  guides the solver to a solution, but these clauses may not (and generally are not) implied by the formula. Streamlining constraints are
  usually patterns observed in solutions of a given problem, potentially of  smaller sizes. We experimented with various
  streamlining constraints, such as enforcing that each type 0, type 1, and type 2 term occurs either zero times or twice in a scheme
  (instead of an even number of times).  The most effective streamlining constraint that we came up with was observed in the Smirnov scheme:
  for each summand that is assigned a single type 3 term, enforce that (i) one matrix has either two rows, two columns
  or a row and a column fully assigned to zero and (ii) another matrix has two rows and two columns assigned to zero,
  i.e., the matrix has a single nonzero entry. This streamlining constraint reduced the runtime from minutes to seconds. Yet
  some random pairings may only be extended to a valid scheme that does not satisfy the streamlining constraint.



 §.§ Neighborhood Search


The second method is based on neighborhood search: we select a scheme, randomly fix some the corresponding base variables,
and search for an assignment for the remaining base variables. This simple method turned out to be remarkably effective to find
new schemes. The only parameter for this method is the number of base variables that will be fixed. The lower the number of
fixed base variables, the higher the probability to find a different scheme and the higher the costs to find an assignment for the remaining base variables.
We experimented with various values and it turned out that fixing 2/3 of the 621 base variables (414) is effective to
find many new schemes in a reasonable time.

The neighborhood search is able to find many new schemes, but in almost all cases they have the same pairing
of type 3 terms. Only in some rare cases the pairing of type 3 terms is different. Figure <ref>
shows such an example: scheme A has term a_13b_31c_11 in summand 22 and
term a_23b_33c_32 in summand 23, while the neighboring scheme B has term a_13b_33c_31 in summand 22
and terms a_13b_31c_11, a_23b_33c_32, and a_13b_33c_31 in summand 23.



 

§ EVALUATION AND ANALYSIS


The methods presented in Section <ref> enabled us to find several hundreds of solutions individually, but they were particularly
effective when combined.
The first method allows finding schemes that can be quite different compared to the known schemes. However, finding
a scheme using that method may require a few CPU hours as most pairings of type 3 terms cannot be extended to a valid scheme that
satisfies the streamlining constraints. The second method can find schemes that are very similar to known ones with a second. The
neighborhood of known solutions turned out to be limited to a few hundred of new schemes.

In contrast, some of the schemes found using the first method have a large neighborhood. We approximated the size of the neighborhood
of a scheme using the following experiment: Start with a given scheme S and find a neighboring scheme by randomly fixing 2/3 of the base variables. Once
a neighboring scheme S' is found, find a neighboring scheme of S', etc. We ran this experiment on a machine with 48 cores of the Lonestar 5 cluster of
Texas Advanced Computing Center. We started each experiment using 48 threads with each thread assigned a different seed.
Figure <ref> shows the number of different schemes (after sorting) found in 1000 seconds when
starting with one of the four known schemes and a scheme that was computed from the streamlining method.
Some of these different schemes are new, while others are equivalent to each other or known ones. We only assure here that they are not identical after sorting the summands.

Observe that the number of different schemes found in 1000 seconds depends a lot on the starting scheme. No different neighboring scheme was found for Laderman's scheme, only 9 different schemes were found for the scheme of Courtois et al., 94 different schemes were found for the scheme of Oh et al., 561 new schemes were found for Smirnov's scheme, and 3359 schemes were found using a randomly selected new scheme obtained with the streamlining method.



In view of the large number of solutions we found, it is also interesting to compare them with each other. For example,
if we define the support of a solution as the number of base variables set to 1, we observe that the support
seems to follow a normal distribution with mean around 160, see Fig. <ref> for a histogram. We can also see
that the Laderman scheme differs in many ways from all the other solutions. It is, for example, the only scheme whose
core consists of four quadruples of type 3 terms. In 89% of the solutions, the core consists of four pairs of type 3 terms,
about 10% of the solution have three pairs and one quadrupel, and less than 1% of the schemes have cores of the form 2-2-2-2-3
or 2-2-2-3-4.




 

§ CHALLENGES



The many thousands of new schemes
 that we found may still be just the tip of the iceberg. However, we also observed that the state-of-the-art SAT solving techniques are unable to answer
 several other questions. This section provides four challenges for SAT solvers with increasing difficulty. For each challenge we constructed one or
 more formulas that are available at 


<https://github.com/marijnheule/matrix-challenges>. 


The challenges are hard, but they may be
 doable in the coming years.



  
Challenge 1: Local search without streamlining.

Our first method combines randomly pairing the type 3 terms with streamlining constraints. The latter was required to
limit the search. We expect that local search solvers can be optimized to efficiently solve the formulas without the streamlining
constraints. This may result in schemes that are significantly different compared to ones we found. We prepared ten
satisfiable formulas with hardcoded pairings of type 3 terms. Five of these formulas can be solved using yalsat in a few
minutes. All of these formulas appear hard for CDCL solvers (and many local search solvers).



  
Challenge 2: Prove unsatisfiability of subproblems.

We observed that complete SAT solvers performed weakly on our
matrix multiplication instances. It seems therefore unlikely that one could prove any optimality results for
the product of two 3 × 3 matrices using SAT solvers in the near future. A more realistic challenge concerns proving
unsatisfiability of some subproblems. We prepared ten formulas with 23 multiplications and hardcoded pairings of type 3 terms.
We expect that these formulas are unsatisfiable.



  
Challenge 3: Avoiding a type 3 term in a summand.

All known schemes have the following property: each summand has at least one type 3 term.
We do not know whether there exists a scheme with 23 multiplications such that one of the
summands contains no type 3 term. The challenge problem blocks the existence of a type 3 term
in the last summand and does not have any additional (streamlining) constraints.



  
Challenge 4: Existence of a scheme with 22 multiplications.

The main challenge concerns finding a scheme with only 22 multiplications. The hardness
of this challenge strongly depends on whether there exists such a scheme. The repository
contains a plain formula for a scheme with 22 multiplications.



  
*Acknowledgments.
The authors acknowledge the Texas Advanced Computing Center at The University of Texas at Austin for providing HPC
resources that have contributed to the research results reported within this paper.

plain
 

