GAP Instructional Material

Return to Index

GAP code is shown like this.

[Press ?? to see GAP output]

You can also choose to see a page showing ALL GAP output.

Vectors and Matrices in GAP

Vectors over the rational numbers

We begin by looking at ways to define a vector.

v1:= [ 1/2, 3/4, -2/3, 22/7 ];      ??

We can see if what we have defined is a vector

IsVector( v1 );      ??

We define a vector space over the rational numbers

Suppose we want a three dimensional subspace of R4

v2:= [ 1, 3, 2, 4 ];      ??

v3:= [ 1/2, 1/4, 1/3, 3/4 ];      ??

V:= VectorSpace( Rationals, [ v1, v2, v3 ] );      ??

It is easy to test if vectors lie in the vector space V

[ 1, 1, 1, 1 ] in V;      ??

[ 28, 70, 84, 45 ] in V;      ??

Scalar multiplication of vectors is easy

1/2*v1;      ??

As is the dot product

v1*v2;      ??

Of course linear combinations of vectors from V lie in V. We try an example

1/2*v1 - 3/4*v2 + 5/4*v3 in V;      ??

Vectors over a finite field

Let us construct the 3 dimensional vector space over GF(3)

V:= FullRowSpace( GF( 3 ), 3 );      ??

Let us check whether a vector is in V

[ 1, 1, 1 ] in V;      ??

This gives false since 1 is not in GF(3). We need to use the identity of GF(3)

o:= One( GF( 3 ) );      ??

[ 1, 1, 1 ]*o in V;      ??

Let us look at the sum of all the vectors of V

s:= [ 0, 0, 0 ]*o;;
for v in V do
  s:= s + v;
od;
Print(s, "\n");      ??

It is easy to see why one gets 0 since for each v in V, -v is also in V. What about a field of characteristic 2?

V:= FullRowSpace( GF( 4 ), 4 );;
o:= One( GF( 4 ) );;
s:= [ 0, 0, 0, 0 ]*o;;
for v in V do
  s:= s + v;
od;
Print(s, "\n");      ??

So the answer is still 0. This is not quite so obvious - what about GF(2) with dimension 1?

Bases

We set up a vector space over the rationals spanned by 4 vectors

v1:= [ 2, 2, 1, 3 ];;
v2:= [ 7, 5, 5, 5 ];;
v3:= [ 3, 2, 2, 1 ];;
v4:= [ 2, 1, 2, 1 ];;

V:= VectorSpace( Rationals, [ v1, v2, v3, v4 ] );      ??

We compute its basis

B:= Basis( V );      ??

What is the dimension of V?

Length( B );      ??

Let us see the basis vectors GAP has computed

BasisVectors( B );      ??

We can express vectors in V as a linear combination of the basis vectors B

Coefficients( B, [2,1,2,1]);      ??

Since the coefficients are [ 2, -1, 1/4 ] we should have

[ 2, 1, 2, 1 ] = 2*B[1] - 1*B[2] + 1/4*B[3];      ??

What if we try to express a vector not in the space as a linear combination?

Coefficients( B, [1, 0, 0, 0]);      ??

Not surprisingly we get "fail".

We can take linear combinations of the vectors of B with the LinearCombination command. For example

LinearCombination( B, [ 1/2, 1/3, 1/4 ] );      ??

produces the same result as

1/2*B[1] + 1/3*B[2] + 1/4*B[3];      ??

We next look at ways to define a matrix.

Defining matrices

The following command defines a matrix as a list of lists.

m1:= [ [ 1, 2, 3 ], [ 2, 3, -1 ], [ 1, -2, 5] ];      ??

We can display this on the screen as a more standard looking matrix.

Display( m1 );      ??

To print the (2, 3)-entry of m1 we use

m1[2][3];      ??

To change an entry use

m1[2][2]:=300;      ??

Display( m1 );      ??

Next we define a matrix over GF(5)
We let o be the identity of GF(5) and z be the zero of GF(5)

o:= One( GF( 5 ) );      ??

