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`;
    With[{ext = 12},
               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!


2 Answers 2


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

 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]];


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

(* True *)

--- end edit ---

  • 1
    $\begingroup$ Thank you for great references! It seems that I has to reimplement finite fields stuff in Mathematica to improve performance. This is bad news for me :( $\endgroup$ Commented Nov 25, 2012 at 11:28
  • $\begingroup$ I think the key is Algebra`PolynomialExtendedGCDModList because it allows you to find field inverses. You invoke it using polynomials represented as their coefficient lists (constant term on the left and ascending from there). If defpoly is the list for your field defining polynomial, then to invert the field element r1 you would do r1inv=Algebra`PolynomialExtendedGCDModList[r1,defpoly,p][[2,1]]. In this notation p is the prime charatceristic, so it would be 2 for the cases of apparent interest. $\endgroup$ Commented Nov 26, 2012 at 19:04
  • $\begingroup$ @Daniel Lichtblau Thank you for the point! I see, the easiest ways are to implement finite field extensions myself or do the patch for FiniteFields. If I get some progress, I will share it in the community. $\endgroup$ Commented Nov 27, 2012 at 10:11
  • $\begingroup$ Yes, please do share with community. Here and/or library.wolfram.com would be good places. Is your work only with linear algebra on field elements, or does it involve e.g. polynomials over extension fields? if the latter, there may be some existing functionality of interest I could locate. $\endgroup$ Commented Nov 27, 2012 at 15:35
  • $\begingroup$ @DanielLichtblau My work includes both, the linear algebra and polynomials over extension fields. For example, I have to compute the determinant of big matrices over field extensions and to search roots of polynomials by their factorization. Sure, it will be great! Thanks. $\endgroup$ Commented Dec 3, 2012 at 12:32

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}

  • $\begingroup$ Works well over prime fields but not extensions thereof. $\endgroup$ Commented May 27, 2013 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.