Algebra Lineal en Sage
system:sage


<p style="text-align: center;"><span style="font-size: xx-large;">Tutorial de &Aacute;lgebra Lineal<br />Sage Days 15<br />Mayo 16, 2009<br /><br />Robert A. Beezer</span></p>
<p style="text-align: center;"><span style="font-size: xx-large;">Universidad de 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;">Gr&aacute;fico: Gidon Eshel, Univ. de Chicago</p>
<p><span style="font-size: x-large;">Resolviendo Sistema de Ecuaciones</span></p>
<p>Crear una matriz y un vector.</p>
<p>$A$ es una matriz no singular $3\times 3$</p>
<p>$b$ es un 3-slot vector.</p>

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

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

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

{{{id=4|
A.solve_right(b)
///
(2, 1, -1)
}}}

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

{{{id=5|
A.inverse()
///
[   6/13   -1/13       0]
[ 27/299  95/598   -1/23]
[-12/299  -9/598    3/23]
}}}

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

{{{id=8|
A.inverse()*b
///
(2, 1, -1)
}}}

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

{{{id=9|
b/A
///
(770/299, 44/897, -20/23)
}}}

{{{id=12|
b/A.transpose()
///
(2, 1, -1)
}}}

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

{{{id=16|
R = A.augment(matrix(b).transpose())
R
///
[   2    1  1/3 14/3]
[  -1    6    2    2]
[ 1/2    1    8   -6]
}}}

{{{id=13|
R.echelon_form()
///
[ 1  0  0  2]
[ 0  1  0  1]
[ 0  0  1 -1]
}}}

<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
///
[ 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]
}}}

{{{id=22|
B.right_kernel()
///
Vector space of degree 5 and dimension 1 over Rational Field
Basis matrix:
[   1 -3/2  1/2   -1 -1/2]
}}}

{{{id=23|
B.right_kernel().basis()
///
[
(1, -3/2, 1/2, -1, -1/2)
]
}}}

{{{id=25|
B.left_kernel().basis()
///
[
(1, 0, 3, -1, 3, 1),
(0, 1, -2, 1, 1, -1)
]
}}}

{{{id=26|
B.row_space()
///
Vector space of degree 5 and dimension 4 over Rational Field
Basis matrix:
[ 1  0  0  0  2]
[ 0  1  0  0 -3]
[ 0  0  1  0  1]
[ 0  0  0  1 -2]
}}}

{{{id=27|
B.column_space()
///
Vector space of degree 6 and dimension 4 over Rational Field
Basis matrix:
[   1    0    0    0 -1/4 -1/4]
[   0    1    0    0 -1/4  3/4]
[   0    0    1    0 -1/4 -9/4]
[   0    0    0    1    0    1]
}}}

{{{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()
///
Rank 4
Right Nullity 1
Columns 5

Rank 4
Left Nullity 2
Rows 6
}}}

<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
///
50 x 50 dense matrix over Rational Field (type 'print C.str()' to see all of the entries)
}}}

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

{{{id=39|
C[4]
///
(3/2, 2/5, 0, 8/5, -2/5, -4/7, 7/4, -3/5, -2/3, 1/9, 7, 2, 0, 5/4, -7/2, -10/9, 5/7, 2, 0, 1, -9/2, 1/5, 7/2, 1, -10/3, -1, -9, -10, -5/3, 0, -1/3, -9/10, -2/5, -7/6, 5/8, 0, 3/4, -2, -3/2, -5/3, -1/3, -2, -1, 1, 1/10, 1/5, 5/4, 9/10, -8, 9/7)
}}}

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

