Sage 6b525r
This document was ed by and they confirmed that they have the permission to share it. If you are author or own the copyright of this book, please report to us by using this report form. Report r6l17
Overview 4q3b3c
& View Sage as PDF for free.
More details 26j3b
- Words: 2,214
- Pages: 14
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
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
