Authors: Judith Ludwig and Florent Schaffhauser, Uni-Heidelberg, Summer Semester 2024
Matrices in Sage¶
Basic operations¶
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 determinant 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
Reduced row-echelon form¶
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
Exercises¶
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. $$
Invertible matrices¶
Computing the inverse¶
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
Row operations and elementary matrices¶
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 A 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?
Exercises¶
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$.
Projects¶
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.
Elementary matrices¶
In this project, you are asked to write a short program that lets you:
- Input a matrix $A$.
- Ask the user which row operation they want to apply to that matrix: this row operation should be encoded under the form $\lambda R_i + R_j \rightarrow R_j$ or $R_i \leftrightarrow R_j$ or $\lambda R_i \rightarrow R_i$.
- Output the elementary matrix corresponding to the row operation entered by the user and the resulting new matrix $A$.
PLU Decomposition¶
In this project, you are asked to write a program that lets you:
- Input a matrix $A$.
- Determine whether $A$ is invertible and, if so, compute its inverse and write it down as a product of elementary matrices.
This entails:
- Recalling why the inverse of an invertible matrix can be expressed as a product of elementary matrices.
- Looking into LU and PLU decompositions in Sage.
Cramer's rule¶
In this project, you are asked to write a program that lets you:
- Input a matrix $A$.
- Input a vector $Y$ such that the equality $AX=Y$ makes sense.
- Determine whether $A$ is invertible and, if so, compute the solution $X$ of the system $AX=Y$ using Cramer's rule.
This entails:
- Recalling the mathematical proof of Cramer's rule.
- In the program, substitute a given column of $A$ with $Y$ and compute the determinant of the matrix thus obtained.
Gauß's algorithm¶
This is the most programmation-oriented project. You are asked to write a program that lets you:
- Input a matrix.
- Put it in reduced row-echelon form.
- (If possible) Indicate the list of elementary operations that has been performed, along with the corresponding elementary matrices.