Prime slicing for matrix multiplication

[simon-king-jena]

In Sage, matrix arithmetic over finite fields is fast in the following cases:

• Prime fields, using linbox (#4260) resp. M4RI over GF(2)
• GF(2^e), using M4RIE (#9562)

In all other cases, it sucks. There is the suggestion to use a wrapper of a fork of C-MeatAxe (#12103), but this would only work up to field size 255 and might have further disadvantages.

Martin Albrecht suggested to use "prime slicing" instead:

• Represent matrices over GF(p^n) by a list of n matrices over GF(p) (with Linbox as backend)
• Matrix multiplication over GF(p^n) is implemented by a series of multiplications over GF(p)
• Karatsuba type formulae reduce the number of "prime" multiplications involved.

On sage-devel, Martin gave a couple of references:

• Boothby and Bradshaw Bitslicing and the Method of Four Russians Over Larger Finite Fields
• Montgomery Five, Six, and Seven-Term Karatsuba-Like Formulae

On another occasion, Martin also pointed out how to compute echelon forms in that setting.

The purpose of this ticket is to make that idea real.

TODO

• addition, subtraction (easy)
• scalar multiplication (medium)
• row swaps, column swaps (easy)
• C = C + A*B (easy)
• fast randomize (for testing, easy)
• apply LAPACK permutations (easy)
• matrix windows (cheap submatrices, arder)
• TRSM lower left / TRSM upper left (medium)
• PLE (medium)
• Karatsuba-like for up to degree 10 (harder)

CC: @malb @burcin @robertwb @boothby

Component: linear algebra

Keywords: prime slicing, Karatsuba

Issue created by migration from https://trac.sagemath.org/ticket/12177

[simon-king-jena]

Description changed:

--- 
+++ 
@@ -14,8 +14,7 @@
 On [sage-devel](http://groups.google.com/group/sage-devel/browse_thread/thread/8a988a2f5d997a5a/ebcde393f8ef5240), Martin gave a couple of references:

 * [Boothby and Bradshaw](http://www.google.com/url?sa=D&q=http://arxiv.org/abs/0901.1413&usg=AFQjCNFg5cUrXjBTqkon1V6n4Yl_CVQ9yA) Bitslicing and the Method of Four Russians Over Larger Finite Fields
-* [Montgomery](http://bit.ly/r5PVgv) Five, Six, and Seven-Term
-Karatsuba-Like Formulae
+* [Montgomery](http://bit.ly/r5PVgv) Five, Six, and Seven-Term Karatsuba-Like Formulae

 On another occasion, Martin also pointed out how to compute echelon forms in that setting.

[malb]

comment:2

Hi Simon, great that you're getting the ball rolling. Two comments:

(a) I think the this code shouldn't be in Sage but on a lower level, preferably LinBox. But perhaps we can write proof of concept in Sage first and then port to C++?

(b)I figure we should get an idea about what performance to expect, so:

p = 17
K.<a> = GF(p^2)
A = [random_matrix(GF(p),2000,2000) for _ in range(2)]
B = [random_matrix(GF(p),2000,2000) for _ in range(2)]
C = [matrix(GF(p),2000,2000) for _ in range(2)]

def mul2(K, C,A,B):
   poly = K.polynomial()
   T0 = A[0]*B[0]
   C[0] += T0
   C[1] -= T0
   C[1] += (A[0] + A[1])*(B[0]+B[1])
   T0 = A[1]*B[1]
   C[1] -= T0
   for i,c in enumerate(list(poly)):
       if i == 2:
           break
       C[i] -= c*T0
   return C

mul2 (I didn't check for bugs) with these inputs takes about 3 seconds on my computer, which was to be expected since a product over GF(p) takes about 1 second on my computer. For comparison Magma is a bit better:

Magma V2.15-10    Sun Dec 18 2011 12:53:05 on road     [Seed = 1273554698]
Type ? for help.  Type <Ctrl>-D to quit.
> K:=GF(17^2);
> A:=RandomMatrix(K,2000,2000);
> B:=RandomMatrix(K,2000,2000);
> time C:=A*B;
Time: 2.320

I'm not sure we can do much about it, since the performance essentially depends on the speed of mod-p matrices.

[simon-king-jena]

comment:3

I think we have to do two things:

. Have some code that deals with lists of matrices. I understand that Burcin has such code.

. Find a way to find a Karatsuba type formula that reduces the number of mod-p multiplications.

I think the second point should actually done not by Karatsuba, but by a modification of Level-n Toom multiplication (where n is the degree of our field extension). See Wikipedia for a method to compute the level-n formula.

Ideally, we would also merge the defining polynomial of our field extension into the Toom formula.

[simon-king-jena]

comment:4

Here are some thoughts - I am not sure if this is all correct, but I think it is a starting point:

Elements of K=GF(p^n) are represented by polynomials of degree n-1 over k=GF(p).

Assumption: p>>n.

Each polynomial of degree n-1 (hence, any element of K) is determined by evaluation at n points. When multipliying two polynomials then (modulo the defining polynomial of K) they also give rise to a polynomial of degree n-1. Hence, again, n evaluation points are enough.

I am not sure, though, whether we can choose any n-tuple of elements of k.

First step:

Choose n points from k. Evalutation of a generic polynomial of degree n-1 gives rise to n linear combinations of the polynomial's coefficients, described by some invertible square matrix. Let A be the inverse of that matrix.

Multiplying two polynomials of degree n modulo the defining polynomial now means: Evaluate the two factors at the given n points. For each point, multiply the two evaluations. Multiply this with A, and read off the coefficients of the product.

So, this can be used for our polynomials of matrices. If my arguments work, then we would be down to n multiplications over GF(p) for one multiplication over GF(p^n). And comparing the figures of Linbox over GF(17) and Magma over GF(17^2), we would be competitive with Magma.

[simon-king-jena]

comment:5

Sorry, in my previous post, I totally forgot to include "reduction modulo defining polynomial". Anyway, I am now trying to produce some code.

[malb]

Attachment: toom.py.gz

proof of concept

[malb]

comment:6

I attached a proof of concept that the Toom thingy works.

[burcin]

comment:7

I attached a patch with template classes implementing this polynomial with matrix coefficients representation. There is also an implementation of the naive multiplication algorithm for GF(p^k) to demonstrate FFLAS calls.

Timings for the naive multiplication from Martin's example above (comment:2) for p=17 are:

k     time
2     4.51    
3    10.28
4    19.17

That's n² coefficient matrix multiplications with some overhead to handle the rollover with the minimal polynomial. We should look at a better algorithm. :)

[simon-king-jena]

comment:8

Replying to @burcin:

That's n² coefficient matrix multiplications with some overhead to handle the rollover with the minimal polynomial. We should look at a better algorithm. :)

... which probably is Toom. Martin, do you have code for Toom multiplication?

[malb]

comment:9

I've attached it, but it's going to be slow if you try it because a*A for a in the field there is no special code in Sage yet.

[simon-king-jena]

Attachment: toom_matrix.py.gz

Another proof of concept

[simon-king-jena]

comment:10

Martin, in your Toom proof of concept, aren't you evaluating at too few points? We start with polynomials of degree less than the degree of the field extension. But the product will have (before reduction) twice that degree.

Hence, instead of FWD = Matrix(K, l, l, [K(i**j) for i in range(l) for j in range(l)]), I think we need

    FWD = Matrix(K, l, l, [K(i**j) for i in range(l) for j in range(2*l-1)])

The following lines are to be changed accordingly.

We could of course include the reduction to degree l-1 into the matrix BCK. I did something along these lines in attachment: toom_matrix.py, which returns the equivalent of FWD and BCK - perhaps you can plug it into your code?

[malb]

comment:11

Simon, l = len(A)+len(B)-1 i.e., l is not the degree of the inputs. Merging the modular reduction with the interpolation might be a good idea though.

[simon-king-jena]

comment:12

Replying to @malb:

Simon, l = len(A)+len(B)-1 i.e., l is not the degree of the inputs. Merging the modular reduction with the interpolation might be a good idea though.

Sorry! I was somehow misreading it as l = len(A) - what happened to my eyes?

[burcin]

Attachment: trac_12177-coeff_matrices_template.patch.gz

[burcin]

comment:13

I refreshed my patch with a version that implements the multi point evaluation approach using FFLAS. It can be reached via the _multiply_toom() function:

sage: K.<a> = GF(17^6)
sage: MS = MatrixSpace(K, 2000, 2000)
sage: from sage.matrix.matrix_modq_dense_float import Matrix_modq_dense_float
sage: M = Matrix_modq_dense_float(MS, a^5)
sage: res = M._multiply_toom(M)
<some debugging output>

[malb]

Description changed:

--- 
+++ 
@@ -19,3 +19,16 @@
 On another occasion, Martin also pointed out how to compute echelon forms in that setting.

 The purpose of this ticket is to make that idea real.
+
+**TODO**
+
+* addition, subtraction (easy)
+* scalar multiplication (medium)
+* row swaps, column swaps (easy)
+* C = C + A*B (easy)
+* fast randomize (for testing, easy)
+* apply LAPACK permutations (easy)
+* matrix windows (cheap submatrices, arder)
+* TRSM lower left / TRSM upper left (medium)
+* PLE (medium)
+* Karatsuba-like for up to degree 10 (harder)

[kedlaya]

comment:15

See #9888 for a related discussion.