Download Matrix Operations and Decompositions and more Assignments Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
Maple linalg package
Mth 351 October 27 2002 Maple 6
Bent E. Petersen
Assignment 3 - Due Nov 5 - see below
Introduction
Maple 6 comes with two packages for linear algebra, linalg and LinearAlgebra. In this worksheet we
will only explore the linalg package. Many of the functions in the linalg package take optional
parameters, or have several forms. We will not investigate every, or even many, possibilities.
Note that Maple commands are polymorphic. Their behavior depends on the number and type of
parameters that you feed to them. In some situations you may need to be careful about data types and
the number of parameters passed.
The Assignment
Assignment 3 is included near the end of this worksheet.
> restart;
> with(linalg):
Warning, the protected names norm and trace have been redefined and unprotected
matrix()
Matrices can be entered by listing their rows (lists are denoted by [ ] in Maple) or by giving the size and
then listing the entries in the matrix row-wise. Here I enter the matrix A in both ways.
> A:=matrix([[1,2,3,4],[5,6,7,8],[9,0,1,2],[4,4,4,4]]);
A :=
> A:=matrix(4,4,[1,2,3,4,5,6,7,8,9,0,1,2,4,4,4,4]);
A :=
Personally I prefer the second method for interactive entry. It is also possible to build a matrix by
specifying its columns (as vectors or lists).
Multiplication of matrices, evalm()
Matrix multiplication is denoted by &*. The ampersand serves to mark the multiplication as being
noncommutative so that Maple will not make incorrect assumptions when simplifying expressions
involving matrix multiplication.
> B:=matrix(4,2,[3,2,-4,-5,7,6,-3,-5]);
B :=
> A&*B;
A &* B
Note that Maple often suppresses output of matrices. To force evaluation we use the evalm() function
> evalm(A&B);*
transpose()
> B:=matrix(2,3,[1,2,3,4,5,6]);
B :=
> transpose(B);
Matrices with unspecified entries
We can construct matrices with unspecified entries and then specify the entries later:
> C:=matrix(3,4);
> band([1,3,1],6);
extend()
The extend function enlarges a matrix. The function call extend(B,m,n,x) adds m rows and n columns
to the matrix B. The new entries are initialized to x or left undefined if the value x is not included in the
function call.
> B:=matrix(2,3,[1,2,3,4,5,6]);
B :=
> C0:=extend(B,2,1,0);
C0 :=
> C:=extend(B,2,1);
C :=
1 2 3 C 1 4,
4 5 6 C 2 4,
C
3 1,
C
3 2,
C
3 3,
C
3 4,
C 4 1,
C
4 2,
C
4 3,
C
4 4,
The identity matrix
The expression &*() may be used to designate the identity matrix of unspecified size. The size is
determined from the context. For example
> evalm( &() + matrix(2,2,[1,2,3,4]) );*
To produce an identity matrix of an explicit size we can use the diag() command in conjunction with
the seq() command.
> diag( seq(1,k=1..5) );
The function band() is even more convenient for creating the identity matrix
> band([1],5);
Note the identity matrix is implicitly assumed when scalars are added to matrices:
> evalm( 3 + matrix(2,2,[1,2,3,4]) );
Hilbert matrix, hilbert()
The (i,j)th entry in the Hilbert matrix is 1/(i+j-1).
> hilbert(4);
Zero matrix
A zero matrix of any size may be easily created by using the sequence command, seq():
> matrix(3,5,[seq(0,k=1..35)]);*
> C:=matrix(3,3,[1,2,2,6,5,2,7,8,9]);
C :=
> det(C);
vector()
A vector in Maple is an ordered n-tuple. There is no concept of rows or columns iin connection with a
vector; We can have row or column vectors if we wish, but they are matrices and so have two
subscripts.
> a:=vector([1,2,3,4]);
a :=[ 1 2 3 4, , , ]
> b:=matrix(1,4,[1,2,3,4]);
b :=[ 1 2 3 4 ]
> c:=matrix(4,1,[1,2,3,4]);
c :=
We can convert vectors to matrices and 1-by-n matrices to vectors
> convert(a,matrix);
> convert(c,vector);
[ 1 2 3 4, , , ]
rowspace(), colspace()
The rowspace() function returns a basis for the rowspace of a matrix. The basis is returned as a set of
vectors, that is, unordered.
> basis:=rowspace(A);
basis :={ [ 1 0 0 0, , , ], [ 0 1 0 -1, , , ],[ 0 0 1 2, , , ]}
The quickest way to obtained an ordered basis is to convert to a list. However, the order may not be
what you want.
> convert(basis,list);
[ [ 1 0 0 0, , , ], [ 0 1 0 -1, , , ],[ 0 0 1 2, , , ]]
The colspace() function behaves similarly
> colspace(A);
{ [ 1 0 0 -1, , , ], [ 0 1 0 1, , , ],[ 0 0 1 0, , , ]}
Note both rowspace() and colspace() return vectors.
Augmentation - augment(), concat() and stackmatrix()
We can augment an m-by-n matix by an m-by-q matrix or by a vector or list. The matrices are simply
concatenated horizontally. The vector or list is added as a column. The function concat() is a synonym
for augment(). Recall our matrix A
> evalm(A);
> augment([1,2,3,4],A,[5,6,7,8]);
> augment(A,A^2);
The stackmatrix() command concatenates vertically. A vector or a list is added as a row.
> stackmatrix([0,0,1,2],A,[1,2,0,0]);
> evalm(A);
> rref(A);
Gauss elimination - gausselim()
The function gausselim() does Gauss elimination, or row reduction, but stops when a triangular matrix
is achieved. It does not go all the way to the fully row reduced row echelon form.
> B:=matrix(3,4,[x,x,3,4,x,y,z,8,y+1,y,3,4]);
B :=
x x 3 4
x y z 8
y + 1 y 3 4
> gausselim(B);
x x 3 4
x − y − 1
x
x − y − 1
x
− z x − 6 y x + 3 y + +
2 3 y 3 x
2
x
− 2 x − 2 y x + y + +
2 y x
2
x
Fraction free row reduction - ffgausselim()
The function ffgausselim does Gauss elimination, or row reduction, without introducing fractions
> ffgausselim(B);
x x 3 4
0 − x 3 x − 3 y − 3 4 x − 4 y − 4
0 0 − z x − 6 y x + 3 y + +
2 3 y 3 x
2 − 8 x − 8 y x + 4 y + +
2 4 y 4 x
2
Hermite normal form for integer matrices - ihermite()
The function call ihermite(B) row reduces the matrix B over the integers. The Hermite normal form is
obtained by doing only integer elementary row operations. This includes interchanging rows,
multiplying a row by -1, and adding an integral multiple of one row to another. The goal is to produce
an upper triangular matrix with as many zeros in each pivotal column as possible.
> C1:=matrix(3,4,[13,11,3,0,7,3,5,0,11,2,1,9]);
C1 :=
> ihermite(C1);
The function call ihermite(B,U) in addition returns an integer matrix U such that U &* B = ihermite(B).
> C2:=matrix(3,3,[1,11,3,5,3,5,1,1,1]);
C2 :=
> ihermite(C2,U);
> evalm(U);
> evalm( U & C2 );*
> C3:=matrix(3,3,[2,2,1,2,2,2,0,0,2]);
C3 :=
> ihermite(C3);
LU decomposition - LUdecomp()
The function LUdecomp() returns only the U in the LU decomposition. To get the permutation matrix
> evalm(Lh&Uh);*
We see it checks, B = LU and P is the identity.
Cholesky decomposition - cholesky(), definite()
Here is an example of a positive definite matrix
> B:=hilbert(4);
B :=
> definite(B,positive_def);
true
In the Cholesky decomposition we have a lower traingualr matrix C such the B = C &* transpose(C).
> C:=cholesky(B);
C :=
Let’s very that it works -
> evalm( C & transpose(C) );*
Sure enough!
We can produce a positive definite matrix from any nonsingular matrix by multiplying by the transpose
(check this):
> B0:=matrix(3,3,[1,2,2,3,3,3,2,2,1]); det(B0);
B0 :=
> B:= evalm(transpose(B0) & B0);*
B :=
> C:=cholesky(B);
C :=
A :=
> b:=vector([2,3,1]);
b :=[ 2 3 1, , ]
> soln:=linsolve(A,b);
soln :=
− 1 − _t , , , 1 _t 1
Maple returns a list of values. The underscore t sub 1 is a parameter. You assign it arbitrary values to
obtain all the solutions. If you find the expression difficult to read you can always substitute something
more agreeable, for example
> soln:=linsolve(A,b,rnk,s);
soln :=
− 1 − s 1 , s 1 , ,
Note the parameter name has to appear as the fourth variable, so we had to insert a third variable, rnk.
It turns out that linsolve stuffs the rank of A into the third variable. Thus
> rnk;
One can pick out the individual components easily
> soln[1]; soln[2]; soln[3]; soln[4];
− 1 − s 1
s 1
We can also solve directly by row reduction of course -
> M:=augment(A,b); R:=rref(M);
M :=
R :=
Here we can read the solution easily. Moreover, we get the canonical form of the solution, with the free
variable(s) as the parameter. This is not always the form that Maple returns (though in the present
example it is).
> x[1]:=R[1,5]-R[1,2]s; x[2]:=s; x[3]:=R[2,5]; x[4]:=R[3,5];*
x := 1
− 1 − s
x 2 := s
x := 3
x 4 :=-
Least squares solution - leastsqrs()
The function call leastsqrs(A,b) returns the solutions which minimize the Euclidean norm of Ax-b.
First here is an example that has an actual solution:
> A:=matrix(3,4,[1,1,2,2,3,3,4,4,5,5,6,7]);
A :=
> b:=vector([2,3,1]);
b :=[ 2 3 1, , ]
> linsolve(A,b);
− 1 − _t 1 , _t 1 , ,
> leastsqrs(A,b);
− 1 − _t 1 , _t 1 , ,
Now here is an example which has no actual solution. Note Maple indicates the lack of a solution by
returning the empty set.
> c:=vector([1,2,3,4]);
c :=[ 1 2 3 4, , , ]
M4 :=
> M5:=mulrow(M4,2,-1/4);
M5 :=
> M6:=addrow(M5,2,3,-1);
M6 :=
> M7:=mulrow(M6,3,4);
M7 :=
> M8:=addrow(M7,3,1,-2);
M8 :=
> M9:=addrow(M8,3,2,-3/4);
M9 :=
> M10:=addrow(M9,2,1,-2);
M10 :=
Let’s check our work
> rref(M);
More row reduction - pivot()
The function call pivot(A, i, j) will add multiples of the ith row to every other row in the matrix, with
the result that the jth column of the matrix will be all zeros, except for the (i,j)th element. The function
call pivot(A, i, j, r..s) acts like pivot(A, i, j) except that only rows r through s are set to zero in the jth
column (except row i if it is in the range).
> A:=matrix(3,4,[1,2,2,4,2,3,3,5,6,4,1,2]);
A :=
> pivot(A,2,3);
> pivot(A,2,3,1..2);
GramSchmidt()
The function GramSchmidt() returns a set of orthogonal vectors with the same span as a given list or
set of linearly independent vectors. The vectors are normalized if that option is specified.
> v1:=vector([1,2,3]);
v1 :=[ 1 2 3, , ]
> v2:=vector([1,2,0]);
v2 :=[ 1 2 0, , ]
> v3:=vector([1,0,0]);
v3 :=[ 1 0 0, , ]
> GramSchmidt([v1,v2,v3],normalized);
