The LinearAlgebra:-Generic package provides generic implementations of algorithms for linear algebra over fields, Euclidean domains, integral domains and rings.

We show below an example of how to use the MatrixInverse operation in the package to invert a Matrix A over a field.  The calling sequence is

MatrixInverse[F](A);

The indexed parameter F is called the domain of computation. For the matrix inverse operation, it must define the operations of a field, namely, addition, subtraction, multiplication and division of elements of F. It must also define the 0 and 1 of the field and a boolean procedure for testing if two elements are equal. To achieve this, the parameter F must be a table or module with the following procedures and constants:

F[`0`] : a constant, the additive identity of F

F[`1`] : a constant, the multiplicative identity of F

F[`+`] : a procedure for adding zero or more elements of F

F[`-`] : a procedure for negating or subtracting elements of F

F[`*`] : a procedure for multiplying two elements of F

F[`/`] : a procedure for dividing two elements of F

F[`=`] : a boolean procedure for comparing two elements of F

Implementation of F for Q, the field of rational numbers, is fairly straightforward:

F[`0`] := 0:

F[`1`] := 1:

F[`=`] := (x,y)->evalb( x=y ):

F[`+`] := `+`:

F[`-`] := `-`:

F[`*`] := `*`:

F[`/`] := `/`:

To use F we do

with(LinearAlgebra:-Generic):

A := Matrix([[1,2,3],[2,1,2],[3,2,1]]);

$A≔\left[\begin{array}{ccc}1& 2& 3\\ 2& 1& 2\\ 3& 2& 1\end{array}\right]$

MatrixInverse[F](A);

$\left[\begin{array}{ccc}-\frac{3}{8}& \frac{1}{2}& \frac{1}{8}\\ \frac{1}{2}& \mathrm{-1}& \frac{1}{2}\\ \frac{1}{8}& \frac{1}{2}& -\frac{3}{8}\end{array}\right]$

The code for MatrixInverse is generic.  It can be used to invert a matrix over any field, we just need to define the required operations. Let's do it for a finite field.  We will represent elements of GF(27) as polynomials in z over the integers modulo 3.  That is, we represent elements of the field as polynomials modulo m(z), an irreducible polynomial of degree 3 in z.  We will assume that elements of the field are reduced, i.e., have degree less than 3.

p := 3;

$p≔3$

m := z^3+2*z+1;

$m≔z^3+2z+1$

GF27[`0`] := 0:

GF27[`1`] := 1:

GF27[`=`] := `=`:

GF27[`+`] := ()->modp(`+`(args),p):

GF27[`-`] := (a,b)->`if`(nargs=1,modp(-a,p),modp(a-b,p)):

GF27[`*`] := (a,b)->modp(Rem(a*b,m,z),p):

GF27[`/`] := proc(a,b) local i;
               if b=0 then
                 error "division by zero";
               end if;
               Gcdex(b,m,z,'i') mod p;
               Rem(a*i,m,z) mod p;
             end proc:

A := Matrix([[1,z,z^2],[z,1,z],[z^2,z,1]]);

$A≔\left[\begin{array}{ccc}1& z& z^2\\ z& 1& z\\ z^2& z& 1\end{array}\right]$

B := MatrixInverse[GF27](A);

$B≔\left[\begin{array}{ccc}z& 2z^2& 0\\ 2z^2& 2z+2& 2z^2\\ 0& 2z^2& z\end{array}\right]$

MatrixMatrixMultiply[GF27](A,B);

$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

In the above example we demonstrated the use of the MatrixMatrixMultiply routine to verify that the computed inverse really is the inverse.

List of LinearAlgebra:-Generic Package Commands

The package exports algorithms and commands for linear algebra which are parameterized by the ring in which the matrix entries lie. The package also exports algorithms.  For example, the Bareiss fraction-free algorithm is well known.  It is an algorithm which computes the determinant of a matrix over any integral domain D. It assumes exact division and does O(n^3) operations from D. The Berkowitz algorithm is division free and takes O(n^4) operations. The output Matrix B over D satisfies certain properties, in particular, B[i,i] is the determinant (up to sign) of the principal i x i minor. So one can compute the determinant using this algorithm by doing

Determinant[D](A,method=BareissAlgorithm);

BareissAlgorithm[D](A);

The domain may be a Maple table or a Maple module.

A common mistake is to forget to specify the domain of computation. If you do this you will get an error.  For example:

MatrixInverse(A); # should be MatrixInverse[GF27](A);

