Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
Linear Algebra in Sage Linear Algebra Tutorial Sage Days 15 May 16, 2009 Robert A. Beezer University of Puget Sound
Graphic: Gidon Eshel, Univ. of Chicago
Solving Systems of Equations Create a matrix and a vector.
A is a 3 Â 3 nonsingular matrix. 1 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
b is a 3-slot vector. A = matrix([[2,1,1/3],[-1,6,2],[1/2,1,8]]) A
0
2 B @ À1
1 3
1 6 1
1 2
1
C 2A 8
b = vector([14/3, 2, -6]) b
( 14 3 ; 2; À 6) Solve Ax = b. Solve commands ("right" is location of solution vector).
A.solve_right(b)
(2; 1; À 1)
print A.det() print A.det() == 0 299/3 False A.inverse()
0
B @
6 13 27 299 12 À 299
1 À 13
95 598 9 À 598
0
1
1 C À 23 A 3 23
Vectors are rows or columns as appropriate, compute AÀ1 b A.inverse()*b
(2; 1; À 1) Can "divide" by a matrix b/A 770 ( 299 ;
44 20 897 ; À 23 )
b/A.transpose()
2 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
(2; 1; À 1) Augment and row-reduce R = A.augment(matrix(b).transpose()) R
0
1 3
1 6 1
2 B @ À1
1 2
14 3
2 8
R.echelon_form()
0
1 @0 0
0 1 0
0 0 1
1
C 2A À6
1 2 1A À1
Properties of Matrices B is a 6 Â 5 matrix of rank 4 over the rationals Q B = matrix(QQ,[ [10,0,3,8,7], [-16,-1,-4,-10,-13], [-6,1,-3,-6,-6], [0,2,-2,-3,-2], [3,0,1,2,3], [-1,-1,1,1,0] ]) B
0
10 B B À16 B B À6 B B 0 B B @ 3 À1
0 À1 1 2 0 À1
3 À4 À3 À2 1 1
8 À10 À6 À3 2 1
B.right_kernel()
À RowSpanQ 1
À 23
1 2
B.right_kernel().basis()
3 of 14
1 7 C À13 C C À6 C C À2 C C C 3A 0 À1
À 21
Á
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
[(1; À 23 ; 21 ; À 1; À 21 )]
B.left_kernel().basis()
[(1; 0; 3; À 1; 3; 1) ; (0; 1; À 2; 1; 1; À 1)]
B.row_space()
0
1 0 0 0 B B0 1 0 0 RowSpanQ B @0 0 1 0 0 0 0 1
B.column_space()
0
1 0 0 0 B B0 1 0 0 RowSpanQ B @0 0 1 0 0 0 0 1
print print print print print print print
1 2 C À3 C C 1A À2
À 41 À 41 À 41 0
À 41
3 4 À 49
1
1
C C C A
"Rank", B.rank() "Right Nullity", B.right_nullity() "Columns", B.ncols() "Rank", B.rank() "Left Nullity", B.left_nullity() "Rows", B.nrows()
Rank 4 Right Nullity 1 Columns 5 Rank 4 Left Nullity 2 Rows 6
C is a random 50 Â 50 matrix over Q C = random_matrix(QQ, 50, num_bound=10, den_bound=10) C
4 of 14
06/13/2011 11:13 AM
B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B
Linear Algebra in BSage -- Sage
http://localhost:8000/home//21/print
À7
7 8 À 73 4 3
2 5
10 À 89
À2 3 À4 5 7
3 5 À 25
( 52 ; À 4; À
3 5
2 5 À 45 À 45
À10 2 À 31 6 1 À1 À 94
1 5
À4 À1
7 2 9 À 21 1 4 4
1 8 À 25
1
Column 5 of the matrix. C[4]
À1 À 35 À2
0
3 1 10 ; 5 ; À
À3 À 59 À 53 À 65
3 8
8 4 9
À5 À 53 À 32
9 10
8 À 25 À9
À1 À 34 5 0 À9 À2 À7 À 65 À5 À1
4 3 1 3
À4 À1 2 7 2
À7
10
9 7 8
4 À2 2 À 59 À 43 1 À 25 À8 7 À 10 À5 1 1 6 À1 5 3
1 À 51 À5 À 87 1 0 4 5 5 2
2 À2
2 5
10
À À À
4; À 51 ; 23 ; 54 ; 61 ; À 81 ; 0; 10; 0; 0; 76 ; À 1; À 2; À 25 ; À 6; À 1; À
Python slicing for entries, show() for looks show(C[20:30, 25:35])
5 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
0
B B B B B B B B B B B B B B B B B B @
8 9 7 4 9 2 3 8 1 3 1 4
http://localhost:8000/home//21/print
3 5 4 5
6
8 2
10 À 41 À 67
1 2
1 3 1 2 À 98 À 32 À 61 À 45
0 3 0 À 89
À 25 2 À 71 À 41
4 7 À 91 1 3
1 4 À 21
À 95 À 47 À 41 À9 À 25 1 9 2
3 7 9
4
0
5 2
À 10 7
1 À 21 À10
1 2 2 5
1 5 2 7 2 À 10 7
0 À 21
8 5 5 6
À 34 À 76 0 À6 À 56 À1 À1 À 51 9 10 9
7 6 À 31 À 58
1 4 7
À1
3 5
À1
1 4 4 5
À 91
1 1 5 1 2
3 2
À10
0 7 À 10 2 0 À1 À 76 1
7 6 À 51
0 0 0
2 7 À 21
1
C C C C C C C C C C C C C C C C C C A
C.det()
187558143411975607125910504689634916293080048459801797630721278449241962022353955141284605117745 668930993248348814342449533378685145998341315304441572331718197554296221857483527
C.trace()
À 15937 1260
C.norm()
40:7778652883
EigenStuff D is a 4 Â 4 matrix, not diagonalizable D = matrix(QQ, [ [-2,1,-2,-4], [12,1,4,9], [6,5,-2,-4], [3,-4,5,10] ]) show(D)
0
À2 B B 12 B @ 6 3
6 of 14
1 1 5 À4
À2 4 À2 5
1 À4 C 9C C À4 A 10
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
Build a characteristic polynomial for D (then simplify) var('t') S = D - t*identity_matrix(4) print S print p(t)=S.det() p(t) [ -t - 2 [ 12 [ 6 [ 3
1 -t + 1 5 -4
-2 -4] 4 9] -t - 2 -4] 5 -t + 10]
t4 À 7 t3 + 18 t2 À 20 t + 8
p.find_root(-10, 10)
1:9999948114 show(p.factor())
(t À 2)3 (t À 1) Easier, but less instructive: : : q(t) = D.charpoly('t') q
t 7! (((t À 7)t + 18)t À 20)t + 8
D.eigenvalues()
[1; 2; 2; 2] diag, print print print print print
evecs = D.eigenmatrix_right() "Diagonal matrix with eigenvalues" diag "Matrix with eigenvectors as columns" evecs
Diagonal matrix with eigenvalues [1 0 0 0] [0 2 0 0] [0 0 2 0] [0 0 0 2]
7 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
Matrix with eigenvectors as columns [ 1 1 0 0] [-3 -2 0 0] [-3 1 0 0] [ 0 -2 0 0] Jordan Canonical Form D.jordan_form()
0
Á 1 0 0 0 0 1C 2 1 0 C C @0 2 1AA 0 0 2
À Á 1 B0 1 B 0 B @ @0A 0
À
J, P = D.jordan_form(transformation = True) print J print print P [1|0 0 0] [-+-----] [0|2 1 0] [0|0 2 1] [0|0 0 2] [ 1 5 [ -3 -10 [ -3 5 [ 0 -10
-8 21 2 11
1] 0] 0] 1]
Check the results P^(-1)*D*P
0
1 B B0 B @0 0
0 2 0 0
1 0 0 C 1 0C C 2 1A 0 2
Manipulate the basis for the Jordan form representation Q = P.transpose().gram_schmidt()[0].transpose() Q
8 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
0
http://localhost:8000/home//21/print
75 19 À 130 19 155 19
1 B À3 B B @ À3 0
20 29 70 87 50 À 87 65 À 87
À10
15 14
1
0C C C 5 A 14 5 7
Orthogonal? Q.transpose()*Q
0
19 B 0 B B @ 0 0
4350 19
0 0
1 0 0C C C 0A
0 0
0
175 87
25 14
0
Resulting representation? Q^(-1)*D*Q
0
1 B B0 B @0 0
20 19
2 0 0
1 535 À 133 698 C À 609 C C 1A
2915 À 1653 1 2 0
2
Decompositions LU, QR, SVD, Jordan Canonical Form, Smith Normal Form, Cholesky Decomposition, : : : Convert D to a matrix E over the reals (double precision), R, to obtain QR decomposition E = D.change_ring(RDF) E
0
À2:0 B B 12:0 B @ 6:0 3:0
1:0 1:0 5:0 À4:0
À2:0 4:0 À2:0 5:0
1 À4:0 C 9:0 C C À4:0 A 10:0
ortho, triangular = E.QR() print "Orthogonal"
9 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
print ortho print print "Upper Triangular" print triangular print Orthogonal [ -0.14396315015 -0.206755085285 0.83647166171 0.486664263392] [ 0.863778900898 0.11873886424 0.366775046298 -0.324442842262] [ 0.431889450449 -0.661782341253 -0.372388950068 0.486664263392] [ 0.215944725225 0.710772502024 -0.164674510583 0.648885684523] Upper Triangular [ 13.8924439894 [ -0.0 [ -0.0 [ -0.0
2.0154841021 3.95898662912 8.78175215913] -6.24001793541 5.76589282015 11.6505245045] 0.0 -0.284437791007 -0.202100535715] -0.0 -0.0 -0.324442842262]
Checks (ortho.transpose()*ortho - identity_matrix(4)).norm()
5:57101754397 Â 10À16
(ortho*triangular - E).norm()
5:59659361322 Â 10À15
Vector Spaces Sage is so much more than numerical computation. Can work naturally with vector spaces and modules over a variety of fields and rings. F.
= FiniteField(3^2) F
F32 V is a 3-dimensional vector space over F V=F^3 V
F332 10 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
V.list()
[(0; 0; 0) ; (2a; 0; 0) ; (a + 1; 0; 0) ; (a + 2; 0; 0) ; (2; 0; 0) ; (a; 0; 0) ; (2a + 2; 0 "Generator" of all 2-D subspaces of V subs = V.subspaces(2)
Python "list comprehension" all_subs = [U for U in subs]
Grab one of the subspaces, #42
all_subs[42]
How many such subspaces? How many basis matrices in echelon_form? all_subs[86] len(all_subs)
= (32 )2 + (32 )1 + 1
11 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
Accuracy Octave and Matlab emphasize numerical results - everything is a floating point number. Their rational form is deceptive. Example by William Stein: octave:1> format rat; octave:2> a = [-86/17,40/29,-68/43,-20/11;-24/17,-1/38,-2/25,49/17] a= -86/17 40/29 -68/43 -20/11 -24/17 -1/38 -2/25 49/17 octave:3> rref(a) ans = 1 0 155/2122 -725/384 0 1 -152/173 -6553/795 and in Matlab: >> format rat; >> a = [-86/17,40/29,-68/43,-20/11;-24/17,-1/38,-2/25,49/17] a= -86/17 -24/17
40/29 -1/38
-68/43 -2/25
-20/11 49/17
>> rref(a) ans = 1 0
0 1
13/178 -152/173
-725/384 -1426/173
Now in Sage: F = matrix(2,[-86/17, 40/29, -68/43, -20/11, -24/17, -1/38, -2/25, 49/17]) show(F.echelon_form())
Entry in lower right corner:
12 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
print N(-6553/795, digits = 9), " Octave" print N(-1426/173, digits=9), " Matlab" print N(-30037214/3644069, digits = 9), " Sage"
Speed Matrices with symbolic entries var('x y') n=6 entries = [x^i-y^j+i+j for i in range(1,n+1) for j in range(1,n+1)] G = matrix(SR, n, entries) G
0
xÀy+2 B 2 Bx Ày+3 B 3 Bx Ày+4 B B x4 À y + 5 B B 5 @ x Ày+6 x6 À y + 7
G.det()
Ày 2 + x + 3 x2 À y 2 + 4 x3 À y 2 + 5 x4 À y 2 + 6 x5 À y 2 + 7 x6 À y 2 + 8
Ày 3 + x + 4 x2 À y 3 + 5 x3 À y 3 + 6 x4 À y 3 + 7 x5 À y 3 + 8 x6 À y 3 + 9
Ày 4 + x + 5 x2 À y 4 + 6 x3 À y 4 + 7 x4 À y 4 + 8 x5 À y 4 + 9 x6 À y 4 + 10
Ày 5 + x + 6 x2 À y 5 + 7 x3 À y 5 + 8 x4 À y 5 + 9 x5 À y 5 + 10 x6 À y 5 + 11
0 G.det().simplify_full()
0 time G.det()
0 Time: U 0.00 s, Wall: 0.00 s Reals, 53-bit precision H = random_matrix(RR, 10) time H.det()
À3:07894731367045
Time: U 0.12 s, Wall: 0.15 s Reals, 200-bit precision
13 of 14
06/13/2011 11:13 AM
Linear Algebra in Sage -- Sage
http://localhost:8000/home//21/print
J = random_matrix(RealField(200), 10) time J.det()
9:1432891763792026925313110748493443825766112531483542162293 Time: U 0.14 s, Wall: 0.14 s Reals, Interval Arithmetic K = (1/17.0)*random_matrix(RIF, 10, bound = 10) time K.det() Traceback (click to the left of this block for traceback) ... TypeError: random_element() got an unexpected keyword argument 'bound' Reals, Double Precision L = random_matrix(RDF, 300) timeit("L.det()") 5 loops, best of 3: 61.3 ms per loop Rationals (Exact) M = random_matrix(QQ, 800, num_bound = 10, den_bound = 10) time M.det()
À
Time: U 145.71 s, Wall: 146.92 s
14 of 14
06/13/2011 11:13 AM