{{{id=29|
show(C[20:30, 25:35])
///
<html><div class="math">\left(\begin{array}{rrrrrrrrrr}
-\frac{7}{6} & 1 & \frac{8}{5} & \frac{5}{3} & -\frac{4}{7} & \frac{5}{3} & -\frac{9}{7} & -\frac{7}{4} & -5 & -\frac{5}{4} \\
-\frac{7}{8} & -\frac{4}{3} & -7 & -\frac{9}{4} & \frac{3}{5} & -\frac{2}{7} & \frac{10}{3} & -1 & 6 & -\frac{7}{9} \\
0 & -2 & 1 & -\frac{5}{4} & \frac{1}{10} & 0 & 5 & 1 & -\frac{4}{3} & -\frac{7}{6} \\
\frac{5}{4} & 0 & \frac{4}{9} & \frac{5}{4} & 10 & -\frac{3}{8} & \frac{3}{4} & \frac{3}{2} & \frac{5}{9} & -2 \\
\frac{5}{4} & -7 & -\frac{2}{9} & 1 & -\frac{8}{9} & \frac{10}{3} & -5 & -\frac{7}{10} & 1 & -\frac{10}{3} \\
3 & \frac{2}{5} & \frac{5}{6} & 0 & \frac{4}{9} & -\frac{1}{4} & -\frac{3}{5} & -\frac{1}{3} & 7 & \frac{10}{3} \\
-\frac{7}{9} & 5 & 2 & \frac{6}{7} & -\frac{5}{8} & -\frac{5}{4} & \frac{3}{2} & \frac{1}{2} & \frac{10}{7} & \frac{8}{9} \\
1 & \frac{3}{4} & \frac{3}{10} & 4 & 0 & -\frac{1}{2} & \frac{3}{4} & 5 & 5 & 0 \\
1 & 0 & \frac{1}{2} & \frac{7}{10} & 5 & -\frac{3}{2} & \frac{3}{2} & 0 & 0 & -\frac{5}{2} \\
-\frac{7}{5} & \frac{3}{4} & -\frac{7}{4} & 0 & -\frac{2}{3} & -\frac{9}{4} & 1 & \frac{10}{9} & -\frac{3}{5} & -\frac{5}{7}
\end{array}\right)</div></html>
}}}

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

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

{{{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)
///
<html><div class="math">\left(\begin{array}{rrrr}
-2 & 1 & -2 & -4 \\
12 & 1 & 4 & 9 \\
6 & 5 & -2 & -4 \\
3 & -4 & 5 & 10
\end{array}\right)</div></html>
}}}

<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)
///
[ -t - 2       1      -2      -4]
[     12  -t + 1       4       9]
[      6       5  -t - 2      -4]
[      3      -4       5 -t + 10]

t^4 - 7*t^3 + 18*t^2 - 20*t + 8
}}}

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

{{{id=52|
show(p.factor())
///
<html><div class="math">{\left(t - 2\right)}^{3} {\left(t - 1\right)}</div></html>
}}}

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

{{{id=53|
q(t) = D.charpoly('t')
q
///
t |--> (((t - 7)*t + 18)*t - 20)*t + 8
}}}

{{{id=55|
D.eigenvalues()
///
[1, 2, 2, 2]
}}}

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

Matrix with eigenvectors as columns
[ 1  1  0  0]
[-3 -2  0  0]
[-3  1  0  0]
[ 0 -2  0  0]
}}}

<p>Jordan Canonical Form</p>

{{{id=57|
D.jordan_form()
///
[1|0 0 0]
[-+-----]
[0|2 1 0]
[0|0 2 1]
[0|0 0 2]
}}}

{{{id=60|
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  -8   1]
[ -3 -10  21   0]
[ -3   5   2   0]
[  0 -10  11   1]
}}}

<p>Check the results</p>

{{{id=61|
P^(-1)*D*P
///
[1 0 0 0]
[0 2 1 0]
[0 0 2 1]
[0 0 0 2]
}}}

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

{{{id=66|
Q = P.transpose().gram_schmidt()[0].transpose()
Q
///
[      1   75/19   20/29   15/14]
[     -3 -130/19   70/87       0]
[     -3  155/19  -50/87    5/14]
[      0     -10  -65/87     5/7]
}}}

<p>Orthogonal?</p>

{{{id=67|
Q.transpose()*Q
///
[     19       0       0       0]
[      0 4350/19       0       0]
[      0       0  175/87       0]
[      0       0       0   25/14]
}}}

<p>Resulting representation?</p>

{{{id=68|
Q^(-1)*D*Q
///
[         1      20/19 -2915/1653   -535/133]
[         0          2          1   -698/609]
[         0          0          2          1]
[         0          0          0          2]
}}}

<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
///
[-2.0  1.0 -2.0 -4.0]
[12.0  1.0  4.0  9.0]
[ 6.0  5.0 -2.0 -4.0]
[ 3.0 -4.0  5.0 10.0]
}}}

{{{id=81|
ortho, triangular = E.QR()
print "Orthogonal"
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    2.0154841021   3.95898662912   8.78175215913]
[           -0.0  -6.24001793541   5.76589282015   11.6505245045]
[           -0.0             0.0 -0.284437791007 -0.202100535715]
[           -0.0            -0.0            -0.0 -0.324442842262]
}}}

<p>Checks</p>

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

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

<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
///
Finite Field in a of size 3^2
}}}

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