Error, (in LinearAlgebra:-Generic:-GenericCheck) LinearAlgebra:-Generic:-MatrixInverse is not indexed by a domain

Another common error is to omit one of the operations in the domain. The package checks that the required operations are defined and of type procedure (or constant).  Here is an example of an error

HermiteForm[GF27](A);

Error, (in LinearAlgebra:-Generic:-GenericCheck) missing operation: Gcdex

The required operations for the domain of computation F for each algorithm are specified in the help page for that algorithm. The package been designed to keep this list minimal. Note, even if the algorithm requires a field, it may not use all of the basic operations from a field, for example, the algorithm for Gaussian elimination does not use addition.  Nevertheless, addition is a required operation.

The additive identity 0 may or may not be unique.  The algorithms will always use the following to test for zero:

if F[`=`](x,F[`0`]) then ... else ... end if;

Other elements of a ring, including 1, do not need a unique representation.

$\mathrm{with}\left(\mathrm{LinearAlgebra}:-\mathrm{Generic}\right)\:$

$R≔\mathbf{module}\left(\right)\mathrm{\_export}\left(\mathrm{`0`}\,\mathrm{`1`}\,\mathrm{`=`}\,\mathrm{`+`}\,\mathrm{`-`}\,\mathrm{`*`}\right)\;\mathrm{`0`}≔0\;\mathrm{`1`}≔1\;\mathrm{`=`}≔\left(x\,y\right)\→\mathrm{evalb}\left(x\=y\right)\;\mathrm{`+`}≔\left(\right)\→\mathrm{modp}\left(:-\mathrm{`+`}\left(\mathrm{args}\right)\,n\right)\;\mathrm{`*`}≔\left(a\,b\right)\→\mathrm{modp}\left(a\*b\,n\right)\;\mathrm{`-`}≔\left(a\,b\right)\→\mathrm{`if`}\left(\mathrm{nargs}\=1\,\mathrm{modp}\left(\−a\,n\right)\,\mathrm{modp}\left(a-b\,n\right)\right)\mathbf{end\; module}$

$R≔\mathbf{module}\left(\right)\mathbf{export}0\,1\,\mathrm{`=`}\,\mathrm{`+`}\,\mathrm{`-`}\,\mathrm{`*`}\;\mathbf{end\; module}$

$n≔6\:$

$R\left[\mathrm{`*`}\right]\left(2\,4\right)$

$2$

$A≔\mathrm{Matrix}\left(\left[\left[2\,3\,4\,5\right]\,\left[3\,4\,1\,2\right]\,\left[4\,1\,1\,5\right]\,\left[3\,1\,5\,4\right]\right]\right)\:$

$\mathrm{Determinant}\left[R\right]\left(A\,\mathrm{method}=\mathrm{BerkowitzAlgorithm}\right)$

$3$

$Z\left[\mathrm{`0`}\right],Z\left[\mathrm{`1`}\right]≔0,1$

$Z_0,Z_1≔0,1$

$Z\left[\mathrm{`+`}\right],Z\left[\mathrm{`-`}\right],Z\left[\mathrm{`*`}\right],Z\left[\mathrm{`=`}\right]≔\mathrm{`+`},\mathrm{`-`},\mathrm{`*`},\mathrm{`=`}$

$Z_{\mathrm{`+`}},Z_{\mathrm{`-`}},Z_{\mathrm{`*`}},Z_{\mathrm{`=`}}≔\mathrm{`+`},\mathrm{`-`},\mathrm{`*`},\mathrm{`=`}$

$Z\left[\mathrm{Divide}\right]≔\left(x\,y\right)\mapsto \mathrm{evalb}\left(\mathrm{irem}\left(x\,y\,\mathrm{args}\left[3..\mathrm{nargs}\right]\right)=0\right)\:$

$A≔\mathrm{Matrix}\left(\left[\left[2\,3\,4\,5\right]\,\left[3\,4\,1\,2\right]\,\left[4\,1\,1\,5\right]\,\left[3\,1\,5\,4\right]\right]\right)\:$

$\mathrm{Determinant}\left[Z\right]\left(A\,\mathrm{method}=\mathrm{BareissAlgorithm}\right)$

$195$

$\mathrm{BareissAlgorithm}\left[Z\right]\left(A\right)$

$\left[\begin{array}{cccc}2& 3& 4& 5\\ 0& \mathrm{-1}& \mathrm{-10}& \mathrm{-11}\\ 0& 0& \mathrm{-43}& \mathrm{-50}\\ 0& 0& 0& 195\end{array}\right]$

