Calculating the Inverse of a Matrix

Okay, we've established that if a square matrix $A$ has an inverse $A^{-1}$, we can solve the system $Ax = b$ by calculating $x = A^{-1}b$. But how do we actually find this matrix $A^{-1}$? The method depends heavily on the size of the matrix and whether you're doing it by hand or using a computer.

The 2x2 Case: A Simple Formula

For a 2x2 matrix, there's a straightforward formula. Let's say our matrix is:

$$A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$$

Its inverse, $A^{-1}$, is given by:

$$A^{-1} = \frac{1}{ad - bc} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix}$$

Notice the term $ad - bc$ in the denominator. This is the determinant of the 2x2 matrix $A$, which we'll discuss more in the next section. For the inverse to exist, this determinant must not be zero. If $ad - bc = 0$, the matrix is singular (or non-invertible), and $A^{-1}$ does not exist.

The formula essentially involves swapping the diagonal elements ($a$ and $d$), negating the off-diagonal elements ($-b$ and $-c$), and dividing the entire resulting matrix by the determinant.

Example: Find the inverse of $A = \begin{bmatrix} 4 & 7 \\ 2 & 6 \end{bmatrix}$.

1. Calculate the determinant: $ad - bc = (4)(6) - (7)(2) = 24 - 14 = 10$. Since $10 \neq 0$, the inverse exists.
2. Apply the formula: $$A^{-1} = \frac{1}{10} \begin{bmatrix} 6 & -7 \\ -2 & 4 \end{bmatrix} = \begin{bmatrix} 6/10 & -7/10 \\ -2/10 & 4/10 \end{bmatrix} = \begin{bmatrix} 0.6 & -0.7 \\ -0.2 & 0.4 \end{bmatrix}$$

You can verify this by multiplying $A$ and $A^{-1}$ (in either order) to see if you get the identity matrix $I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$.

Larger Matrices: Gauss-Jordan Elimination

For matrices larger than 2x2 (or sometimes even 3x3), calculating the inverse by hand using formulas (like the adjugate method sometimes taught in pure math courses) becomes extremely tedious and prone to errors. The standard algorithmic approach, and the one typically implemented in software, is based on Gauss-Jordan elimination.

The idea is to augment the matrix $A$ with the identity matrix $I$ of the same size, creating $[A | I]$. Then, you apply elementary row operations (the same ones used in Gaussian elimination to solve systems) to transform the left side ($A$) into the identity matrix $I$. The magic is that the same sequence of row operations will transform the right side ($I$) into the inverse matrix $A^{-1}$. The final form will be $[I | A^{-1}]$.

Elementary Row Operations include:

1. Swapping two rows.
2. Multiplying a row by a non-zero scalar.
3. Adding a multiple of one row to another row.

While understanding this process is valuable conceptually, you will almost never perform it manually in a machine learning context. It's computationally intensive and best left to optimized libraries.

Using NumPy for Matrix Inversion

In practice, especially within machine learning workflows using Python, you'll rely on numerical libraries like NumPy to compute matrix inverses. NumPy's linear algebra module, numpy.linalg, provides the inv() function.

Here's how you'd compute the inverse of our previous 2x2 example matrix using NumPy:

import numpy as np

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

# Calculate the inverse
try:
    A_inv = np.linalg.inv(A)
    print("Matrix A:")
    print(A)
    print("\nInverse of A:")
    print(A_inv)

    # Verification (optional): Multiply A by A_inv
    # Should result in the identity matrix (within floating point precision)
    identity_check = np.dot(A, A_inv)
    print("\nVerification (A @ A_inv):")
    # Use np.allclose to handle potential tiny floating-point errors
    print(np.round(identity_check, decimals=5))
    print("\nIs close to identity:", np.allclose(identity_check, np.eye(2)))

except np.linalg.LinAlgError:
    print("\nMatrix A is singular and cannot be inverted.")

Output:

Matrix A:
[[4 7]
 [2 6]]

Inverse of A:
[[ 0.6 -0.7]
 [-0.2  0.4]]

Verification (A @ A_inv):
[[1. 0.]
 [0. 1.]]

Is close to identity: True

This matches our manual calculation.

For a larger matrix:

# Example of a 3x3 matrix
B = np.array([[1, 2, 3],
              [0, 1, 4],
              [5, 6, 0]])

# Calculate the inverse
try:
    B_inv = np.linalg.inv(B)
    print("\nMatrix B:")
    print(B)
    print("\nInverse of B:")
    print(B_inv)

    # Verification
    identity_check_B = np.dot(B, B_inv)
    print("\nVerification (B @ B_inv close to identity):", np.allclose(identity_check_B, np.eye(3)))

except np.linalg.LinAlgError:
    print("\nMatrix B is singular and cannot be inverted.")

Output:

Matrix B:
[[1 2 3]
 [0 1 4]
 [5 6 0]]

Inverse of B:
[[-24.  18.   5.]
 [ 20. -15.  -4.]
 [ -5.   4.   1.]]

Verification (B @ B_inv close to identity): True

NumPy handles the underlying Gauss-Jordan elimination (or more numerically stable variants) efficiently. If you pass a singular matrix (one whose determinant is zero) to np.linalg.inv(), it will raise a LinAlgError.

Before attempting to calculate an inverse, it's often wise to check if the matrix is invertible in the first place. This brings us to the determinant.