Inverting matrices
The first topic we discussed this week is inverting matrices. This works exactly as you would expect:
M = matrix([[1,2,3],[1,1,1],[3,1,4]])
M.inverse()
As a sanity check:
M * M.inverse()
(Note that the .inverse() comes first in the order of operations, as it should. If we didn’t know Python would do this way, we would add parenthesis — M*(M.inverse()).)
You can also invert a matrix by raising it to the -1’st power:
M^(-1)
Of course, if the matrix is not invertible Sage raises and exception (i.e., gives an error):
M = matrix([[1,2],[4,8]])
M.inverse()
N = matrix([[1,2,3],[4,5,6]])
N.inverse()
In fact, if you look at the last line of the error message, it tells you briefly why the matrix was not invertible (“singular” meaning square but not invertible in the first example, not square in the second example).
LU factorizations
Sage also has built in code for computing LU factorizations. As the textbook notes:
- Not every matrix has an LU factorization, unless one allows L to be permuted lower triangular, i.e., unless one allows pivoting.
- Even if you only consider matrices which have LU factorizations, roundoff error causes problems unless you use partial pivoting when computing the factorization.
The Sage LU factorization code has several options, depending on how much pivoting you want to allow. The function that returns an LU factorization of M is M.LU(). To see its options, type “M.LU?” and press shift-return. (There’s a lot of documentation there, so don’t print it all out.) The option that matches class most closely is M.LU(pivot=”nonzero”), which means “only do row swapping when absolutely necessary”.
I find the output of M.LU illegible unless we turn on pretty printing. So, try:
%typeset_mode True
M.LU(pivot='nonzero')
Hmm…why did we get three matrices? The first matrix is a permutation matrix: it says how you need to permute (swap) the rows in the second matrix L so that LU is right. If no pivoting was needed (i.e., there is an honest LU factorization) then the first matrix will be the identity matrix (like here). By contrast, if we try:
N = matrix([[1,2,3],[2,4,7],[1,1,1]])
N.LU(pivot='nonzero')
Let’s give names to the matrices in the LU factorization:
(P,L,U) = M.LU(pivot='nonzero')
P
L
U
(P stands for permutation, though you could of course use any names you like.) Next, let’s make sure that LU really is the right thing:
L*U
L*U == M
For N, L*U will not be N:
(P2, L2, U2) = N.LU(pivot='nonzero')
(P2, L2, U2)
L2 * U2
L2 * U2 == N
P2 * L2
P2 * L2 * U2
P2 * L2 * U2 == N
It works…