z:= Zero( GF( 5 ) );      ??

m2:= [ [ o, z, z ], [ 2*o, 3*o, z ], [ z, o, 4*o ] ];      ??

Now Display gives a more readable form

Display( m2 );      ??

We can easily define matrices whose (i, j)th entry is given by a function of i and j.

m3:= List( [ 1, 2, 3 ], i -> List( [ 1, 2, 3 ], j -> i*j ) );      ??

Display(m3);      ??

Of course it is easy to define large matrices this way

m4:= List( [ 1 .. 12 ], i -> List( [ 1 .. 12 ], j -> i*j ) );      ??

Display( m4 );      ??

Defining the same matrix, but this time over Z(5) can be done by either

m5:= List( [ 1 .. 12 ], i -> List( [ 1 .. 12 ], j -> i*j*o ) );      ??

or

m5:= List( [1 .. 12], i -> List( [1 .. 12], j -> i*j ) )*o;      ??

Try Display with these matrices.

Special types of matrices

The following need little explanation
The identity matrix must be square.

id:= IdentityMat( 5 );;
Display( id );      ??

The zero matrix can be rectangular

zm:= NullMat(6, 8);;
Display( zm );      ??

zm:= NullMat(6, 7, GF(2));;
Display( zm );      ??

A diagonal matrix is square. We give it the diagonal entries.

dm:= DiagonalMat( [ 1, 1, 2, 2, 3, 4 ] );;
Display( dm );      ??

We can define a similar matrix over a finite field.

o:= One(GF( 3 ));;
dm:= DiagonalMat( [ 1, 1, 2, 2, 3, 4 ]*o );;
Display( dm );      ??

GAP will produce a random matrix whose entries are integers in the range [-10..10] with a binomial distribution.

rm:= RandomMat(10, 10);;
Display( rm );      ??

Of course, running this code again will produce a different matrix.

GAP will also produce a random matrix with rational entries -- useful for testing other routines.

Display(RandomMat(10, 10,Rationals));      ??

If GAP is asked to produce a random matrix with entries in a finite field, the entries will be distributed uniformly over the field elements.

rm1:= RandomMat(10, 10, GF(5));;
Display( rm1 );      ??

Matrix algebra

Addition, multiplication, etc. for matrices is straightforward. We give examples

m:= [1..8];;
m[1]:= [ [1,2,3], [2,3,-1], [1,-2,5] ];;
m[2]:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ];;
m[3]:= m[1] + m[2];;
m[4]:= m[1]*m[2];;
m[5]:= 3*m[1] - 7*m[2];;
m[6]:= m[1]^3;;
m[7]:= m[1]^(-1);;
m[8]:= m[6]*(5*m[5] + 6*m[2]*m[4])^(-1);;

for i in [1..8] do
  Display(m[i]);
  Print("\n");
od;      ??

To find the transpose of a matrix use TransposedMat.

m:=[ [1,2,3], [2,3,-1], [1,-2,5] ];;
mdash:=TransposedMat(m);      ??

Display(m);      ??

Display(mdash);      ??


Some determinants

We use the Determinant function in the obvious way

m1:= [ [1,2,3], [2,3,-1], [1,-2,5] ];;
m2:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ];;
Determinant( m1 );      ??

Determinant( m2 );      ??

Of course the determinant of a product is the product of the determinants

Determinant(m1*m2);      ??

If we look at matrices over GF(2) then their determinants will be either 0 or 1. What proportion would one expect for each with 2 cross 2 matrices? Let us examine the situation.

o:= One( GF(2) );;
z:= Zero( GF(2) );;

countOne:= 0;;
countZero:= 0;;
for i in [1..1000] do
  mat:= RandomMat(2, 2, GF(2));
   if Determinant(mat) = o then
     countOne:= countOne + 1;
   else
     countZero:= countZero + 1;
   fi;
od;
Print(countOne, " ", countZero, "\n");      ??

