Opened 7 years ago

Last modified 3 years ago

#12177 new enhancement

Prime slicing for matrix multiplication

Reported by: SimonKing Owned by: jason, was
Priority: major Milestone: sage-feature
Component: linear algebra Keywords: prime slicing, Karatsuba
Cc: malb, burcin, robertwb, boothby Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by malb)

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

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

Attachments (3)

toom.py (779 bytes) - added by malb 7 years ago.
proof of concept
toom_matrix.py (1.2 KB) - added by SimonKing 7 years ago.
Another proof of concept
trac_12177-coeff_matrices_template.patch (29.6 KB) - added by burcin 7 years ago.

Download all attachments as: .zip

Change History (18)

comment:1 Changed 7 years ago by SimonKing

  • Description modified (diff)

comment:2 Changed 7 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 7 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 7 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 7 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.

Changed 7 years ago by malb

proof of concept

comment:6 Changed 7 years ago by malb

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

comment:7 follow-up: Changed 7 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 7 years ago by SimonKing

Replying to burcin:

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 7 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 7 years ago by SimonKing

Another proof of concept

comment:10 Changed 7 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: Changed 7 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 7 years ago by SimonKing

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?

Changed 7 years ago by burcin

comment:13 Changed 7 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 7 years ago by malb

  • Description modified (diff)

comment:15 Changed 3 years ago by kedlaya

See #9888 for a related discussion.

Note: See TracTickets for help on using tickets.