Linear Algebra in Sage
system:sage

<p style="text-align: center;"><span style="font-size: xx-large;">Linear Algebra Tutorial<br />Sage Days 15<br />May 16, 2009<br /><br />Robert A. Beezer</span></p>
<p style="text-align: center;"><span style="font-size: xx-large;">University of Puget Sound<br /></span></p>
<p style="text-align: center;"><img src="linear-algebra.gif" alt="" width="557" height="359" /></p>
<p style="text-align: center;">Graphic: Gidon Eshel, Univ. of Chicago</p>
<p><span style="font-size: x-large;">Solving Systems of Equations</span></p>
<p>Create a matrix and a vector.</p>
<p>$A$ is a $3\times 3$ nonsingular matrix.</p>
<p>$b$ is a 3-slot vector.</p>

{{{id=0|
A = matrix([[2,1,1/3],[-1,6,2],[1/2,1,8]])
A
///
}}}

{{{id=3|
b = vector([14/3, 2, -6])
b
///
}}}

<p>Solve $Ax=b$.</p>
<p>Solve commands ("right" is location of solution vector).</p>

{{{id=4|
A.solve_right(b)
///
}}}

{{{id=7|
print A.det()
print A.det() == 0
///
}}}

{{{id=5|
A.inverse()
///
}}}

<p>Vectors are rows or columns as appropriate, compute $A^{-1}b$</p>

{{{id=8|
A.inverse()*b
///
}}}

<p>Can "divide" by a matrix</p>

{{{id=9|
b/A
///
}}}

{{{id=12|
b/A.transpose()
///
}}}

<p>Augment and row-reduce</p>

{{{id=16|
R = A.augment(matrix(b).transpose())
R
///
}}}

{{{id=13|
R.echelon_form()
///
}}}

<p><span style="font-size: x-large;">Properties of Matrices</span></p>

<p>$B$ is a $6\times 5$ matrix of rank 4 over the rationals $\QQ$</p>

{{{id=19|
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
///
}}}

{{{id=22|
B.right_kernel()
///
}}}

{{{id=23|
B.right_kernel().basis()
///
}}}

{{{id=25|
B.left_kernel().basis()
///
}}}

{{{id=26|
B.row_space()
///
}}}

{{{id=27|
B.column_space()
///
}}}

{{{id=28|
print "Rank", B.rank()
print "Right Nullity", B.right_nullity()
print "Columns", B.ncols()
print
print "Rank", B.rank()
print "Left Nullity", B.left_nullity()
print "Rows", B.nrows()
///
}}}

<p>$C$ is a random $50\times 50$ matrix over $\QQ$</p>

{{{id=35|
C = random_matrix(QQ, 50, num_bound=10, den_bound=10)
C
///
}}}

<p>Column 5 of the matrix.</p>

{{{id=39|
C[4]
///
}}}

<p>&nbsp; Python slicing for entries, show() for looks</p>

{{{id=29|
show(C[20:30, 25:35])
///
}}}

{{{id=32|
C.det()
///
}}}

{{{id=33|
C.trace()
///
}}}

{{{id=44|
C.norm()
///
}}}

<p><span style="font-size: x-large;">EigenStuff</span></p>

<p>$D$ is a $4\times 4$ matrix, not diagonalizable</p>

{{{id=45|
D = matrix(QQ, [
[-2,1,-2,-4],
[12,1,4,9],
[6,5,-2,-4],
[3,-4,5,10]
])
show(D)
///
}}}

<p>Build a characteristic polynomial for $D$ (then simplify)</p>

{{{id=48|
var('t')
S = D - t*identity_matrix(4)
print S
print
p(t)=S.det()
p(t)
///
}}}

{{{id=49|
p.find_root(-10, 10)
///
}}}

{{{id=52|
show(p.factor())
///
}}}

<p>Easier, but less instructive$\dots$</p>

{{{id=53|
q(t) = D.charpoly('t')
q
///
}}}

{{{id=55|
D.eigenvalues()
///
}}}

{{{id=56|
diag, evecs = D.eigenmatrix_right()
print "Diagonal matrix with eigenvalues"
print diag
print
print "Matrix with eigenvectors as columns"
print evecs
///
}}}

<p>Jordan Canonical Form</p>

{{{id=57|
D.jordan_form()
///
}}}

{{{id=60|
J, P = D.jordan_form(transformation = True)
print J
print
print P
///
}}}

<p>Check the results</p>

{{{id=61|
P^(-1)*D*P
///
}}}

<p>Manipulate the basis for the Jordan form representation</p>

{{{id=66|
Q = P.transpose().gram_schmidt()[0].transpose()
Q
///
}}}

<p>Orthogonal?</p>

{{{id=67|
Q.transpose()*Q
///
}}}

<p>Resulting representation?</p>

{{{id=68|
Q^(-1)*D*Q
///
}}}

<p><span style="font-size: x-large;">Decompositions</span></p>

<p>LU, QR, SVD, Jordan Canonical Form, Smith Normal Form, Cholesky Decomposition, $\dots$</p>

<p>Convert $D$ to a matrix $E$ over the reals (double precision), $\RR$, to obtain QR decomposition</p>

{{{id=69|
E = D.change_ring(RDF)
E
///
}}}

{{{id=81|
ortho, triangular = E.QR()
print "Orthogonal"
print ortho
print
print "Upper Triangular"
print triangular
print
///
}}}

<p>Checks</p>

