AIM Basic Sage Linear Algebra

1184 days ago by pub

Sage Tutorial

By Jason Grout and Leslie Hogben

 
       
# Any line that starts with a number sign is a comment - this is an example of a comment. # Comments are purely explanatory; they are not executed 
       

If you can't see everything in a cell on a Mac, drag the lower right corner to make it bigger. Generally the cells should automatically resize to fit the contents, though.

The Sage notebook is based on cells (those are the boxes on the page that contain code) and chronological. This means that Sage only cares about the order in which you execute things, not the order on the screen. You do not have to execute things in order on the screen.

To execute (evaluate) a cell, click in the cell and either click the "execute" link that appears just below the cell or press the SHIFT and RETURN keys at the same time.

For example, evaluate the following cell to ask Sage what 2+3*4 is.

2+3*4 
       
14
14

We can assign these values to variables using a single equals sign. Evaluate the following cell to give y the value 2+3*4.

y = 2+3*4 
       

To see the value of a variable, just type the variable alone.

       
14
14

You can combine multiple computations in the same cell by putting each statement on a new line. Here, we assign y to be a value, then print ask Sage to print out the value. In general, Sage will print out whatever the last command in the cell returns.

Go ahead and change the values to have Sage compute with some of your favorite numbers. You can evaluate the cell over and over again.

y=6*2-5 y 
       
7
7

To make a new cell, move the mouse between cells. When you see a horizontal blue line, click and a new cell will appear. Try making a new cell and doing a computation.

 
       

Linear Algebra

Now we will enter some matrices and do some arithmetic with them. First, we'll create some matrices and assign them to variables. We create a matrix by giving the matrix command a list of rows. Lists in Sage and Python are enclosed in square brackets, so we specify a list of rows as follows.

m1 = matrix([[1,2],[3,4]]) m1 
       
[1 2]
[3 4]
[1 2]
[3 4]
m2=matrix([[7,-3],[0,5]]) m2 
       
[ 7 -3]
[ 0  5]
[ 7 -3]
[ 0  5]

Matrix arithmetic is pretty self-explanatory.

m1+m2 
       
[ 8 -1]
[ 3  9]
[ 8 -1]
[ 3  9]
m1*m2 
       
[ 7  7]
[21 11]
[ 7  7]
[21 11]
m1^3 
       
[ 37  54]
[ 81 118]
[ 37  54]
[ 81 118]

There are several ways to get the inverse.

m1^(-1) 
       
[  -2    1]
[ 3/2 -1/2]
[  -2    1]
[ 3/2 -1/2]
~m1 
       
[  -2    1]
[ 3/2 -1/2]
[  -2    1]
[ 3/2 -1/2]

Let's check that that is actually the inverse.

m1*(~m1) 
       
[1 0]
[0 1]
[1 0]
[0 1]
m1*(~m1)==identity_matrix(2) 
       
True
True
 
       
 
       

We can get the matrix entries as follows. Following a convention in computers, Sage starts its numbering at 0, so the upper left corner of the matrix is the (0,0) entry.

m1[0,0] 
       
1
1

We can also set a matrix entry using this notation.

m1[0,0] = 2 m1 
       
[2 2]
[3 4]
[2 2]
[3 4]

And we can get not just matrix entries, but whole submatrices using this convention. To get a submatrix, just give a list instead of a single row or column number. In the following we get entries in the 0th or 1st rows and the 0th column.

m1[[0,1],0] 
       
[2]
[3]
[2]
[3]

We can also do more sophisticated matrix calculations, like rank, transpose, etc.

rank(m1) 
       
2
2
transpose(m1) 
       
[2 3]
[2 4]
[2 3]
[2 4]
det(m1) 
       
2
2

We can get the characteristic polynomial as well.

characteristic_polynomial(m1) 
       
x^2 - 6*x + 2
x^2 - 6*x + 2

Or you can use the shortened form "charpoly"

