Calculating Inverses with NumPy

Now that we understand the concept of the matrix inverse ($A^{-1}$) and its theoretical role in solving systems of linear equations of the form $Ax = b$, let's look at how to actually compute the inverse using NumPy. While directly calculating the inverse isn't always the preferred method for solving systems in practice (as we'll see in the next section), knowing how to compute it is a fundamental skill.

The numpy.linalg Module

NumPy, the cornerstone library for numerical computing in Python, includes a submodule dedicated to linear algebra operations: numpy.linalg. This module contains functions for matrix decomposition, eigenvalue calculation, solving linear systems, and, importantly for us right now, calculating matrix inverses.

Calculating the Inverse with inv()

To compute the inverse of a square matrix, we use the numpy.linalg.inv() function. It takes a square matrix (represented as a NumPy array) as input and returns its inverse, also as a NumPy array.

Remember, only square matrices can have inverses, and even then, only if they are non-singular (invertible).

Let's try an example. Consider the following 2x2 matrix $A$:

$$A = \begin{pmatrix} 4 & 7 \\ 2 & 6 \end{pmatrix}$$

We can represent this in NumPy and calculate its inverse:

import numpy as np

# Define matrix A
A = np.array([[4, 7],
              [2, 6]])

print("Matrix A:")
print(A)

# Calculate the inverse of A
try:
    A_inv = np.linalg.inv(A)
    print("\nInverse of A (A_inv):")
    print(A_inv)
except np.linalg.LinAlgError as e:
    print(f"\nCould not compute inverse: {e}")

Executing this code will output the original matrix $A$ and its calculated inverse $A_{inv}$.

Verifying the Inverse

A defining property of the matrix inverse $A^{-1}$ is that when multiplied by the original matrix $A$, it yields the identity matrix $I$. That is:

$$A A^{-1} = A^{-1} A = I$$

We can verify this using NumPy's matrix multiplication capabilities. Recall that the @ operator performs matrix multiplication (or you can use np.dot()):

# Verify A * A_inv
identity_check_1 = A @ A_inv
print("\nVerification (A @ A_inv):")
print(identity_check_1)

# Verify A_inv * A
identity_check_2 = A_inv @ A
print("\nVerification (A_inv @ A):")
print(identity_check_2)

# Create the expected identity matrix
I = np.identity(A.shape[0]) # A.shape[0] gives the number of rows (which is 2 here)
print("\nIdentity matrix I:")
print(I)

You should see that both identity_check_1 and identity_check_2 produce a matrix very close to the 2x2 identity matrix:

$$I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$$

A Note on Floating-Point Precision

You might notice that the resulting matrices aren't exactly the identity matrix. Instead of perfect zeros, you might see very small numbers like 2.22044605e-16 (which is $2.22 \times 10^{-16}$). This is a normal consequence of how computers handle calculations with non-integer numbers, known as floating-point arithmetic. These tiny differences are usually negligible.

If you need to programmatically check if two matrices are equal within a certain tolerance, NumPy provides the np.allclose() function:

# Check if A @ A_inv is close to the identity matrix I
are_close = np.allclose(A @ A_inv, np.identity(A.shape[0]))
print(f"\nIs A @ A_inv numerically close to I? {are_close}") # Output should be True

What About Non-Invertible Matrices?

What happens if we try to compute the inverse of a matrix that doesn't have one (a singular matrix)? Let's try with a matrix where one column is a multiple of another, making it singular:

$$B = \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix}$$

# Define a singular matrix B
B = np.array([[1, 2],
              [2, 4]])

print("\nMatrix B (Singular):")
print(B)

# Attempt to calculate the inverse of B
try:
    B_inv = np.linalg.inv(B)
    print("\nInverse of B:")
    print(B_inv)
except np.linalg.LinAlgError as e:
    print(f"\nCould not compute inverse of B: {e}")

When you run this code, NumPy will detect that the matrix is singular and raise a LinAlgError: Singular matrix. This is NumPy's way of telling you that the inverse does not exist for the given matrix B. This behavior correctly reflects the mathematical properties discussed in the "Conditions for Invertibility" section.

Calculating the inverse is a fundamental operation, and np.linalg.inv() provides a straightforward way to do it. However, keep in mind that for solving systems of equations $Ax=b$, using the inverse directly ($x = A^{-1}b$) can sometimes be less numerically stable and efficient than using dedicated solvers, which we will explore next.