{{{id=85|
(ortho.transpose()*ortho - identity_matrix(4)).norm()
///
}}}

{{{id=82|
(ortho*triangular - E).norm()
///
}}}

<p><span style="font-size: x-large;">Vector Spaces</span></p>
<p>Sage is so much more than numerical computation.</p>
<p>Can work naturally with vector spaces and modules over a variety of fields and rings.</p>

{{{id=84|
F.<a> = FiniteField(3^2)
///
}}}

{{{id=89|
F
///
}}}

<p>$V$ is a 3-dimensional vector space over $F$</p>

{{{id=90|
V=F^3
V
///
}}}

{{{id=91|
V.list()
///
}}}

<p>"Generator" of all 2-D subspaces of $V$</p>

{{{id=93|
subs = V.subspaces(2)
///
}}}

<p>Python "list comprehension"</p>

{{{id=95|
all_subs = [U for U in subs]
///
}}}

<p>Grab one of the subspaces, #42</p>

{{{id=101|
all_subs[42]
///
}}}

<p>How many such subspaces?</p>
<p>How many basis matrices in echelon_form?</p>

{{{id=104|
all_subs[86]
///
}}}

{{{id=97|
len(all_subs)
///
}}}

<p>&nbsp;</p>
<p>&nbsp;</p>
<p>&nbsp;</p>
<p>&nbsp;</p>
<p>&nbsp;</p>
<p>$=(3^2)^2 + (3^2)^1 + 1$</p>
<p>&nbsp;</p>
<p>&nbsp;</p>
<p><span style="font-size: x-large;">Accuracy</span></p>
<p>Octave and Matlab emphasize numerical results - everything is a floating point number.</p>
<p>Their rational form is deceptive.&nbsp; Example by William Stein:</p>
<p>octave:1&gt; format <strong style="color: black; background-color: #a0ffff;">rat</strong>; <br /> octave:2&gt; a = [-86/17,40/29,-68/43,-20/11;-24/17,-1/38,-2/25,49/17] <br /> a = <br /> &nbsp; &nbsp; &nbsp;-86/17 &nbsp; &nbsp; &nbsp;40/29 &nbsp; &nbsp; -68/43 &nbsp; &nbsp; -20/11 <br /> &nbsp; &nbsp; &nbsp;-24/17 &nbsp; &nbsp; &nbsp;-1/38 &nbsp; &nbsp; &nbsp;-2/25 &nbsp; &nbsp; &nbsp;49/17 <br /> octave:3&gt; rref(a) <br /> ans = <br /> &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 1 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0 &nbsp; 155/2122 &nbsp; -725/384 <br /> &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 0 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;1 &nbsp; -152/173 &nbsp;-6553/795</p>
<p>and in<span style="background-color: #ffffff;"> Matlab:</span></p>
<div id="qhide_150749" class="qt" style="display: block;">&gt;&gt; format <strong style="color: black; background-color: #a0ffff;">rat</strong>; <br /> &gt;&gt; a = [-86/17,40/29,-68/43,-20/11;-24/17,-1/38,-2/25,49/17] <br /> <br /></div>
<p>a = <br /> &nbsp; &nbsp; &nbsp;-86/17 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;40/29 &nbsp; &nbsp; &nbsp; &nbsp; -68/43 &nbsp; &nbsp; &nbsp; &nbsp; -20/11 <br /> &nbsp; &nbsp; &nbsp;-24/17 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;-1/38 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;-2/25 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;49/17</p>
<div id="qhide_150750" class="qt" style="display: block;">&gt;&gt; rref(a) <br /> <br /></div>
<p>ans = <br /> &nbsp; &nbsp; &nbsp; &nbsp;1 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;0 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 13/178 &nbsp; &nbsp; &nbsp; -725/384 <br /> &nbsp; &nbsp; &nbsp; &nbsp;0 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;1 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; -152/173 &nbsp; &nbsp; &nbsp;-1426/173</p>
<p>Now in Sage:</p>

{{{id=107|
F = matrix(2,[-86/17, 40/29, -68/43, -20/11, -24/17, -1/38, -2/25, 49/17]) 
show(F.echelon_form())
///
}}}

<p>Entry in lower right corner:</p>

{{{id=130|
print N(-6553/795, digits = 9), "  Octave"
print N(-1426/173, digits=9), "  Matlab"
print N(-30037214/3644069, digits = 9), "  Sage"
///
}}}

<p><span style="font-size: x-large;">Speed</span></p>
<p>Matrices with symbolic entries</p>

{{{id=54|
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
///
}}}

{{{id=114|
G.det()
///
}}}

{{{id=115|
G.det().simplify_full()
///
}}}

{{{id=116|
time G.det()
///
}}}

<p>Reals, 53-bit precision</p>

{{{id=117|
H = random_matrix(RR, 10)
time H.det()
///
}}}

<p>Reals, 200-bit precision</p>

{{{id=119|
J = random_matrix(RealField(200), 10)
time J.det()
///
}}}

<p>Reals, Interval Arithmetic</p>

{{{id=123|
K = (1/17.0)*random_matrix(RIF, 10, bound = 10)
time K.det()
///
}}}

<p>Reals, Double Precision</p>

{{{id=126|
L = random_matrix(RDF, 300)
timeit("L.det()")
///
}}}

<p>Rationals (Exact)</p>

{{{id=128|
M = random_matrix(QQ, 800, num_bound = 10, den_bound = 10)
time M.det()
///
}}}

{{{id=129|

///
}}}