{{{id=90|
V=F^3
V
///
Vector space of dimension 3 over Finite Field in a of size 3^2
}}}

{{{id=91|
V.list()
///
[(0, 0, 0), (2*a, 0, 0), (a + 1, 0, 0), (a + 2, 0, 0), (2, 0, 0), (a, 0, 0), (2*a + 2, 0, 0), (2*a + 1, 0, 0), (1, 0, 0), (0, 2*a, 0), (2*a, 2*a, 0), (a + 1, 2*a, 0), (a + 2, 2*a, 0), (2, 2*a, 0), (a, 2*a, 0), (2*a + 2, 2*a, 0), (2*a + 1, 2*a, 0), (1, 2*a, 0), (0, a + 1, 0), (2*a, a + 1, 0), (a + 1, a + 1, 0), (a + 2, a + 1, 0), (2, a + 1, 0), (a, a + 1, 0), (2*a + 2, a + 1, 0), (2*a + 1, a + 1, 0), (1, a + 1, 0), (0, a + 2, 0), (2*a, a + 2, 0), (a + 1, a + 2, 0), (a + 2, a + 2, 0), (2, a + 2, 0), (a, a + 2, 0), (2*a + 2, a + 2, 0), (2*a + 1, a + 2, 0), (1, a + 2, 0), (0, 2, 0), (2*a, 2, 0), (a + 1, 2, 0), (a + 2, 2, 0), (2, 2, 0), (a, 2, 0), (2*a + 2, 2, 0), (2*a + 1, 2, 0), (1, 2, 0), (0, a, 0), (2*a, a, 0), (a + 1, a, 0), (a + 2, a, 0), (2, a, 0), (a, a, 0), (2*a + 2, a, 0), (2*a + 1, a, 0), (1, a, 0), (0, 2*a + 2, 0), (2*a, 2*a + 2, 0), (a + 1, 2*a + 2, 0), (a + 2, 2*a + 2, 0), (2, 2*a + 2, 0), (a, 2*a + 2, 0), (2*a + 2, 2*a + 2, 0), (2*a + 1, 2*a + 2, 0), (1, 2*a + 2, 0), (0, 2*a + 1, 0), (2*a, 2*a + 1, 0), (a + 1, 2*a + 1, 0), (a + 2, 2*a + 1, 0), (2, 2*a + 1, 0), (a, 2*a + 1, 0), (2*a + 2, 2*a + 1, 0), (2*a + 1, 2*a + 1, 0), (1, 2*a + 1, 0), (0, 1, 0), (2*a, 1, 0), (a + 1, 1, 0), (a + 2, 1, 0), (2, 1, 0), (a, 1, 0), (2*a + 2, 1, 0), (2*a + 1, 1, 0), (1, 1, 0), (0, 0, 2*a), (2*a, 0, 2*a), (a + 1, 0, 2*a), (a + 2, 0, 2*a), (2, 0, 2*a), (a, 0, 2*a), (2*a + 2, 0, 2*a), (2*a + 1, 0, 2*a), (1, 0, 2*a), (0, 2*a, 2*a), (2*a, 2*a, 2*a), (a + 1, 2*a, 2*a), (a + 2, 2*a, 2*a), (2, 2*a, 2*a), (a, 2*a, 2*a), (2*a + 2, 2*a, 2*a), (2*a + 1, 2*a, 2*a), (1, 2*a, 2*a), (0, a + 1, 2*a), (2*a, a + 1, 2*a), (a + 1, a + 1, 2*a), (a + 2, a + 1, 2*a), (2, a + 1, 2*a), (a, a + 1, 2*a), (2*a + 2, a + 1, 2*a), (2*a + 1, a + 1, 2*a), (1, a + 1, 2*a), (0, a + 2, 2*a), (2*a, a + 2, 2*a), (a + 1, a + 2, 2*a), (a + 2, a + 2, 2*a), (2, a + 2, 2*a), (a, a + 2, 2*a), (2*a + 2, a + 2, 2*a), (2*a + 1, a + 2, 2*a), (1, a + 2, 2*a), (0, 2, 2*a), (2*a, 2, 2*a), (a + 1, 2, 2*a), (a + 2, 2, 2*a), (2, 2, 2*a), (a, 2, 2*a), (2*a + 2, 2, 2*a), (2*a + 1, 2, 2*a), (1, 2, 2*a), (0, a, 2*a), (2*a, a, 2*a), (a + 1, a, 2*a), (a + 2, a, 2*a), (2, a, 2*a), (a, a, 2*a), (2*a + 2, a, 2*a), (2*a + 1, a, 2*a), (1, a, 2*a), (0, 2*a + 2, 2*a), (2*a, 2*a + 2, 2*a), (a + 1, 2*a + 2, 2*a), (a + 2, 2*a + 2, 2*a), (2, 2*a + 2, 2*a), (a, 2*a + 2, 2*a), (2*a + 2, 2*a + 2, 2*a), (2*a + 1, 2*a + 2, 2*a), (1, 2*a + 2, 2*a), (0, 2*a + 1, 2*a), (2*a, 2*a + 1, 2*a), (a + 1, 2*a + 1, 2*a), (a + 2, 2*a + 1, 2*a), (2, 2*a + 1, 2*a), (a, 2*a + 1, 2*a), (2*a + 2, 2*a + 1, 2*a), (2*a + 1, 2*a + 1, 2*a), (1, 2*a + 1, 2*a), (0, 1, 2*a), (2*a, 1, 2*a), (a + 1, 1, 2*a), (a + 2, 1, 2*a), (2, 1, 2*a), (a, 1, 2*a), (2*a + 2, 1, 2*a), (2*a + 1, 1, 2*a), (1, 1, 2*a), (0, 0, a + 1), (2*a, 0, a + 1), (a + 1, 0, a + 1), (a + 2, 0, a + 1), (2, 0, a + 1), (a, 0, a + 1), (2*a + 2, 0, a + 1), (2*a + 1, 0, a + 1), (1, 0, a + 1), (0, 2*a, a + 1), (2*a, 2*a, a + 1), (a + 1, 2*a, a + 1), (a + 2, 2*a, a + 1), (2, 2*a, a + 1), (a, 2*a, a + 1), (2*a + 2, 2*a, a + 1), (2*a + 1, 2*a, a + 1), (1, 2*a, a + 1), (0, a + 1, a + 1), (2*a, a + 1, a + 1), (a + 1, a + 1, a + 1), (a + 2, a + 1, a + 1), (2, a + 1, a + 1), (a, a + 1, a + 1), (2*a + 2, a + 1, a + 1), (2*a + 1, a + 1, a + 1), (1, a + 1, a + 1), (0, a + 2, a + 1), (2*a, a + 2, a + 1), (a + 1, a + 2, a + 1), (a + 2, a + 2, a + 1), (2, a + 2, a + 1), (a, a + 2, a + 1), (2*a + 2, a + 2, a + 1), (2*a + 1, a + 2, a + 1), (1, a + 2, a + 1), (0, 2, a + 1), (2*a, 2, a + 1), (a + 1, 2, a + 1), (a + 2, 2, a + 1), (2, 2, a + 1), (a, 2, a + 1), (2*a + 2, 2, a + 1), (2*a + 1, 2, a + 1), (1, 2, a + 1), (0, a, a + 1), (2*a, a, a + 1), (a + 1, a, a + 1), (a + 2, a, a + 1), (2, a, a + 1), (a, a, a + 1), (2*a + 2, a, a + 1), (2*a + 1, a, a + 1), (1, a, a + 1), (0, 2*a + 2, a + 1), (2*a, 2*a + 2, a + 1), (a + 1, 2*a + 2, a + 1), (a + 2, 2*a + 2, a + 1), (2, 2*a + 2, a + 1), (a, 2*a + 2, a + 1), (2*a + 2, 2*a + 2, a + 1), (2*a + 1, 2*a + 2, a + 1), (1, 2*a + 2, a + 1), (0, 2*a + 1, a + 1), (2*a, 2*a + 1, a + 1), (a + 1, 2*a + 1, a + 1), (a + 2, 2*a + 1, a + 1), (2, 2*a + 1, a + 1), (a, 2*a + 1, a + 1), (2*a + 2, 2*a + 1, a + 1), (2*a + 1, 2*a + 1, a + 1), (1, 2*a + 1, a + 1), (0, 1, a + 1), (2*a, 1, a + 1), (a + 1, 1, a + 1), (a + 2, 1, a + 1), (2, 1, a + 1), (a, 1, a + 1), (2*a + 2, 1, a + 1), (2*a + 1, 1, a + 1), (1, 1, a + 1), (0, 0, a + 2), (2*a, 0, a + 2), (a + 1, 0, a + 2), (a + 2, 0, a + 2), (2, 0, a + 2), (a, 0, a + 2), (2*a + 2, 0, a + 2), (2*a + 1, 0, a + 2), (1, 0, a + 2), (0, 2*a, a + 2), (2*a, 2*a, a + 2), (a + 1, 2*a, a + 2), (a + 2, 2*a, a + 2), (2, 2*a, a + 2), (a, 2*a, a + 2), (2*a + 2, 2*a, a + 2), (2*a + 1, 2*a, a + 2), (1, 2*a, a + 2), (0, a + 1, a + 2), (2*a, a + 1, a + 2), (a + 1, a + 1, a + 2), (a + 2, a + 1, a + 2), (2, a + 1, a + 2), (a, a + 1, a + 2), (2*a + 2, a + 1, a + 2), (2*a + 1, a + 1, a + 2), (1, a + 1, a + 2), (0, a + 2, a + 2), (2*a, a + 2, a + 2), (a + 1, a + 2, a + 2), (a + 2, a + 2, a + 2), (2, a + 2, a + 2), (a, a + 2, a + 2), (2*a + 2, a + 2, a + 2), (2*a + 1, a + 2, a + 2), (1, a + 2, a + 2), (0, 2, a + 2), (2*a, 2, a + 2), (a + 1, 2, a + 2), (a + 2, 2, a + 2), (2, 2, a + 2), (a, 2, a + 2), (2*a + 2, 2, a + 2), (2*a + 1, 2, a + 2), (1, 2, a + 2), (0, a, a + 2), (2*a, a, a + 2), (a + 1, a, a + 2), (a + 2, a, a + 2), (2, a, a + 2), (a, a, a + 2), (2*a + 2, a, a + 2), (2*a + 1, a, a + 2), (1, a, a + 2), (0, 2*a + 2, a + 2), (2*a, 2*a + 2, a + 2), (a + 1, 2*a + 2, a + 2), (a + 2, 2*a + 2, a + 2), (2, 2*a + 2, a + 2), (a, 2*a + 2, a + 2), (2*a + 2, 2*a + 2, a + 2), (2*a + 1, 2*a + 2, a + 2), (1, 2*a + 2, a + 2), (0, 2*a + 1, a + 2), (2*a, 2*a + 1, a + 2), (a + 1, 2*a + 1, a + 2), (a + 2, 2*a + 1, a + 2), (2, 2*a + 1, a + 2), (a, 2*a + 1, a + 2), (2*a + 2, 2*a + 1, a + 2), (2*a + 1, 2*a + 1, a + 2), (1, 2*a + 1, a + 2), (0, 1, a + 2), (2*a, 1, a + 2), (a + 1, 1, a + 2), (a + 2, 1, a + 2), (2, 1, a + 2), (a, 1, a + 2), (2*a + 2, 1, a + 2), (2*a + 1, 1, a + 2), (1, 1, a + 2), (0, 0, 2), (2*a, 0, 2), (a + 1, 0, 2), (a + 2, 0, 2), (2, 0, 2), (a, 0, 2), (2*a + 2, 0, 2), (2*a + 1, 0, 2), (1, 0, 2), (0, 2*a, 2), (2*a, 2*a, 2), (a + 1, 2*a, 2), (a + 2, 2*a, 2), (2, 2*a, 2), (a, 2*a, 2), (2*a + 2, 2*a, 2), (2*a + 1, 2*a, 2), (1, 2*a, 2), (0, a + 1, 2), (2*a, a + 1, 2), (a + 1, a + 1, 2), (a + 2, a + 1, 2), (2, a + 1, 2), (a, a + 1, 2), (2*a + 2, a + 1, 2), (2*a + 1, a + 1, 2), (1, a + 1, 2), (0, a + 2, 2), (2*a, a + 2, 2), (a + 1, a + 2, 2), (a + 2, a + 2, 2), (2, a + 2, 2), (a, a + 2, 2), (2*a + 2, a + 2, 2), (2*a + 1, a + 2, 2), (1, a + 2, 2), (0, 2, 2), (2*a, 2, 2), (a + 1, 2, 2), (a + 2, 2, 2), (2, 2, 2), (a, 2, 2), (2*a + 2, 2, 2), (2*a + 1, 2, 2), (1, 2, 2), (0, a, 2), (2*a, a, 2), (a + 1, a, 2), (a + 2, a, 2), (2, a, 2), (a, a, 2), (2*a + 2, a, 2), (2*a + 1, a, 2), (1, a, 2), (0, 2*a + 2, 2), (2*a, 2*a + 2, 2), (a + 1, 2*a + 2, 2), (a + 2, 2*a + 2, 2), (2, 2*a + 2, 2), (a, 2*a + 2, 2), (2*a + 2, 2*a + 2, 2), (2*a + 1, 2*a + 2, 2), (1, 2*a + 2, 2), (0, 2*a + 1, 2), (2*a, 2*a + 1, 2), (a + 1, 2*a + 1, 2), (a + 2, 2*a + 1, 2), (2, 2*a + 1, 2), (a, 2*a + 1, 2), (2*a + 2, 2*a + 1, 2), (2*a + 1, 2*a + 1, 2), (1, 2*a + 1, 2), (0, 1, 2), (2*a, 1, 2), (a + 1, 1, 2), (a + 2, 1, 2), (2, 1, 2), (a, 1, 2), (2*a + 2, 1, 2), (2*a + 1, 1, 2), (1, 1, 2), (0, 0, a), (2*a, 0, a), (a + 1, 0, a), (a + 2, 0, a), (2, 0, a), (a, 0, a), (2*a + 2, 0, a), (2*a + 1, 0, a), (1, 0, a), (0, 2*a, a), (2*a, 2*a, a), (a + 1, 2*a, a), (a + 2, 2*a, a), (2, 2*a, a), (a, 2*a, a), (2*a + 2, 2*a, a), (2*a + 1, 2*a, a), (1, 2*a, a), (0, a + 1, a), (2*a, a + 1, a), (a + 1, a + 1, a), (a + 2, a + 1, a), (2, a + 1, a), (a, a + 1, a), (2*a + 2, a + 1, a), (2*a + 1, a + 1, a), (1, a + 1, a), (0, a + 2, a), (2*a, a + 2, a), (a + 1, a + 2, a), (a + 2, a + 2, a), (2, a + 2, a), (a, a + 2, a), (2*a + 2, a + 2, a), (2*a + 1, a + 2, a), (1, a + 2, a), (0, 2, a), (2*a, 2, a), (a + 1, 2, a), (a + 2, 2, a), (2, 2, a), (a, 2, a), (2*a + 2, 2, a), (2*a + 1, 2, a), (1, 2, a), (0, a, a), (2*a, a, a), (a + 1, a, a), (a + 2, a, a), (2, a, a), (a, a, a), (2*a + 2, a, a), (2*a + 1, a, a), (1, a, a), (0, 2*a + 2, a), (2*a, 2*a + 2, a), (a + 1, 2*a + 2, a), (a + 2, 2*a + 2, a), (2, 2*a + 2, a), (a, 2*a + 2, a), (2*a + 2, 2*a + 2, a), (2*a + 1, 2*a + 2, a), (1, 2*a + 2, a), (0, 2*a + 1, a), (2*a, 2*a + 1, a), (a + 1, 2*a + 1, a), (a + 2, 2*a + 1, a), (2, 2*a + 1, a), (a, 2*a + 1, a), (2*a + 2, 2*a + 1, a), (2*a + 1, 2*a + 1, a), (1, 2*a + 1, a), (0, 1, a), (2*a, 1, a), (a + 1, 1, a), (a + 2, 1, a), (2, 1, a), (a, 1, a), (2*a + 2, 1, a), (2*a + 1, 1, a), (1, 1, a), (0, 0, 2*a + 2), (2*a, 0, 2*a + 2), (a + 1, 0, 2*a + 2), (a + 2, 0, 2*a + 2), (2, 0, 2*a + 2), (a, 0, 2*a + 2), (2*a + 2, 0, 2*a + 2), (2*a + 1, 0, 2*a + 2), (1, 0, 2*a + 2), (0, 2*a, 2*a + 2), (2*a, 2*a, 2*a + 2), (a + 1, 2*a, 2*a + 2), (a + 2, 2*a, 2*a + 2), (2, 2*a, 2*a + 2), (a, 2*a, 2*a + 2), (2*a + 2, 2*a, 2*a + 2), (2*a + 1, 2*a, 2*a + 2), (1, 2*a, 2*a + 2), (0, a + 1, 2*a + 2), (2*a, a + 1, 2*a + 2), (a + 1, a + 1, 2*a + 2), (a + 2, a + 1, 2*a + 2), (2, a + 1, 2*a + 2), (a, a + 1, 2*a + 2), (2*a + 2, a + 1, 2*a + 2), (2*a + 1, a + 1, 2*a + 2), (1, a + 1, 2*a + 2), (0, a + 2, 2*a + 2), (2*a, a + 2, 2*a + 2), (a + 1, a + 2, 2*a + 2), (a + 2, a + 2, 2*a + 2), (2, a + 2, 2*a + 2), (a, a + 2, 2*a + 2), (2*a + 2, a + 2, 2*a + 2), (2*a + 1, a + 2, 2*a + 2), (1, a + 2, 2*a + 2), (0, 2, 2*a + 2), (2*a, 2, 2*a + 2), (a + 1, 2, 2*a + 2), (a + 2, 2, 2*a + 2), (2, 2, 2*a + 2), (a, 2, 2*a + 2), (2*a + 2, 2, 2*a + 2), (2*a + 1, 2, 2*a + 2), (1, 2, 2*a + 2), (0, a, 2*a + 2), (2*a, a, 2*a + 2), (a + 1, a, 2*a + 2), (a + 2, a, 2*a + 2), (2, a, 2*a + 2), (a, a, 2*a + 2), (2*a + 2, a, 2*a + 2), (2*a + 1, a, 2*a + 2), (1, a, 2*a + 2), (0, 2*a + 2, 2*a + 2), (2*a, 2*a + 2, 2*a + 2), (a + 1, 2*a + 2, 2*a + 2), (a + 2, 2*a + 2, 2*a + 2), (2, 2*a + 2, 2*a + 2), (a, 2*a + 2, 2*a + 2), (2*a + 2, 2*a + 2, 2*a + 2), (2*a + 1, 2*a + 2, 2*a + 2), (1, 2*a + 2, 2*a + 2), (0, 2*a + 1, 2*a + 2), (2*a, 2*a + 1, 2*a + 2), (a + 1, 2*a + 1, 2*a + 2), (a + 2, 2*a + 1, 2*a + 2), (2, 2*a + 1, 2*a + 2), (a, 2*a + 1, 2*a + 2), (2*a + 2, 2*a + 1, 2*a + 2), (2*a + 1, 2*a + 1, 2*a + 2), (1, 2*a + 1, 2*a + 2), (0, 1, 2*a + 2), (2*a, 1, 2*a + 2), (a + 1, 1, 2*a + 2), (a + 2, 1, 2*a + 2), (2, 1, 2*a + 2), (a, 1, 2*a + 2), (2*a + 2, 1, 2*a + 2), (2*a + 1, 1, 2*a + 2), (1, 1, 2*a + 2), (0, 0, 2*a + 1), (2*a, 0, 2*a + 1), (a + 1, 0, 2*a + 1), (a + 2, 0, 2*a + 1), (2, 0, 2*a + 1), (a, 0, 2*a + 1), (2*a + 2, 0, 2*a + 1), (2*a + 1, 0, 2*a + 1), (1, 0, 2*a + 1), (0, 2*a, 2*a + 1), (2*a, 2*a, 2*a + 1), (a + 1, 2*a, 2*a + 1), (a + 2, 2*a, 2*a + 1), (2, 2*a, 2*a + 1), (a, 2*a, 2*a + 1), (2*a + 2, 2*a, 2*a + 1), (2*a + 1, 2*a, 2*a + 1), (1, 2*a, 2*a + 1), (0, a + 1, 2*a + 1), (2*a, a + 1, 2*a + 1), (a + 1, a + 1, 2*a + 1), (a + 2, a + 1, 2*a + 1), (2, a + 1, 2*a + 1), (a, a + 1, 2*a + 1), (2*a + 2, a + 1, 2*a + 1), (2*a + 1, a + 1, 2*a + 1), (1, a + 1, 2*a + 1), (0, a + 2, 2*a + 1), (2*a, a + 2, 2*a + 1), (a + 1, a + 2, 2*a + 1), (a + 2, a + 2, 2*a + 1), (2, a + 2, 2*a + 1), (a, a + 2, 2*a + 1), (2*a + 2, a + 2, 2*a + 1), (2*a + 1, a + 2, 2*a + 1), (1, a + 2, 2*a + 1), (0, 2, 2*a + 1), (2*a, 2, 2*a + 1), (a + 1, 2, 2*a + 1), (a + 2, 2, 2*a + 1), (2, 2, 2*a + 1), (a, 2, 2*a + 1), (2*a + 2, 2, 2*a + 1), (2*a + 1, 2, 2*a + 1), (1, 2, 2*a + 1), (0, a, 2*a + 1), (2*a, a, 2*a + 1), (a + 1, a, 2*a + 1), (a + 2, a, 2*a + 1), (2, a, 2*a + 1), (a, a, 2*a + 1), (2*a + 2, a, 2*a + 1), (2*a + 1, a, 2*a + 1), (1, a, 2*a + 1), (0, 2*a + 2, 2*a + 1), (2*a, 2*a + 2, 2*a + 1), (a + 1, 2*a + 2, 2*a + 1), (a + 2, 2*a + 2, 2*a + 1), (2, 2*a + 2, 2*a + 1), (a, 2*a + 2, 2*a + 1), (2*a + 2, 2*a + 2, 2*a + 1), (2*a + 1, 2*a + 2, 2*a + 1), (1, 2*a + 2, 2*a + 1), (0, 2*a + 1, 2*a + 1), (2*a, 2*a + 1, 2*a + 1), (a + 1, 2*a + 1, 2*a + 1), (a + 2, 2*a + 1, 2*a + 1), (2, 2*a + 1, 2*a + 1), (a, 2*a + 1, 2*a + 1), (2*a + 2, 2*a + 1, 2*a + 1), (2*a + 1, 2*a + 1, 2*a + 1), (1, 2*a + 1, 2*a + 1), (0, 1, 2*a + 1), (2*a, 1, 2*a + 1), (a + 1, 1, 2*a + 1), (a + 2, 1, 2*a + 1), (2, 1, 2*a + 1), (a, 1, 2*a + 1), (2*a + 2, 1, 2*a + 1), (2*a + 1, 1, 2*a + 1), (1, 1, 2*a + 1), (0, 0, 1), (2*a, 0, 1), (a + 1, 0, 1), (a + 2, 0, 1), (2, 0, 1), (a, 0, 1), (2*a + 2, 0, 1), (2*a + 1, 0, 1), (1, 0, 1), (0, 2*a, 1), (2*a, 2*a, 1), (a + 1, 2*a, 1), (a + 2, 2*a, 1), (2, 2*a, 1), (a, 2*a, 1), (2*a + 2, 2*a, 1), (2*a + 1, 2*a, 1), (1, 2*a, 1), (0, a + 1, 1), (2*a, a + 1, 1), (a + 1, a + 1, 1), (a + 2, a + 1, 1), (2, a + 1, 1), (a, a + 1, 1), (2*a + 2, a + 1, 1), (2*a + 1, a + 1, 1), (1, a + 1, 1), (0, a + 2, 1), (2*a, a + 2, 1), (a + 1, a + 2, 1), (a + 2, a + 2, 1), (2, a + 2, 1), (a, a + 2, 1), (2*a + 2, a + 2, 1), (2*a + 1, a + 2, 1), (1, a + 2, 1), (0, 2, 1), (2*a, 2, 1), (a + 1, 2, 1), (a + 2, 2, 1), (2, 2, 1), (a, 2, 1), (2*a + 2, 2, 1), (2*a + 1, 2, 1), (1, 2, 1), (0, a, 1), (2*a, a, 1), (a + 1, a, 1), (a + 2, a, 1), (2, a, 1), (a, a, 1), (2*a + 2, a, 1), (2*a + 1, a, 1), (1, a, 1), (0, 2*a + 2, 1), (2*a, 2*a + 2, 1), (a + 1, 2*a + 2, 1), (a + 2, 2*a + 2, 1), (2, 2*a + 2, 1), (a, 2*a + 2, 1), (2*a + 2, 2*a + 2, 1), (2*a + 1, 2*a + 2, 1), (1, 2*a + 2, 1), (0, 2*a + 1, 1), (2*a, 2*a + 1, 1), (a + 1, 2*a + 1, 1), (a + 2, 2*a + 1, 1), (2, 2*a + 1, 1), (a, 2*a + 1, 1), (2*a + 2, 2*a + 1, 1), (2*a + 1, 2*a + 1, 1), (1, 2*a + 1, 1), (0, 1, 1), (2*a, 1, 1), (a + 1, 1, 1), (a + 2, 1, 1), (2, 1, 1), (a, 1, 1), (2*a + 2, 1, 1), (2*a + 1, 1, 1), (1, 1, 1)]
}}}

<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]
///
Vector space of degree 3 and dimension 2 over Finite Field in a of size 3^2
Basis matrix:
[      1       0 2*a + 2]
[      0       1       2]
}}}

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

{{{id=104|
all_subs[86]
///
Vector space of degree 3 and dimension 2 over Finite Field in a of size 3^2
Basis matrix:
[1 a 0]
[0 0 1]
}}}

{{{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|

///
}}}