Obviously more matrices have determinant 0 than have determinant 1. Can you see why this is so? Show that the expected number with determinant 0 is 5/8 of the total, so 625 in this case.

What happens with 3 cross 3 matrices, 4 cross 4 matrices, 5 cross 5 matrices?

countOne:= 0;;
countZero:= 0;;
for i in [1..1000] do
  mat:=RandomMat(5, 5, GF(2));
    if Determinant(mat) = o then
      countOne:= countOne + 1;
    else
      countZero:= countZero + 1;
    fi;
od;
Print(countOne, " ", countZero, "\n");      ??

Before leaving this problem let us examine the case of matrices over GF(3), where of course there are three possible determinants

o:= One( GF(3) );;
z:= Zero( GF(3) );;

countTwo:= 0;;
countOne:= 0;;
countZero:= 0;;
for i in [1..1000] do
  mat:= RandomMat(2, 2, GF(3));
    if Determinant(mat) = o then
      countOne:= countOne+1;
    elif Determinant(mat) = z then
      countZero:= countZero + 1;
    else
      countTwo:= countTwo + 1;
    fi;
od;
Print(countTwo, " ",countOne, " ", countZero, "\n");      ??

Eigenvalues and eigenvectors

To compute the eigenvalues and eigenvectors of a matrix we need to specify the field

m:= [1..3];;
m[1]:= [ [1,2,3], [2,3,-1], [1,-2,5] ];;
m[2]:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ];;
m[3]:= [ [-2,-3,-3], [-1,0,-1], [5,5,6] ];;
for i in [1,2,3] do
  Print(Eigenvalues(Rationals, m[i]), "\n");
  Print(Eigenvectors(Rationals, m[i]), "\n");
  Print(CharacteristicPolynomial(m[i]), "\n");
  Print("\n");
od;      ??

Let us do the same calculation over GF(7)

o:= One( GF(7) );;
m:= [1..3];;
m[1]:= [ [1,2,3], [2,3,-1], [1,-2,5] ]*o;;
m[2]:= [ [-1,4,-3], [1,2,-1], [-1,-2,3] ]*o;;
m[3]:= [ [-2,-3,-3], [-1,0,-1], [5,5,6] ]*o;;
for i in [1, 2, 3] do
  Print(Eigenvalues(GF(7), m[i]), "\n");
  Print(Eigenvectors(GF(7), m[i]), "\n");
  Print(CharacteristicPolynomial(m[i]), "\n");
  Print("\n");
od;      ??


Systems of linear equations

If we are given a system of linear equations XA = B, where B is a row vector, then we can find a solution to the equations (of course only if one exists).

A:= [[1,2,1],[1,-1,2],[1,2,1]];      ??

B:=[1,1,4/3];      ??

SolutionMat(A,B);      ??

Now this only gives one solution, even when the system of equations has infinitely many solutions. The general solution is then SolutionMat(A, B) + v for any v in NullspaceMat(A). Let us check that fact for our particular example.

NullspaceMat(A);      ??

Since the null space has dimension 1 the general solution will be the particular solution SolutionMat(A, B) plus any multiple of the basis vector for NullspaceMat(A). For example

M:= SolutionMat(A,B) + 5/11*NullspaceMat(A)[1];      ??

M*A = B;      ??

Here is another example.

A:= [[1,2,1,1],[-1,1,1,4],[0,-1,1,2],[-1,-4,-4,-8],[-5,1,-6,-5]];      ??

B:=[-1,-1,3,8];      ??

If we compute

NullspaceMat(A);      ??

we see that this time it is two dimensional. The general solution will be the particular solution SolutionMat(A, B) plus any linear combination of the basis vector for NullspaceMat(A). For example

M:= SolutionMat(A,B) + 5/11*NullspaceMat(A)[1] - 3/4* NullspaceMat(A)[2];      ??

will satisfy MA = B.

M*A = B;      ??

In these examples the system of equations has a solution. We could have checked this by seeing that the rank of A was equal to the rank of the augmented matrix

Rank(A);      ??

