Opened 10 years ago

# Prime slicing for matrix multiplication — at Version 14

Reported by: Owned by: SimonKing jason, was major sage-feature linear algebra prime slicing, Karatsuba malb, burcin, robertwb, boothby N/A

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:

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

• 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)

### comment:1 Changed 10 years ago by SimonKing

• Description modified (diff)

### comment:2 Changed 10 years ago by malb

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.

### comment:3 Changed 10 years ago by SimonKing

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.

### comment:4 Changed 10 years ago by SimonKing

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.

### comment:5 Changed 10 years ago by SimonKing

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

proof of concept

### comment:6 Changed 10 years ago by malb

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

### comment:7 follow-up: ↓ 8 Changed 10 years ago by burcin

• Milestone changed from sage-4.8 to sage-feature

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(pk) 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 n2 coefficient matrix multiplications with some overhead to handle the rollover with the minimal polynomial. We should look at a better algorithm. :)

### comment:8 in reply to: ↑ 7 Changed 10 years ago by SimonKing

That's n2 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?

### comment:9 Changed 10 years ago by malb

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.

### Changed 10 years ago by SimonKing

Another proof of concept

### comment:10 Changed 10 years ago by SimonKing

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 toom_matrix.py, which returns the equivalent of `FWD` and `BCK` - perhaps you can plug it into your code?

### comment:11 follow-up: ↓ 12 Changed 10 years ago by 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.

### comment:12 in reply to: ↑ 11 Changed 10 years ago by SimonKing

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?

### comment:13 Changed 10 years ago by burcin

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>
```

### comment:14 Changed 10 years ago by malb

• Description modified (diff)
Note: See TracTickets for help on using tickets.