charpoly(m1) 
       
x^2 - 6*x + 2
x^2 - 6*x + 2

There are tens of thousands of possible functions and hundreds of mathematical structures that Sage understands. However, any specific mathematical structure only makes sense to a handful of functions. For example, it hardly makes sense to talk about the transpose of a finite field, does it?

To deal with this, Sage adopts a convention in computer science that the mathematical objects in Sage are smart; they all come with the functions that make sense attached to them. You can call the functions that are attached to an object by typing the object, then a period, then the function name, followed by parentheses. For example, to compute some of the previous results using this, we do the following computations.

m1.determinant() 
       
2
2
m1.charpoly() 
       
x^2 - 6*x + 2
x^2 - 6*x + 2
m1.inverse() 
       
[   2   -1]
[-3/2    1]
[   2   -1]
[-3/2    1]

The way to think about this notation is to imagine yourself asking the mathematical object a question or requesting that it perform some function. For example, in the previous example, we ask the matrix m1 to tell us what its inverse is.

This convention helps us understand what sorts of things we can do to an object we have. To get a list of all the functions attached to an object, just type the object, a period, and then press tab. If you start typing a function name and then press tab, it will just give you the functions that start with the letters you've already typed.

Try seeing what functions Sage offers on integer matrices by typing "m1." and then pressing tab in the following cell.

 
       

You'll notice that there are lots and lots of functions attached to integer matrices.  Some of these functions are only available using this dot notation.

You can see the documentation for a function by following it with a questionmark and then evaluating the cell.

m1.determinant? 
       

Putting two question marks in gives the actual code to the function.  You can examine the code to any function in Sage.  You can also modify it and improve it, if you'd like.  Many undergraduate and graduate students are active developers of Sage.

m1.determinant?? 
       
 
       

Sage can deal with matrices over different types of numbers.  By default, Sage tries to guess from the entries of a matrix, so if you enter only integers, Sage treats the matrix as a matrix over integers. 

To specifically tell Sage to make the matrix over the rational numbers, for example, use QQ.

mQ = matrix(QQ, [[1,2,3], [4,5,6], [7,8,9]]) mQ 
       
[1 2 3]
[4 5 6]
[7 8 9]
[1 2 3]
[4 5 6]
[7 8 9]

Now we compute the row-reduced echelon form. For right now, it is important that we defined the matrix to be over the rationals; otherwise Sage would have computed a special type of echelon form without using division (i.e., just using integer operations).

mQ.echelon_form() 
       
[ 1  0 -1]
[ 0  1  2]
[ 0  0  0]
[ 1  0 -1]
[ 0  1  2]
[ 0  0  0]

In order to compute the row-reduced echelon form of a matrix that is over the integers (for now), like our matrix m1 above, we can change the ring of the matrix. Notice that we can "chain" functions together, calling them one right after another.

m1.change_ring(QQ).echelon_form() 
       
[1 0]
[0 1]
[1 0]
[0 1]

To print an intermediate result in the middle of a cell, use the print command like the following. A print command by itself just puts in a blank line. You can also use the "show" command to display a nice, mathematical-looking form of your object.

print mQ print show(mQ[[0,2],[0,2]]) 
       
[1 2 3]
[4 5 6]
[7 8 9]

\left(\begin{array}{rr} 1 & 3 \\ 7 & 9 \end{array}\right)
[1 2 3]
[4 5 6]
[7 8 9]

\left(\begin{array}{rr} 1 & 3 \\ 7 & 9 \end{array}\right)

Sage can also compute eigenvalues.

print m2 print m2.eigenvalues() 
       
[ 7 -3]
[ 0  5]

[7, 5]
[ 7 -3]
[ 0  5]

[7, 5]
print m1 print m1.eigenvalues() 
       
[2 2]
[3 4]

[0.3542486889354094?, 5.645751311064590?]
[2 2]
[3 4]

[0.3542486889354094?, 5.645751311064590?]