To create the augmented matrix of the system we take a copy of A, then add the vector B as the final row.

augA:= ShallowCopy(A);      ??

Add(augA,B);

Now augA and A should have the same rank since our system had a solution.

Rank(augA) = Rank(A);      ??

Of course if the rank of the augmented matrix is different from the rank of A then the system has no solution. Here is an example.

B:=[1,-1,3,8];      ??

augA:= ShallowCopy(A);      ??

Add(augA,B);
Rank(augA) = Rank(A);      ??

If we tried to solve XA = B, for this different row vector B, then we will find no solution to the equations and SolutionMat(A, B) will return fail.

SolutionMat(A,B);      ??

Matrices over polynomial rings

We begin by defining R to be the ring of polynomials over the rational numbers in the indeterminate x. First we set up a list of names that GAP will use for the indeterminates; in this case it only contains one entry.

indetnames:= ["x"];      ??

Now define R to be the polynomial ring in a single indeterminate which GAP will write as we have specified.

R:= PolynomialRing(Rationals, indetnames);      ??

Now we want to use "x" for the name of the indeterminate so we set this up.

indets:= IndeterminatesOfPolynomialRing(R);      ??

x:= indets[1];      ??

We let o be the identity of R.

o:= One( R );      ??

We must be careful when we define matrices over R to make sure we use the notation 2*o for the integer 2 in R. The element 2 itself is an integer, but not an element of R. Let us define a matrix over R.

mat:= [[o, x],[2*o + x, x^2]];      ??

To check that mat is a matrix use the IsMatrix command.

IsMatrix( mat );      ??

Note that the following will not define a matrix as the elements do not come from a ring (remember that 1 and 2 are not in R but x is in R).

mat1:= [ [1, x], [2+x, x^2] ];      ??

IsMatrix(mat1);      ??

We can carry out the usual matrix operations with mat.

mat^2;      ??

mat^(-1);      ??

To convert mat1 to a matrix we could multiply it by o.

mat1:= mat1*o;      ??

IsMatrix(mat1);      ??

mat1^(-1);      ??

Determinant(mat1);      ??

Determinants of matrices over polynomial rings

We set up a polynomial ring over the rationals with 4 indeterminates. As before we first set up the names that GAP will use for the indeterminates, then set up the names that we can use for the indeterminates (which for convenience we make the same).

indetnames:= ["t","x","y","z"];;
R:= PolynomialRing(Rationals, indetnames);;
indets:= IndeterminatesOfPolynomialRing(R);      ??

t:= indets[1];;
x:= indets[2];;
y:= indets[3];;
z:= indets[4];;

Here is a 4 cross 4 matrix over R in one indeterminate

o:= One( R );;
mat:= List([1..4], i -> [2^(i-1), 5^(i-1), 17^(i-1), x^(i-1)])*o;;
Display(mat);      ??

Factors(Determinant(mat));      ??

GAP will factorise polynomials in one indeterminate but not polynomials with more than one indeterminate

o:= One(R);;
mat:= List([1..4], i -> [t^(i-1), x^(i-1), y^(i-1), z^(i-1)])*o;;
Display(mat);      ??

Determinant(mat);      ??

We define a matrix using a function

f:=function(i, j)
  if j < i then return x;
   else
     return -2;
  fi;
end;      ??

n:= 6;;
mat:= List([1 .. n], i -> List([1 .. n], j -> f(i,j)))*o;;
Display(mat);      ??

Determinant(mat);      ??

Factors(Determinant(mat));      ??

Clearly there is a theorem here!

for n in [2..12] do
  mat:= List([1 .. n], i -> List([1 .. n], j -> f(i,j)))*o;;
  Display(mat);
  Print(Determinant(mat), "\n");
  Print(Factors(Determinant(mat)), "\n");
od;      ??

JOC/EFR January 2003



Home   |   Members   |   Pre-prints   |   Library   |   Reports   |   Seminars   |   Conferences   |   Research   |   Publicity   |  Pictures   |   GAP   |   Webmail