The Inverse Matrix Isn't as Great as You Think

Today, let’s talk about a familiar term: the inverse matrix. We’ll explore why inverse matrices aren’t as great as you might think, and why we should avoid using them in practical computation.

Solving Linear Equations

As mentioned in the previous article, we can directly find the inverse of matrix A and then calculate $A^{-1}b$ to solve for x.

That sounds effortless. However, in practical applications, solving linear equations by directly computing the inverse matrix is relatively rare. In real-world scenarios, matrix elements aren’t all neat integers; they are usually filled with various floating-point numbers, and the dimensions are much larger. They aren’t like college exam problems, which are at most 3x3 matrices with numbers carefully crafted to avoid decimal calculations.

The Inverse Matrix

There are two classic approaches to finding an inverse matrix:

• Using the adjugate matrix (adj matrix) and determinant to solve for the inverse matrix

• Using Gauss-Jordan elimination to solve for the inverse matrix

Intuitively, finding the inverse matrix using the adjugate matrix and determinant is far more troublesome than using Gauss-Jordan elimination. In fact, adjugate matrices and determinants are indeed ill-suited for computer computation.

Since the formula contains

when the determinant approaches zero, this value becomes extremely large, leading to numerical errors or even overflow. In addition, computing the adjugate matrix requires breaking matrix A down into smaller submatrices and calculating their determinants, which has a complexity of O(n!). Therefore, even though the formula looks simple, it is practically terrible to compute.

A more common way to find the inverse matrix is using Gauss-Jordan elimination. We take the augmented matrix

and reduce it to an upper triangular matrix using row operations, then eliminate all off-diagonal elements to zero. At this point, the augmented matrix becomes

Here, B is the inverse matrix of A. If written in Python, it would look like this (assuming the matrix is invertible):

import numpy as np

def gauss_jordan(m):
    n = len(m)
    identity = np.eye(n)

		# 增廣矩陣
    augmented = np.hstack((m, identity))

    for i in range(n):
        # checking if diagonal element is zero
        if augmented[i][i] == 0:
            # swapping with row below
            for j in range(i+1, n):
                if augmented[j][i] != 0:
                    augmented[[i, j]] = augmented[[j, i]]
                    break
        augmented[i] = augmented[i] / augmented[i, i]
        for j in range(0, n):
            if i != j:
                augmented[j] -= augmented[j, i] * augmented[i]

    # return the right side of the augmented matrix
    return augmented[:, n:]


m = np.array([[1., 2.],
              [3., 4.]])
inv_m = gauss_jordan(m)

NumPy provides the np.linalg.inv function; I wrote a simple version here for demonstration purposes. Looking at its time complexity, we can see that it is O(n^3).

LU Decomposition

Solving linear equations by computing the inverse matrix on a computer is actually not very common. It’s just that under the hand-calculation-oriented mindset from college, one might wonder: “Why go through all the trouble of matrix decomposition when we could just use the inverse matrix directly?”

LU decomposition breaks down matrix A into two matrices, L and U—a lower triangular matrix and an upper triangular matrix, respectively. Specifically, L is a matrix with diagonal elements of 1, elements below the diagonal may be non-zero, and elements above the diagonal are all 0; U is a matrix where elements below the diagonal are all 0.

After computing LU, we can then solve via:

An implementation using NumPy looks roughly like this:

import numpy as np

def lu_decomposition(A):
    n = len(A)
    L = np.zeros((n, n))
    U = np.zeros((n, n))

    for i in range(n):
        L[i][i] = 1  # Diagonal elements of L are 1

        for j in range(i, n):
            # Compute the upper triangular matrix U
            summation = sum(L[i][k] * U[k][j] for k in range(i))
            U[i][j] = A[i][j] - summation


        for j in range(i + 1, n):
            # Compute the lower triangular matrix L
            summation = sum(L[j][k] * U[k][i] for k in range(i))
            L[j][i] = (A[j][i] - summation) / U[i][i]

    return L, U

def lu_solve(L, U, b):
    n = len(L)
    y = np.zeros(n)
    x = np.zeros(n)

    # Solving L * y = b
    for i in range(n):
        summation = sum(L[i][j] * y[j] for j in range(i))
        y[i] = (b[i] - summation) / L[i][i]

    # Solving U * x = y
    for i in range(n - 1, -1, -1):
        summation = sum(U[i][j] * x[j] for j in range(i + 1, n))
        x[i] = (y[i] - summation) / U[i][i]

    return x

A = np.array([[1, 1, 1], [2, 3, 1], [7,4,2]])
b = np.array([3, 6, 13])

L, U = lu_decomposition(A)
x = lu_solve(L, U, b)

Compared to Gauss-Jordan elimination, the computational cost of LU decomposition is roughly:

Compared to the O(n^3) of Gauss-Jordan elimination, this is about 3 times faster.

The non-intuitive part here is: LU decomposition adds the steps of solving L * y and U * x, which seems more troublesome, yet it is still better than using an inverse matrix? Due to the properties of triangular matrices, solving Ly and Ux only requires roughly O(n^2) operations.

Conclusion

In terms of raw computational complexity, the difference compared to finding the inverse with Gauss-Jordan elimination isn’t massive. However, in practical applications, the systems of equations we need to solve often involve sparse matrices—that is, matrices where 0 makes up the vast majority of elements.

For this type of problem, if you compute the inverse matrix directly, not only will the resulting inverse matrix be messy, but most of its elements will no longer be zero. On the other hand, if you use LU decomposition, most of the zero elements in matrix A can be preserved in L and U, making computations much simpler.

Afterword

Many solutions that seem intuitive at first glance are not necessarily the optimal ones. And the way to find the optimal solution often comes from “looking at the world from a different angle.” For instance, techniques like the Fast Fourier Transform (FFT) and Discrete Cosine Transform (DCT), which we’ll discuss later—while they may seem a bit roundabout—are ingenious strategies for solving problems.