# Matrix Inversion
Matrix inversion is the process of finding a matrix \(A^{-1}\) such that
where \(\mathbb{1}\) is the identity matrix (a matrix of all zeros except with ones along the diagonal).
Consider
The inverse of \(A\) is
such that
Let us check this in Python using the np.linalg.inv function. Here is the code.
import numpy as np
# Define a 2x2 matrix
A = np.array([[4, 7],
[2, 6]])
# Compute the inverse
A_inv = np.linalg.inv(A)
print("Inverse of A:")
print(A_inv)
print("Check that multiplication gives the identitiy.")
print(A@A_inv)
Inverse of A:
[[ 0.6 -0.7]
[-0.2 0.4]]
Check that multiplication gives the identitiy.
[[ 1.00000000e+00 -1.11022302e-16]
[ 1.11022302e-16 1.00000000e+00]]
The code output is as expected. Note that the off-diagonal elements above are not exactly zero. This is due to round-off error. More on round-off error in a future lesson. Let us check a larger matrix. To do so, we will make use of the np.random.rand function to generate a matrix composed of random numbers from 0 to 1.
A=np.random.rand(5,5)
A_inv=np.linalg.inv(A)
A_inv@A
array([[ 1.00000000e+00, -7.05024015e-16, -5.43283954e-16,
-9.77184073e-16, -7.98461246e-16],
[ 3.81205529e-16, 1.00000000e+00, 2.32739549e-16,
4.70686067e-16, 1.84345679e-16],
[-2.59402170e-17, -3.40972115e-16, 1.00000000e+00,
-3.51699704e-16, -1.52644921e-16],
[ 9.19667391e-17, 5.21793986e-17, -5.17185710e-17,
1.00000000e+00, 1.99604795e-17],
[-1.28449376e-16, -1.50622332e-16, 7.77983740e-17,
-5.78610138e-17, 1.00000000e+00]])
Yes, the output is a bit messy. But if you inspect it, you will see that it is the identity matrix up to round off error.
3.5. Solving Systems of Equations#
Matrix inversion appears equations of the form
where \(A\) and \(\mathbf{b}\) are a known matrix and vector and \(\mathbf{x}\) is unknown.
The solution can be solved by multiplying both side by \(A^{-1}\)
A typical use case is solving systems of equations. For example consider the following
This can be written as
We can solve this system in two ways with Python.
Using the np.linalg.inv function.
Using the np.linalg.solve function.
Here is code for both methods.
A = np.array([[4, 7], [2, 6]])
b = np.array([1, 3])
# method 1
A_inv=np.linalg.inv(A)
x = A_inv@b
print(x)
[-1.5 1. ]
# method 2
x=np.linalg.solve(A, b)
print(x)
[-1.5 1. ]
The two methods are equivalent. But it is wise not to directly invert matrices as in method 1 and instead use methods similar to method 2. This topic will be discussed in more detail in latter lectures.
3.6. Should You Directly Invert a Matrix in Numerical Methods? (Hint: almost always No)#
The short answer is: No, you should generally avoid directly inverting a matrix when solving systems of equations.
In numerical computing, using np.linalg.inv(A) and then multiplying with b like x = A⁻¹b can lead to problems:
Problem |
Why it’s bad |
|---|---|
Numerically unstable |
Floating-point errors can result in inaccurate solutions |
Computationally costly |
Inverting an \(n \times n\) matrix takes ~(O(n^3)) operations |
Unnecessary |
You almost never need the explicit inverse to solve a system |
Instead of this:
# Not recommended
x = np.linalg.inv(A) @ b
Do this:
# Recommended
x = np.linalg.solve(A, b)
Before moving on, I must say that studing matrix methods will be a major componet of this course. Be sure are familiar with basics outlines in these notebooks. Don’t forget, that in most cases, it is best not to directly invert a matrix.
3.7. Inversion for \(3\times3\) Matrices#
Here I will present code to invert a \(3\times 3\) matrix using the determinate and the adjugate. This is not the type of coding you would expect to do when doing real scientific computing, but it still helps the student gain familiarity with indexing and array manipulations.
Let \(A\) be a \(3 \times 3\) matrix given by:
The inverse of \(A\), denoted by \(A^{-1}\), is given by:
where \(\det(A)\) is the determinant of \(A\) and \(\text{adj}(A)\) is the adjugate of \(A\).
Step 1: Calculate the Determinant of A (\(\det(A)\)) The determinant of a \(3 \times 3\) matrix is calculated as:
If \(\det(A) = 0\), the matrix \(A\) is singular and has no inverse.
Step 2: Find the Matrix of Minors (\(M\)) The matrix of minors \(M\) is formed by taking the determinant of the \(2 \times 2\) submatrix obtained by removing the \(i\)-th row and \(j\)-th column for each element of \(A\):
Step 3: Find the Matrix of Cofactors (\(C\)) The matrix of cofactors \(C\) is obtained by applying a checkerboard pattern of signs to the matrix of minors \(M\):
where \(\odot\) denotes the element-wise product.
Step 4: Find the Adjugate (Adjoint) of A (\(\text{adj}(A)\)) The adjugate of \(A\) is the transpose of the cofactor matrix \(C\):
Step 5: Calculate the Inverse \(A^{-1}\) Finally, the inverse of \(A\) is obtained by multiplying the adjugate of \(A\) by \(\frac{1}{\det(A)}\):
def determinant_3x3(matrix):
"""
Computes the determinant of a 3x3 NumPy matrix.
Args:
matrix (numpy.ndarray): A 3x3 NumPy matrix.
Returns:
float: The determinant of the input matrix.
Returns None if the input is not a 3x3 matrix.
"""
if matrix.shape != (3, 3):
print("Error: Input matrix must be 3x3.")
return None
return (matrix[0, 0] * (matrix[1, 1] * matrix[2, 2] - matrix[1, 2] * matrix[2, 1]) -
matrix[0, 1] * (matrix[1, 0] * matrix[2, 2] - matrix[1, 2] * matrix[2, 0]) +
matrix[0, 2] * (matrix[1, 0] * matrix[2, 1] - matrix[1, 1] * matrix[2, 0]))
def adjugate_3x3(matrix):
"""
Computes the adjugate (adjoint) of a 3x3 NumPy matrix.
Args:
matrix (numpy.ndarray): A 3x3 NumPy matrix.
Returns:
numpy.ndarray: The adjugate of the input matrix.
Returns None if the input is not a 3x3 matrix.
"""
if matrix.shape != (3, 3):
print("Error: Input matrix must be 3x3.")
return None
# Calculate the matrix of minors
minors = np.zeros((3, 3))
for i in range(3):
for j in range(3):
submatrix = np.delete(np.delete(matrix, i, axis=0), j, axis=1)
minors[i, j] = np.linalg.det(submatrix)
# Calculate the matrix of cofactors
cofactors = np.zeros((3, 3))
for i in range(3):
for j in range(3):
cofactors[i, j] = minors[i, j] * ((-1)**(i + j))
# The adjugate is the transpose of the cofactor matrix
adjugate = cofactors.T
return adjugate
# test
C=np.random.rand(3,3)
detC=determinant_3x3(C)
adj=adjugate_3x3(C)
C_inv=1/detC*adj
print("Inverse of C", C_inv)
print("Check that C^{-1}C, results in identity matrix", C_inv@C)
Inverse of C [[ 1.50329814 1.51865638 -2.47739043]
[ 0.101985 -1.24087009 1.8631654 ]
[-1.40636296 1.16184191 0.78517277]]
Check that C^{-1}C, results in identity matrix [[ 1.00000000e+00 4.58218250e-17 1.53679958e-16]
[-1.69074210e-16 1.00000000e+00 -1.88051858e-16]
[ 2.04348762e-16 4.27450856e-16 1.00000000e+00]]