Do you notice the question marks at the end of the numbers above? These indicate that Sage is treating these numbers in a rather sophisticated way. Even though Sage is printing out a numerical approximation of the eigenvalue, Sage knows what the exact value is and uses that exact value in computations. In effect, while Sage is only showing you a numeric approximation, Sage is using infinite precision in the computations.

Let's check that those two eigenvalues above multiply to be exactly the determinant, with no roundoff error. The prod command takes a list of numbers and returns their product. The double equals sign is an equality test and returns True if the two things are exactly equal and False otherwise.

print m1.det() m1.det()==prod(m1.eigenvalues()) 
       
2
True
2
True

Of course, Sage can compute both right and left eigenvectors. The eigenvectors_right command gives the eigenvalue, a basis of eigenvectors for the eigenvalue, and algebraic multiplicity for the eigenvalue.

show(m2.eigenvectors_right()) 
       
\begin{array}{l}[\left(7, \left[\left(1,0\right)\right], 1\right),\\ \left(5, \left[\left(1,\frac{2}{3}\right)\right], 1\right)]\end{array}
\begin{array}{l}[\left(7, \left[\left(1,0\right)\right], 1\right),\\ \left(5, \left[\left(1,\frac{2}{3}\right)\right], 1\right)]\end{array}

You can also ask for the matrices of eigenvectors and a diagonal matrix of eigenvalues.

show(m2.eigenmatrix_right()) 
       
\left(\left(\begin{array}{rr} 7 & 0 \\ 0 & 5 \end{array}\right), \left(\begin{array}{rr} 1 & 1 \\ 0 & \frac{2}{3} \end{array}\right)\right)
\left(\left(\begin{array}{rr} 7 & 0 \\ 0 & 5 \end{array}\right), \left(\begin{array}{rr} 1 & 1 \\ 0 & \frac{2}{3} \end{array}\right)\right)

To get the list of eigenvectors just ask for the columns of the eigenvector matrix.

evals, evecs = m2.eigenmatrix_right() evecs.columns() 
       
[(1, 0), (1, 2/3)]
[(1, 0), (1, 2/3)]

Let's check that these are eigenvectors by doing the appropriate multiplication.

evec1, evec2 = evecs.columns() eval1, eval2 = m2.eigenvalues() print "First eigenvector: ", evec1 print "First eigenvalue: ", eval1 print m2*evec1 print eval1*evec1 
       
First eigenvector:  (1, 0)
First eigenvalue:  7
(7, 0)
(7, 0)
First eigenvector:  (1, 0)
First eigenvalue:  7
(7, 0)
(7, 0)

Let's create a random matrix over the rationals and see the general structure of the matrix.

a=random_matrix(QQ,100) a.visualize_structure() 
       
a.echelon_form().visualize_structure() 
       

Graphs

# graphs g=Graph({1:[3,4,5,6],2:[8,10,12],6:[7,8],8:[9],10:[11]}) show(g) 
       
m = g.adjacency_matrix() m 
       
[0 0 1 1 1 1 0 0 0 0 0 0]
[0 0 0 0 0 0 0 1 0 1 0 1]
[1 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 1 1 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0 0 0]
[0 1 0 0 0 1 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0 0 0]
[0 1 0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 1 0 0]
[0 1 0 0 0 0 0 0 0 0 0 0]
[0 0 1 1 1 1 0 0 0 0 0 0]
[0 0 0 0 0 0 0 1 0 1 0 1]
[1 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 0 0 0 0 0 0]
[1 0 0 0 0 0 1 1 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0 0 0]
[0 1 0 0 0 1 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0 0 0]
[0 1 0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 1 0 0]
[0 1 0 0 0 0 0 0 0 0 0 0]
gg=Graph(m) show(gg) 
       
# Example of a for loop. Note you must use indentation or it won't work. for i in [1..3]: print m1-i*identity_matrix(2) print