Using Sage for Computations in Linear Algebra¶

Sage is a more generic math package, made to replace Mathematica, Maple, MatLab, etc. It is built on top of Python.

You can use it online with Cocalc. (You can use a free account.)

For simple computations, you can use the Sage Cell Server.

Notebook: If you have how to run Jupyter Notebooks with Sage, you can download the notebook here.

Matrices¶

Let's create the matrix $$ A = \begin{bmatrix} 2 & 8 & 4 \\ 2 & 5 & 1 \\ 4 & 10 & -1 \end{bmatrix} $$

In [1]:
A = matrix([[2, 8, 4], [2, 5, 1], [4, 10, -1]])  # enter a list of rows (each row also a list)
A
Out[1]:
[ 2  8  4]
[ 2  5  1]
[ 4 10 -1]

We can also do:

In [2]:
A = matrix(3, 3, [2, 8, 4, 2, 5, 1, 4, 10, -1])  # enter the dimensions and list
A
Out[2]:
[ 2  8  4]
[ 2  5  1]
[ 4 10 -1]

Or even

In [3]:
A = matrix(3, [2, 8, 4, 2, 5, 1, 4, 10, -1])  # for square matrices we can give only one dimension
A
Out[3]:
[ 2  8  4]
[ 2  5  1]
[ 4 10 -1]

We can also specify where the coefficients lie. For instance, to have the same matrix over $\mathbb{F}_5 = \mathbb{Z}/5\mathbb{Z}$ (i.e., the field with $5$ elements), we can do:

In [4]:
A = matrix(FiniteField(5), 3, [2, 8, 4, 2, 5, 1, 4, 10, -1])  # for square matrices we can give only one dimension
A
Out[4]:
[2 3 4]
[2 0 1]
[4 0 4]

(Note how the entries are reduced modulo $5$!)

We can call the function latex to produce a $\LaTeX$ representation of our matrix:

In [5]:
latex(A)
Out[5]:
\left(\begin{array}{rrr}
2 & 3 & 4 \\
2 & 0 & 1 \\
4 & 0 & 4
\end{array}\right)

Special Matrices¶

We can create the identity matrix of any dimension:

In [6]:
identity_matrix(4)
Out[6]:
[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]

We can also specify the type of coefficients. For instance, for coefficients in $\mathbb{F}_3 = \mathbb{Z} / 3 \mathbb{Z}$:

In [7]:
identity_matrix(FiniteField(3), 2)
Out[7]:
[1 0]
[0 1]

Similarly for zero matrices:

In [8]:
zero_matrix(3, 4)
Out[8]:
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]
In [9]:
zero_matrix(CC, 3, 2)  # CC = complex numbers (using floats)
Out[9]:
[0.000000000000000 0.000000000000000]
[0.000000000000000 0.000000000000000]
[0.000000000000000 0.000000000000000]

We can create random matrices. For finite entry ring (or field), it is as you'd expect:

In [10]:
random_matrix(FiniteField(11), 3, 4)
Out[10]:
[1 7 0 9]
[6 6 1 9]
[1 0 2 8]

(Run the cell above a few times to see different random results.)

For $\mathbb{Q}$ we can specifiy bounds for the numerator and denominator:

In [11]:
random_matrix(QQ, 4, 5, num_bound=10, den_bound=10)
Out[11]:
[ 3/4 -7/2 -5/2  -10   -3]
[-7/5    5  5/8    0   -2]
[-1/3  3/7 -2/3  3/8 -1/9]
[-9/7 -2/5 -3/2 -8/5 -5/7]

For $\mathbb{R}$ (RR in Sage), we can specify the maximum and minimum:

In [12]:
random_matrix(RR, 2, 2, min=-10, max=10)
Out[12]:
[0.891317069574020  8.43768596822291]
[-1.60966517026439 -9.00857316074933]

Complex Matrices¶

There a few implementations of complex matrices in Sage: some "numeric" (using floats/decimals) and one "symbolic" (with exact rational numbers and $i$). The latter uses I for the imaginary number $i$:

In [13]:
I^2
Out[13]:
-1

So, we can make matrices such as:

In [14]:
A = matrix(2, 2, [I/2, 0, 1/4 + 3*I, -3])
A
Out[14]:
[    1/2*I         0]
[3*I + 1/4        -3]

We can use CC to use floats:

In [15]:
B = matrix(CC, 2, 2, [I/2, 0, 1/4 + 3*I, -3])
B
Out[15]:
[                   0.500000000000000*I                      0.000000000000000]
[0.250000000000000 + 3.00000000000000*I                      -3.00000000000000]

If one wants to specify the "symbolic" complex numbers, we can save if in a variable and use:

In [16]:
CCs = I.parent()  # save the "symbolic complex field" in CCs
CCs
Out[16]:
Number Field in I with defining polynomial x^2 + 1 with I = 1*I
In [17]:
C = matrix(CCs, 2, 2, [1, 0, 1/3, I])
C
Out[17]:
[  1   0]
[1/3   I]

Symbolic Matrices¶

If we want to use symbolic computations, we need to introduce the variables:

In [18]:
var("a11, a12, a21, a22")
A = matrix(2, 2, [a11, a12, a21, a22])
A
Out[18]:
[a11 a12]
[a21 a22]
In [19]:
A.det()
Out[19]:
-a12*a21 + a11*a22
In [20]:
A^2
Out[20]:
[  a11^2 + a12*a21 a11*a12 + a12*a22]
[a11*a21 + a21*a22   a12*a21 + a22^2]
In [21]:
A^(-1)
Out[21]:
[1/a11 - a12*a21/(a11^2*(a12*a21/a11 - a22))               a12/(a11*(a12*a21/a11 - a22))]
[              a21/(a11*(a12*a21/a11 - a22))                      -1/(a12*a21/a11 - a22)]

Matrix Operations¶

We can add matrices:

In [22]:
A = random_matrix(QQ, 3, 2, num_bound=5, den_bound=5)
B = random_matrix(QQ, 3, 2, num_bound=5, den_bound=5)
In [23]:
A
Out[23]:
[ 0 -1]
[-3 -1]
[ 2  1]
In [24]:
B
Out[24]:
[-3/5  1/2]
[ 1/4 -1/4]
[   0    1]
In [25]:
A + B
Out[25]:
[ -3/5  -1/2]
[-11/4  -5/4]
[    2     2]

We can also do scalar multiplication:

In [26]:
1 / 3 * A
Out[26]:
[   0 -1/3]
[  -1 -1/3]
[ 2/3  1/3]

We can multiply matrices as expected:

In [27]:
A = matrix(3, [2, 8, 4, 2, 5, 1, 4, 10, -1])  # for square matrices we can give only one dimension
B = matrix(3, 2, [1, 2, 0, 1, 1, 1])
In [28]:
A
Out[28]:
[ 2  8  4]
[ 2  5  1]
[ 4 10 -1]
In [29]:
B
Out[29]:
[1 2]
[0 1]
[1 1]
In [30]:
A * B
Out[30]:
[ 6 16]
[ 3 10]
[ 3 17]

We can transpose:

In [31]:
A.T
Out[31]:
[ 2  2  4]
[ 8  5 10]
[ 4  1 -1]

We can take determinants:

In [32]:
A.det()
Out[32]:
18

We can invert, if possible:

In [33]:
A^(-1)
Out[33]:
[-5/6  8/3 -2/3]
[ 1/3   -1  1/3]
[   0  2/3 -1/3]

(Note that the result is made of rational numbers, not floats (i.e., decimals).)

Echelon Form¶

We can get the reduced row echelon form:

In [34]:
A.echelon_form()
Out[34]:
[2 2 1]
[0 3 0]
[0 0 3]

This result seems unexpected, since we don't have leading ones! This is because Sage is staying with the integers for entries! If we want to allow fractions, we must declare the matrix over the rationals $\mathbb{Q}$, which in Sage is QQ:

In [35]:
A = matrix(QQ, 3, [2, 8, 4, 2, 5, 1, 4, 10, -1])  # declare rational entries!!
A
Out[35]:
[ 2  8  4]
[ 2  5  1]
[ 4 10 -1]
In [36]:
A.echelon_form()
Out[36]:
[1 0 0]
[0 1 0]
[0 0 1]

Vectors¶

We also have vectors:

In [37]:
v = vector([1, 2, 3])
v
Out[37]:
(1, 2, 3)

We can then use if for matrix multiplication:

In [38]:
A * v
Out[38]:
(30, 15, 21)

Note that it knows that the vector must be interpreted as a column matrix!

But we can also do:

In [39]:
v * A
Out[39]:
(18, 48, 3)

Now, it is interpreted as a row matrix.

Row Operations¶

We can perform row operations.

Remember that Python starts indexing from $0$, not $1$

In [40]:
A = matrix(QQ, 3, range(9))
A
Out[40]:
[0 1 2]
[3 4 5]
[6 7 8]

Let's add $-3$ times the first row to the second:

In [41]:
A.add_multiple_of_row(1, 0, -3)
A
Out[41]:
[ 0  1  2]
[ 3  1 -1]
[ 6  7  8]

Note that this changes the original matrix!

Let's multiply the first row by $5$:

In [42]:
A.rescale_row(0, 5)
A
Out[42]:
[ 0  5 10]
[ 3  1 -1]
[ 6  7  8]

Let's swap the first and third rows:

In [43]:
A.swap_rows(0, 2)
A
Out[43]:
[ 6  7  8]
[ 3  1 -1]
[ 0  5 10]

Row Equivalency¶

There is no way to check directly if two matrices are row equivalent, i.e., one can be obtained from the other by row operations. But we can still do it: $A$ and $B$ are row equivalent if and only if they have the same reduced row echelon form.

In [44]:
A = matrix(QQ, 3, 4, range(12))
A
Out[44]:
[ 0  1  2  3]
[ 4  5  6  7]
[ 8  9 10 11]
In [45]:
B = matrix(QQ, 3, 4, range(2, 14))
B
Out[45]:
[ 2  3  4  5]
[ 6  7  8  9]
[10 11 12 13]
In [46]:
A.echelon_form() == B.echelon_form()
Out[46]:
True

So, $A$ and $B$ are row equivalent.

Let's check another matrix:

In [47]:
C = matrix(QQ, 3, 4, [i^2 for i in range(12)])
C
Out[47]:
[  0   1   4   9]
[ 16  25  36  49]
[ 64  81 100 121]
In [48]:
A.echelon_form() == C.echelon_form()
Out[48]:
False

So $A$ and $C$ are not row equivalent.

Matrix Spaces¶

Let's start with a $4 \times 3$ matrix $A$:

In [49]:
A = matrix(QQ, 4, 3, range(12))
A
Out[49]:
[ 0  1  2]
[ 3  4  5]
[ 6  7  8]
[ 9 10 11]

We can ask for the row space:

In [50]:
A.row_space()
Out[50]:
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -1]
[ 0  1  2]

We can ask for the column space:

In [51]:
A.column_space()
Out[51]:
Vector space of degree 4 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -1 -2]
[ 0  1  2  3]

We can ask for the rank:

In [52]:
A.rank()
Out[52]:
2

For the nullspace, we use right_kernel:

In [53]:
A.right_kernel()
Out[53]:
Vector space of degree 3 and dimension 1 over Rational Field
Basis matrix:
[ 1 -2  1]

Note that the method kernel gives the left kernel, meaning the space of solutions of $\vec{x} \cdot A = \vec{0}$!

In [54]:
A.kernel()
Out[54]:
Vector space of degree 4 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -3  2]
[ 0  1 -2  1]
In [55]:
A.left_kernel()
Out[55]:
Vector space of degree 4 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -3  2]
[ 0  1 -2  1]

(Note how left_kernel gives a subspace in $4$ dimensions, while right_kernel gives a subspace in $3$ dimensions, since $A$ is $4 \times 3$!)

So, for the nullity, we need to ask for right_nullity:

In [56]:
A.right_nullity()
Out[56]:
1
In [57]:
(A.T).right_kernel()
Out[57]:
Vector space of degree 4 and dimension 2 over Rational Field
Basis matrix:
[ 1  0 -3  2]
[ 0  1 -2  1]

Solving Numerical Linear Systems¶

Solving Numerical Linear Symbolically¶

Suppose we want to solve:

\begin{align*} 2x + 8y + 4z &= 2 \\ 2x + 5y + z &= 5 \\ 4x + 10y - z &= 1 \end{align*}

We can do:

In [58]:
var("x y z")  # declare the symbolic variables
eqs = [2*x + 8*y + 4*z == 2, 2*x + 5*y + z == 5, 4*x + 10*y -z == 1]  # make a list of equations
solve(eqs, (x, y, z))  # call solve
Out[58]:
[[x == 11, y == -4, z == 3]]

If the system has no solution, say

\begin{align*} x + y &= 1\\ x + y & = 2 \end{align*}
In [59]:
eqs = [x + y == 1, x + y == 2]
solve(eqs, (x, y))
Out[59]:
[]

It may have infinitely many, as in

\begin{align*} x + y &= 1\\ x + y & = 1 \end{align*}
In [60]:
eqs = [x + y == 1, x + y == 1]
solve(eqs, (x, y))
Out[60]:
[[x == -r1 + 1, y == r1]]

As you can see, it gives us a free parameter r1.

Solving Numerical Linear Using Matrices¶

We can also solve systems using matrices:

In [61]:
A = matrix(QQ, 3, 3, [2, 8, 4, 2, 5, 1, 4, 10, -1])
v = vector(QQ, [2, 5, 1])
A.solve_right(v)
Out[61]:
(11, -4, 3)

For the other examples:

In [62]:
B = matrix(QQ, 2, 4 * [1])
B
Out[62]:
[1 1]
[1 1]

We get an error when there is no solution:

In [63]:
v = vector(QQ, [1, 2])
B.solve_right(v)

