 LinearAlgebra - Maple Programming Help

Home : Support : Online Help : Mathematics : Linear Algebra : LinearAlgebra Package : Generic Subpackage : LinearAlgebra/Generic

LinearAlgebra

 Generic
 Overview of the LinearAlgebra[Generic] Package

Description

 • 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 : a constant, the additive identity of F
 F : 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:
 > F := 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]$ (1)
 > MatrixInverse[F](A);
 $\left[\begin{array}{ccc}{-}\frac{{3}}{{8}}& \frac{{1}}{{2}}& \frac{{1}}{{8}}\\ \frac{{1}}{{2}}& {-1}& \frac{{1}}{{2}}\\ \frac{{1}}{{8}}& \frac{{1}}{{2}}& {-}\frac{{3}}{{8}}\end{array}\right]$ (2)
 • 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}$ (3)
 > m := z^3+2*z+1;
 ${m}{≔}{{z}}^{{3}}{+}{2}{}{z}{+}{1}$ (4)
 > GF27 := 0:
 > GF27 := 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]$ (5)
 > B := MatrixInverse[GF27](A);
 ${B}{≔}\left[\begin{array}{ccc}{z}& {2}{}{{z}}^{{2}}& {0}\\ {2}{}{{z}}^{{2}}& {2}{}{z}{+}{2}& {2}{}{{z}}^{{2}}\\ {0}& {2}{}{{z}}^{{2}}& {z}\end{array}\right]$ (6)
 > MatrixMatrixMultiply[GF27](A,B);
 $\left[\begin{array}{ccc}{1}& {0}& {0}\\ {0}& {1}& {0}\\ {0}& {0}& {1}\end{array}\right]$ (7)
 • 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);

and also, one can compute the Matrix B by doing

 > BareissAlgorithm[D](A);

Notes on the definition of the domains

 • 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);
 • 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);
 • 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) then ... else ... end if;
 • Other elements of a ring, including 1, do not need a unique representation.

Examples

Suppose we want to compute the determinant of a matrix of integers modulo a composite integer n using the Berkowitz algorithm.  We need to first define the ring of integers modulo n as a domain.

 > $\mathrm{with}\left(\mathrm{LinearAlgebra}\left[\mathrm{Generic}\right]\right):$
 > $R≔\mathbf{module}\left(\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathrm{_export}\left(\mathrm{0},\mathrm{1},\mathrm{=},\mathrm{+},\mathrm{-},\mathrm{*}\right);\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathrm{0}≔0;\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathrm{1}≔1;\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathrm{=}≔\left(x,y\right)→\mathrm{evalb}\left(x=y\right);\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathrm{+}≔\left(\right)→\mathrm{modp}\left(:-\mathrm{+}\left(\mathrm{args}\right),n\right);\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathrm{*}≔\left(a,b\right)→\mathrm{modp}\left(a*b,n\right);\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\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)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathbf{end module}$
 ${R}{:=}{\mathbf{module}}\left({}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{export}}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{0}{,}{1}{,}{\mathrm{=}}{,}{\mathrm{+}}{,}{\mathrm{-}}{,}{\mathrm{*}}{;}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end module}}$ (8)
 > $n≔6:$
 > $R\left[\mathrm{*}\right]\left(2,4\right)$
 ${2}$ (9)
 > $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}$ (10)

A more efficient algorithm is the Bareiss fraction free algorithm. It assumes only exact division and computes the determinant in O(n^3) operations in R.  Hence it is valid for the integers.  We illustrate this time using a table Z of operations for the ring.

 > $Z\left[\mathrm{0}\right],Z\left[\mathrm{1}\right]≔0,1$
 ${{Z}}_{{0}}{,}{{Z}}_{{1}}{≔}{0}{,}{1}$ (11)
 > $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{=}}$ (12)
 > $Z\left[\mathrm{Divide}\right]≔\left(x,y\right)↦\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}$ (13)
 > $\mathrm{BareissAlgorithm}\left[Z\right]\left(A\right)$
 $\left[\begin{array}{cccc}{2}& {3}& {4}& {5}\\ {0}& {-1}& {-10}& {-11}\\ {0}& {0}& {-43}& {-50}\\ {0}& {0}& {0}& {195}\end{array}\right]$ (14)

Given a Matrix A of algebraic numbers which were input using Maple's RootOf notation, suppose we want to compute the nullspace of A. We will convert to a univariate polynomial representation for the computation to obtain greater efficiency.

 > 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; F := 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)$ (15)
 > $M≔\mathrm{Matrix}\left(\left[\left[1,a,{a}^{2}\right],\left[a+{a}^{2},a+{a}^{2},3\right],\left[{a}^{2},a,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]$ (16)
 > $B≔\mathrm{ANNullSpace}\left(M,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\}$ (17)

For the last example, we compute the nullspace of a matrix of complex rationals.

 > $C\left[\mathrm{0}\right]≔0:$
 > $C\left[\mathrm{1}\right]≔1:$
 > $C\left[\mathrm{=}\right]≔\mathrm{=}:$
 > $C\left[\mathrm{+}\right]≔\mathrm{+}:$
 > $C\left[\mathrm{-}\right]≔\mathrm{-}:$
 > $C\left[\mathrm{*}\right]≔\mathrm{*}:$
 > $C\left[\mathrm{/}\right]≔\mathrm{/}:$
 > $A≔\mathrm{Matrix}\left(\left[\left[1,I,-I\right],\left[3,-1+I,-4I\right],\left[I,1-I,2\right]\right]\right)$
 ${A}{≔}\left[\begin{array}{ccc}{1}& {I}& {-I}\\ {3}& {-1}{+}{I}& {-}{4}{}{I}\\ {I}& {1}{-}{I}& {2}\end{array}\right]$ (18)
 > $B≔\mathrm{NullSpace}\left[C\right]\left(A\right)$
 ${B}{≔}\left\{\left[\begin{array}{c}{-}\frac{{1}}{{5}}{+}\frac{{7}{}{I}}{{5}}\\ {-}\frac{{2}}{{5}}{-}\frac{{I}}{{5}}\\ {1}\end{array}\right]\right\}$ (19)

To check that the answer is correct we need to verify that A.B = 0. We will write a procedure to do this.  This illustrates how the procedures in LinearAlgebra[Generic] are coded. The GenericCheck procedure checks that the procedure was called with a domain (a table or module) and that it has the required values/exports.

 > 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);   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,B\left[1\right]\right)$
 $\left[\begin{array}{c}{0}\\ {0}\\ {0}\end{array}\right]$ (20)