FiniteFields package is very slow. Any fast substitute for Mathematica?

Question

I want to compute the inverse of matrix, say with dimensions $100 \times 100$, defined over a large finite field extension such as $GF(2^{120})$. I am using the package FiniteFields, but Mathematica's computation time is exponential with respect to matrix dimensions. The following code illustrates the problem:

<< FiniteFields`;
Table[
    With[{ext = 12},
       First@AbsoluteTiming@
           Inverse[
               Array[GF[2, ext][RandomInteger[{0, 1}, ext]] &, {n, n}]
           ]
    ],
    {n, 1, 11}
]

I am using an Intel Xeon X5680 @ 3.33GHz (64-bit OS) and Mathematica v8.0.4.0. I have received the following timing results:

{0.0030, 0.0080, 0.0210, 0.0630, 0.1860, 0.5110, 1.3350, 3.3840, 8.9340, 23.0090, 57.4660}

I believe the source of the problem is that the FiniteFields package defines many UpValues, DownValues and SubValues of Times and Plus for head GF and, consequently, the pattern matching of arguments is degraded.

Does anyone know if a patch for the FiniteFields package or a faster substitute providing a similar interface?

Many thanks!

Answer 1

--- edit ---

More recent versions of the work containing the code mentioned below are found at these links: SymbolicFAQ and PCwGB

--- end edit ---

This will not scale to dimension 100 but will be an improvement on what you now have. It is cribbed from the section "Linear Algebra over Galois Fields here as well as the section "Groebner bases over modules and related computations" in this notebook.

deg = 12;
flen = 3;
j = 0;
While[flen > 2 && j++ < 100,
 defpoly = 
  x^deg + 1 + RandomInteger[{0, 1}, deg - 1].x^Range[1, deg - 1];
 flen = Length[FactorList[defpoly, Modulus -> 2]]
 ]

dim = 20;
mat = Table[
   RandomInteger[{0, 1}, deg].x^Range[0, deg - 1], {dim}, {dim}];

newvars = Array[z, 2*dim];
augmat = Transpose[
   ArrayFlatten[{mat, IdentityMatrix[Length[mat]]}, 1]];
auxpolys = Union[Flatten[Outer[Times, newvars, newvars]]];
linpolys = Join[augmat.newvars, {defpoly}, auxpolys];
allvars = Append[newvars, x];

Here is the bulk of the computation.

Timing[gb = GroebnerBasis[linpolys, allvars, Modulus -> 2];]

(* {168.240000, Null} *)

Now extract the inverse matrix.

modulegb = Complement[gb, Join[auxpolys, {defpoly}]];
redmat = Reverse[Sort[Outer[D, modulegb, newvars]]];
invmat = Transpose[redmat[[All, dim + 1 ;; 2*dim]]];

We'll check the result.

PolynomialMod[invmat.mat, {defpoly, 2}] === IdentityMatrix[dim]

(* True *)

A serious bottleneck is the need, using this approach, to have $O(\tt{dim}^2)$ auxiliary polynomials.

One can improve on this by writing a direct row reduction using a list representation form the field elements. This would require working with undocumented functionality in the Algebra` context, such as Algebra`PolynomialTimesModList. This TMJ article may give slightly more information should you opt to go that route.

--- edit ---

Here is an approach that will sometimes work. Treat the elements as algebraics over Q instead of GF[2]. Now you can use some built in functionality. If you are lucky you get substantially the same row reduction and can recover the mod-2 result. I'll show what I mean with the example above.

We first set up a matrix of AlgebraicNumber opjects.

alg = algnum[root[defpoly], CoefficientList[#1, x]] &;
mat2 = Map[alg, mat, {2}] /. x -> # /. root[a_] :> Root[a &, 1] /. 
   algnum -> AlgebraicNumber;
augmat = Transpose[
   ArrayFlatten[{mat2, IdentityMatrix[Length[mat2]]}, 1]];

Now row reduce it and hope we do not get any denominators divisible by 2.

Timing[
 rred = RowReduce[augmat, Method -> "OneStepRowReduction"];]

(* {36.53499999999994, Null} *)

Extract the inverse part.

invmat = Transpose[
   rred[[All, dim + 1 ;; 2*dim]] /. 
    AlgebraicNumber[aa_, bb_] :> 
     PolynomialMod[bb, 2].x^Range[0, deg - 1]];

Check:

PolynomialMod[invmat.mat, {defpoly, 2}] === IdentityMatrix[dim]

(* True *)

--- end edit ---

Answer 2

The textbook algorithm, that is, row reducing {A, I} to {I, A^(-1)} works.

n = 1000;
m = Array[RandomInteger[{0, 1}]&, {n, n}];
m // MatrixPlot

MatrixRank[m, Modulus -> 2] (* see if m is full rank *)

inverse = With[{n = Length@#},
  RowReduce[MapThread[Join, {#, IdentityMatrix[n]}], Modulus -> 2][[;;,n+1;;]]]&;
inverseM = inverse[m]; // AbsoluteTiming

Mod[inverseM. m, 2] // MatrixPlot

On my laptop, computing inverse takes 15.262331 for a full rank matrix.

Often, you use the inverse matrix to solve A x == b, and Mathematica has fast LinearSolve in the FiniteFields package.

n = 1000;
LinearSolve[Array[RandomInteger[{0, 1}]&, {n, n}],
  Array[RandomInteger[{0, 1}]&,n], Modulus -> 2] ;// AbsoluteTiming

On my laptop this last gives {5.110337, Null}