---------------------------------------------------------------------------
NotFullRankError                          Traceback (most recent call last)
File ~/src/sage-10.2/src/sage/matrix/matrix2.pyx:942, in sage.matrix.matrix2.Matrix.solve_right()
    941 try:
--> 942     X = self._solve_right_nonsingular_square(C, check_rank=True)
    943 except NotFullRankError:

File ~/src/sage-10.2/src/sage/matrix/matrix2.pyx:988, in sage.matrix.matrix2.Matrix._solve_right_nonsingular_square()
    987 if check_rank and self.rank() < self.nrows():
--> 988     raise NotFullRankError
    989 D = self.augment(B)

NotFullRankError: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In [63], line 2
      1 v = vector(QQ, [Integer(1), Integer(2)])
----> 2 B.solve_right(v)

File ~/src/sage-10.2/src/sage/matrix/matrix2.pyx:944, in sage.matrix.matrix2.Matrix.solve_right()
    942         X = self._solve_right_nonsingular_square(C, check_rank=True)
    943     except NotFullRankError:
--> 944         X = self._solve_right_general(C, check=check)
    945 
    946 if b_is_vec:

File ~/src/sage-10.2/src/sage/matrix/matrix2.pyx:1062, in sage.matrix.matrix2.Matrix._solve_right_general()
   1060     # Have to check that we actually solved the equation.
   1061     if self*X != B:
-> 1062         raise ValueError("matrix equation has no solutions")
   1063 return X
   1064 

ValueError: matrix equation has no solutions

Here is the one with infinitely many solutions:

In [64]:
v = vector(QQ, [1, 1])
B.solve_right(v)
Out[64]:
(1, 0)

Note that when there are infinitely many solutions, it only gives one particular solution! To see if there are more, we need to inspect the kernel/nullspace:

In [65]:
B.right_kernel()
Out[65]:
Vector space of degree 2 and dimension 1 over Rational Field
Basis matrix:
[ 1 -1]

So, we see that there kernel is not trivial, and thus we have infinitely many solutions.

Let's save the particular solutions in s1:

In [66]:
s1 = B.solve_right(v)

Let's save the basis for the kernel in kernel_basis:

In [67]:
kernel_basis = B.right_kernel().basis()

kernel_basis
Out[67]:
[
(1, -1)
]

We can see the the size of the basis (i.e., the dimension) is $1$:

In [68]:
len(kernel_basis)  # len = length = number of elements
Out[68]:
1

We can now produce new solutions:

In [69]:
B * (s1 + 37/54 * kernel_basis[0]) == v
Out[69]:
True

Solving Systems with the Augmented Matrix¶

In [70]:
A = matrix(QQ, 3, 3, [2, 8, 4, 2, 5, 1, 4, 10, -1])
v = vector(QQ, [2, 5, 1])
In [71]:
B = A.augment(v)
B
Out[71]:
[ 2  8  4  2]
[ 2  5  1  5]
[ 4 10 -1  1]
In [72]:
B.echelon_form()
Out[72]:
[ 1  0  0 11]
[ 0  1  0 -4]
[ 0  0  1  3]
In [73]:
B.echelon_form().column(3)
Out[73]:
(11, -4, 3)

We can solve multiple systems at once! To solve $A \cdot \vec{x} = \vec{u}$, $A \cdot \vec{x} = \vec{v}$, $A \cdot \vec{x} = \vec{w}$, we can do:

In [74]:
A = matrix(QQ, 3, 3, [2, 8, 4, 2, 5, 1, 4, 10, -1])
u = vector(QQ, [1, 1, 0])
v = vector(QQ, [2, 5, 1])
w = vector(QQ, [0, -1, 3])

We augment the matrix by the three vectors:

In [75]:
B = A.augment(matrix([u, v, w]).T)  # note the .T to make the vectors into columns!
B
Out[75]:
[ 2  8  4  1  2  0]
[ 2  5  1  1  5 -1]
[ 4 10 -1  0  1  3]
In [76]:
C = B.echelon_form()
C
Out[76]:
[    1     0     0  11/6    11 -14/3]
[    0     1     0  -2/3    -4     2]
[    0     0     1   2/3     3  -5/3]

The solutions:

In [77]:
C.column(3), C.column(4), C.column(5)
Out[77]:
((11/6, -2/3, 2/3), (11, -4, 3), (-14/3, 2, -5/3))

We can also use this method to see when there are no solutions:

In [78]:
A = matrix(QQ, 2, 2, 4 * [1])
A
Out[78]:
[1 1]
[1 1]
In [79]:
v = vector(QQ, [1, 2])
w = vector(QQ, [1, 1])
B = A.augment(matrix([v, w]).T)
B
Out[79]:
[1 1 1 1]
[1 1 2 1]
In [80]:
B.echelon_form()
Out[80]:
[1 1 0 1]
[0 0 1 0]