Author: Florent Schaffhauser, Uni-Heidelberg, Summer Semester 2023
When we define a matrix in Sage, we can specify the ring or field in which we take the entries.
Let us for instance consider the matrix
$$ \begin{pmatrix} 2 & 4 & 6 \\ 4 & 5 & 6 \\ 3 & 1 & 2 \end{pmatrix} $$and declare it first as a matrix $A$ with entries in $\mathbb{Q}$, then as a matrix $B$ with entries the field with seven elements $\mathbb{F}_7$.
A = matrix( QQ, [[2,4,6],[4,5,6],[3,1,2]] )
show(A)
# The command 'type()' shows what kind of object the argument is
type(A)
<class 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'>
B = matrix( GF(7), [[2,4,6],[4,5,6],[3,1,2]])
show(B)
# Matrices with entries in a finite field are of the "mod n" class
type(B)
<class 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
If we declare our matrix as a matrix with entries in the field with $3$ elements, Sage will automatically reduce the entries modulo $3$:
C = matrix( GF(3), [[2,4,6],[4,5,6],[3,1,2]])
show(C)
# The type of C in this case is the same as the type of B
type(C)
<class 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
As we know, depending on the ring in which the entries are taken to lie, the properties of a "same" matrix may be quite different. For instance, here, $B$ is invertible in $\mathrm{Mat}( 3 \times 3; \mathbb{F}_7 )$, but $C$ is not invertible in $\mathrm{Mat}( 3 \times 3; \mathbb{F}_3)$.
# The determinant function computes the deerminant of square matrix...
A.determinant()
-18
# Recall that -18 mod 7 is 3
B.determinant()
3
# And -18 mod 3 is 0
C.determinant()
0
Also, the inverse of $A$ in $\mathrm{Mat}( 3 \times 3; \mathbb{Q} )$ looks very different from the inverse of $B$ in $\mathrm{Mat}( 3 \times 3; \mathbb{F}_7 )$.
show(A.inverse())
show(B.inverse())
In Sage, we can perform usual matrix operations in a fairly intuitive manner:
show(2*A^(-1)+A^2-3*A)
show(2*B^(-1)+B^2-3*B)
Note that the command "A^(-1)" returns the inverse of $A$.
# We can test in the following way whether an assertion is True or False
A^(-1) == A.inverse()
True
3*B == B^2
False
In order to solve homogeneous linear systems $AX=0$, one reduces $A$ to its reduced row-echelon form (or RREF), which Sage can do.
# Recall our matrix A (with rational entries)
show(A)
# The function echelon_form() returns the RREF of a matrix
show(A.echelon_form())
In this case, the reduced row-echelon form of $A$ is the identity, which shows that $\ker A=\{0\}$. Note that the computation we have done has not affected the variable $A$.
show(A)
In order to solve general linear systems $AX=Y$, we can perform Gaussian elimination on the augmented matrix $(A|Y)$.
For instance, if we consider the linear system $$ \begin{pmatrix} 2 & 4 & 6 \\ 4 & 5 & 6 \\ 3 & 1 & 2 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 18 \\ 24 \\ 4 \end{pmatrix} , $$ we can proceed as follows.
# The function A.augment(Y) is used to define the augmented matrix M = (A|Y)
# We can put a vertical bar to separate A from Y in the resulting matrix M
A = matrix( QQ, [ [2,4,6], [4,5,6], [3,1,2] ] )
# Y = matrix( QQ, [ [18], [24], [4] ] )
Y = column_matrix( QQ, [[ 18, 24, 4 ]])
M = A.augment(Y, subdivide = "True")
# We display the augmented matrix (A|Y)
# Note the vertical bar separating A from Y
show(M)
# The RREF is computed in the exact same way as before
show(M.echelon_form())
The rightmost column of the matrix thus obtained is the solution of the system $AX=Y$, as we now check.
# We first store the reduced row-echelon form of M in a variable N
N = M.echelon_form()
# Then we extract the entries contained in Rows 0,1 and 2 and Column 3
# Note that Row 0 is the first row (in maths, we usually call it Row 1)
X = N.matrix_from_rows_and_columns([0,1,2], [3])
show(X)
# We check that the expected equality indeed holds
A*X == Y
True
# Note that the following also works and displays the solution as a vector
X1 = N.column(3)
show(X1)
type(X1)
<class 'sage.modules.vector_rational_dense.Vector_rational_dense'>
# Since Y was declared as a column matrix, we must make it a vector in order to compare it to A*X1
A*X1 == Y.column(0)
True
Exercise 1. Why does it not work if we try to compute $B+C$ with $B$ and $C$ as above?
Exercise 2. Declare a $2 \times 2$ matrix $A$ with entries in $\mathbb{C}$ and compute its determinant and its square $A^2$. Check the computations manually.
Exercise 3. Use Sage to reduce the following linear system to its reduced row-echelon form and solve manually the equivalent system thus obtained.
$$ \left\{ \begin{array}{rcl} 2x + 4y + 6z & = & 18 \\ 4x + 5y + 6z & = & 24 \\ 2x + 7y + 12z & = & 30 \end{array} \right. $$Given a square matrix $A$ with entries in a ring $R$, we can check whether it is invertible by computing its determinant: if $\det(A)$ is invertible in $R$, then $A$ is invertible.
Zmod(36)
Ring of integers modulo 36
Zmod(24).unit_group_order()
8
A = matrix( Zmod(36), [ [1,-1], [1,42] ])
show(A)
show(det(A))
# Check whether a determinant is invertible in a given ring (here Z/6Z)
det(A).is_unit()
True
# If det(A) is invertible, one can compute its multiplicative inverse easily
det(A)^(-1)
31
If the matrix $A$ is invertible, we can compute its inverse directly, using the inverse() function.
show(A.inverse())
However, if the entries of $A$ are taken to lie in a field, then there is another method, based on the reduction to the RREF of $A$, which tells us if $A$ is invertible and, if so, what $A^{-1}$ is.
If $A$ is of size $n \times n$, we simply reduce the augmented matrix $(A|I_n)$ to its RREF $(H|B)$. Then $A$ is invertible if and only if $H=I_n$ and, in this case $A^{-1}=B$.
Observation: We need the entries of $A$ to lie in a field in order to apply Gaussian elimination. Note that $H$ is the RREF of $A$ in the computation above.
# Let us work in the field with 2^4 = 32 elements
# You can also declare the field GF(2^4) as GF(32) or GF(2,4)
# Whatever you choose, you should use the same convention consistenly throughout the file
A = matrix( GF(2^4), [ [1,-1], [1,42] ])
show(A)
# We construct the identity matrix I and the augmented matrix (A|I)
# For some reason, we cannot use "subdivide" when working with matrices with entries in GF(2^4)
I2 = matrix.identity(GF(2^4),2)
M = A.augment(I2)
show(M)
# We reduce M to its RREF
# The left 2x2 matrix is the identity, so M is invertible
show(M.echelon_form())
# The rightmost 2x2 matrix is the inverse of M and we extract it (directly from M.echelon_form())
B = (M.echelon_form()).matrix_from_rows_and_columns([0,1], [2,3])
show(B)
# We check that B is indeed the inverse of A
B == A.inverse()
True
In Sage, it is easy to manipulate matrices and perform operations on them, for instance on rows.
# We declare a 3x3 matrix (in an alternate way) and perform row operations on it
A = matrix(QQ, 3, 3, [2, 4, 6, 4, 5, 6, 3, 1, -2])
show(A)
# Here is how to modify a coefficient of A (storing it in a new variable first)
B = copy(A)
B[2,0] = 0
B[2,1] = 3 * B[2,1]
show(B)
# Let us start the reduction to row-echelon form of A (storing it in a new variable first)
# We proceed as in Gauß's algorithm
N = copy(A)
# The first operation replaces the first row of N with 1/2 times that row
N[0,:] = 1/2 * N[0,:]
# This replaces the second row with the second row minus 4 times the first row
N[1,:] = N[1,:] - 4 * N[0,:]
# What do you think this does?
N[2,:] = N[2,:] - 3 * N[0,:]
show(N)
# We go on like this
N[1,:] = -1/3 * N[1,:]
N[2,:] = N[2,:] + 5*N[1,:]
N[2:] = -N[2:]
show(N)
As we know, performing row operations on a $m \times n$ matrix is the same as multiplying to the left by an appropriate invertible $m \times m$ matrix. Likewise, performing column operations on a matrix is the same as multiplying to the right by an appropriate invertible $n \times n$ matrix.
In Sage, there is a command to define the elementary matrix corresponding to an elementary operation (e.g. swapping two rows or multiplying a column by a scalar). Unfortunately, the syntax is not very clear, due to Sage numbering method.
# The 4x4 elementary matrix which multiplies row 1 by -1 and adds it to row 4 of a 4xn matrix
E = elementary_matrix( QQ, 4, row1 = 3, row2 = 0, scale= -1 )
show(E)
# Swapping columns 1 and 2 of an mx3 matrix
S = elementary_matrix( QQ, 3, col1 = 0, col2 = 1 )
show(S)
# Scaling the second row of a 3xn matrix
M = elementary_matrix(QQ, 3, row1 = 1, scale = 42)
show(M)
# How would you scale the third column of an mx4 matrix by a factor of -33?
Exercise 1. Retake the matrix $N$ above and put it in RREF by writing down the row operations in Sage. Compare your result with A.echelon_form().
Exercise 2. How would you write down elementary operations on the columns of a matrix $A$?
Exercise 3. Consider the matrix $$ A = \begin{pmatrix} 2 & 4 & 6 \\ 4 & 5 & 6 \\ 3 & 1 & -2 \end{pmatrix} \in \mathrm{Mat}(3\times 3; \mathbb{Q})$$ Define an augmented matrix $M=(A|I_3)$ and reduce it to its RREF, extract the rightmost $3 \times 3$ matrix and check that it is an inverse for $A$.
Choose one the projects below and start working on it. Depending on your degree of proficiency with Sage or Python, you might need to read the documentation to make progress. You may also ask me for guidance and suggestions.
The goal is to present your findings in a short oral presentation later on in the seminar.
In this project, you are asked to write a short program that lets you:
In this project, you are asked to write a program that lets you:
This entails:
In this project, you are asked to write a program that lets you:
This entails:
This is the most programmation-oriented project. You are asked to write a program that lets you:
In this project, you are asked to rewrite this notebook in Python or in R (you may choose whichever you prefer). The commands and the syntax will be different, so you might need to read documentation pages, in particular in order to load the appropriate mathematical packages.
If you manage to solve this quickly, then try to implement either the first or the second project above in Python or in R.