ANNullSpace := proc(AA::Matrix,rof::RootOf)
local M,z,R,m,n,i,j,A,F,B;
   M := subs(_Z=z,op(1,rof));
   R := polynom(rational,z);
   if not type(M,R) then
     error "invalid RootOf";
   end if;
   m,n := LinearAlgebra:-Dimensions(AA);
   A := Matrix(m,n);
   for i to m do
     for j to n do
       A[i,j] := subs(rof=z,AA[i,j]);
       if not type(A[i,j],R) then
         error "invalid matrix entries";
       end if;
       A[i,j] := rem(A[i,j],M,z);
     end do;
   end do;
   F[`0`] := 0; F[`1`] := 1;
   F[`=`] := `=`; F[`+`] := `+`; F[`-`] := `-`;
   F[`*`] := (a,b)->rem(a*b,M,z);
   F[`/`] := proc(a,b)
             local i;
               if b=0 then error "division by zero";
               elif gcdex(b,M,z,'i')<>1 then
                 error "invalid RootOf";
               else
                 rem(a*i,M,z);
               end if;
             end proc:
   B := LinearAlgebra:-Generic:-NullSpace[F](A);
   subs( z=rof, B );
end proc:

$a≔\mathrm{RootOf}\left(z^3-2\right)$

$a≔\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)$

$M≔\mathrm{Matrix}\left(\left[\left[1\&comma;a\&comma;a^2\right]\&comma;\left[a+a^2\&comma;a+a^2\&comma;3\right]\&comma;\left[a^2\&comma;a\&comma;1\right]\right]\right)$

$M≔\left[\begin{array}{ccc}1& \mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)& {\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)}^2\\ {\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)}^2+\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)& {\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)}^2+\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)& 3\\ {\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)}^2& \mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)& 1\end{array}\right]$

$B≔\mathrm{ANNullSpace}\left(M\&comma;a\right)$

$B≔\left\{\left[\begin{array}{c}1\\ -\frac{{\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)}^2}{2}-\mathrm{RootOf}\left({\mathrm{\_Z}}^3-2\right)\\ 1\end{array}\right]\right\}$

$C\left[\mathrm{`0`}\right]≔0\&colon;$

$C\left[\mathrm{`1`}\right]≔1\&colon;$

$C\left[\mathrm{`=`}\right]≔\mathrm{`=`}\&colon;$

$C\left[\mathrm{`+`}\right]≔\mathrm{`+`}\&colon;$

$C\left[\mathrm{`-`}\right]≔\mathrm{`-`}\&colon;$

$C\left[\mathrm{`*`}\right]≔\mathrm{`*`}\&colon;$

$C\left[\mathrm{`/`}\right]≔\mathrm{`/`}\&colon;$

$A≔\mathrm{Matrix}\left(\left[\left[1\&comma;\mathrm{I}\&comma;-\mathrm{I}\right]\&comma;\left[3\&comma;-1+\mathrm{I}\&comma;-4\mathrm{I}\right]\&comma;\left[\mathrm{I}\&comma;1-\mathrm{I}\&comma;2\right]\right]\right)$

$A≔\left[\begin{array}{ccc}1& \mathrm{I}& \mathrm{-I}\\ 3& \mathrm{-1}+\mathrm{I}& -4\mathrm{I}\\ \mathrm{I}& 1-\mathrm{I}& 2\end{array}\right]$

$B≔\mathrm{NullSpace}\left[C\right]\left(A\right)$

$B≔\left\{\left[\begin{array}{c}-\frac{1}{5}+\frac{7\mathrm{I}}{5}\\ -\frac{2}{5}-\frac{\mathrm{I}}{5}\\ 1\end{array}\right]\right\}$

MatrixTimesVector := proc(A::Matrix,v::Vector)
local D,n,p,m,C,i,k;
  D := GenericCheck( procname, [`0`,`+`,`*`] );
  n,p := op(1,A);
  m := op(1,v);
  if p<>m then
    error "incompatible dimensions";
  end if;
  C := Vector(n,'fill'=D[`0`]);
  for i to n do
    for k to p do
      C[i] := D[`+`](C[i],D[`*`](A[i,k],v[k]));
    end do;
  end do;
  C;
end proc:

$\mathrm{MatrixTimesVector}\left[C\right]\left(A\&comma;B\left[1\right]\right)$

$\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$