Opened 10 years ago
Last modified 5 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 )
In Sage, matrix arithmetic over finite fields is fast in the following cases:
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 overGF(p)
(with Linbox as backend) - Matrix multiplication over
GF(p^n)
is implemented by a series of multiplications overGF(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)
Attachments (3)
Change History (18)
comment:1 Changed 10 years ago by
- Description modified (diff)
comment:2 Changed 10 years ago by
comment:3 Changed 10 years ago by
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
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
Sorry, in my previous post, I totally forgot to include "reduction modulo defining polynomial". Anyway, I am now trying to produce some code.
comment:6 Changed 10 years ago by
I attached a proof of concept that the Toom thingy works.
comment:7 follow-up: ↓ 8 Changed 10 years ago by
- 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(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^{2} 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
Replying to burcin:
That's n^{2} 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
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.
comment:10 Changed 10 years ago by
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
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
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 10 years ago by
comment:13 Changed 10 years ago by
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
- Description modified (diff)
comment:15 Changed 5 years ago by
See #9888 for a related discussion.
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:
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:
I'm not sure we can do much about it, since the performance essentially depends on the speed of mod-p matrices.