Python

How to Calculate Matrix Inverse in Python

Jake @Scicoding

Feb 20, 2023 • 7 min read

Matrix inversion is one of the fundamental operations in linear algebra. It's analogous to division in ordinary algebra. Thus, making it crucial, for example, when solving many linear systems.

All essential Python linear algebra libraries support matrix inversion. We will have a look two most popular ones: NumPy and SciPy.

Matrix inversion plays such a fundamental role when working with linear systems that it's all so good to have an idea of how to implement some basic matrix inversion algorithms yourself. For example, when you are working with embedded or some more exotic platforms, where you don't have external libraries at hand. To this end, we also look at a simple implementation of the Gaussian elimination for matrix inversion.

Inverse Matrix Fundamentals

The definition of inverse matrices is fairly straightforward.

We define the inverse matrix \(\mathbf{X}\) of \(n \times n\) square matrix \(\mathbf{A}\) such that \[
\mathbf{XA} = \mathbf{AX} = \mathbf{I},
\] where \(\mathbf{I}\) is the identity matrix. If \(\mathbf{X}\) exists, we call \(\mathbf{A}\) invertible or non-singular and \(\mathbf{X}\) inverse of \(\mathbf {A}\).

While, in theory, the matrix inversion is pretty simple, there are a lot of implications for invertible matrices. Just to name of few

\(\mathbf{A}\) is full rank.
The rows and columns of \(\mathbf{A}\) are linearly independent.
\(\mathbf{A}\) has non-zero determinant.

It's easy to see why invertible matrices play such an important in linear algebra. Not to bore you too much with the fundamentals. Let's get to the topic. How to utilize matrix inversion in Python.

Test data for the matrix inversion

Before jumping into the actual methods for computing the inverse matrices, let's generate a standard 5x5 test matrix, which we use to calculate and verify the results.

import random

random.seed(0) # Set the seed for the RNG

N = 5

test_matrix = [[random.randint(-10, 10) for i in range(N)] for j in range(N)]

This will generate our 5x5 test matrix

[[2, 3, -9, -2, 6],
 [5, 2, -1, 5, 1],
 [8, -4, 6, -6, -1],
 [-6, -7, 9, -2, 7],
 [9, -6, -1, -7, -8]]

Calculating matrix inverse using NumPy

NumPy provides very good support for elementary linear algebra. Interface for the matrix inversion is provided via numpy.linalg.inv. Using this method is very simple. You pass the matrix you want to invert and it returns the inverse. There are no other parameters.

If the given matrix is not square or the matrix inversion fails for numerical reasons the method raises LinAlgError. It should be noted, that for larger matrices the matrix inversion may be somewhat sensitive, and badly scaled matrix can fail due to numerical instability.

import numpy as np

inverse = np.linalg.inv(test_matrix)

print(inverse @ test_matrix)

Calculating matrix inverse using NumPy

If all went well this will print out the identity matrix, which is sufficient for the validation of the results.

Calculating matrix inversion using SciPy

If you're using some more advanced computational science and engineering methods from NumPy that are not readily available in NumPy. You may not want to go back to NumPy just for matrix inversion. Instead, you can use scipy.linalg.inv to calculate matrix inversion.

The interface for scipy.linalg.inv is exactly the same as for the NumPy counterpart. You just pass the matrix to invert and get the inverted matrix as a return

import scipy.linalg

inverse = scipy.linalg.inv(test_matrix)

print(inverse @ test_matrix)

Calculation matrix inverse using SciPy

Again, what this should return is the identity matrix.

Custom algorithm for matrix inversion

Calculating matrix inversion without NumPy or SciPy can become handy in various situations. We will implement a simple Python-only algorithm for matrix inversion using Gaussian elimination.

Gaussian elimination

Finding matrix inverse using Gaussian elimination is really a very simple algorithm.

What we do is we take an input matrix \(\mathbf{A}\), and formulate an augmented matrix \[ \left[ \begin{array}{c | c} \mathbf{A} & \mathbf{I} \end{array} \right]. \]

Then, we start modifying both sides of the augmented matrix with the same operation. We are allowed to do the following:

Scale a row
Add scaled versions of rows together.

With these two operations, our target is to transform \(\mathbf{A}\) into an identity matrix, and what remains on the right side will be the inverse of \(\mathbf{A}\) that is \(\mathbf{A}^{-1}\).

We will accomplish this in three simple steps:

Transform \(\mathbf{A}\) into an upper triangular matrix.
Transform the upper triangular matrix into a diagonal matrix.
Scale the diagonal matrix into the identity matrix.

Example

We will begin with a simple step-by-step example of how to solve matrix inverse using Gaussian elimination. Our example matrix is
\[
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
2 & 0 & 2 & 0 & 1 & 0 \\
0 & 2 & 0 & 0 & 0 & 1
\end{array}
\right].
\]

Step 1: transform to an upper triangular matrix
We begin, by transforming the left side of the augmented matrix into an upper triangular matrix.

\[
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
2 & 0 & 2 & 0 & 1 & 0 \\
0 & 2 & 0 & 0 & 0 & 1
\end{array}
\right]
\Rightarrow
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
\color{red}{2} & 0 & 2 & 0 & 1 & 0 \\
0 & 2 & 0 & 0 & 0 & 1
\end{array}
\right]
\Rightarrow
\]\[
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
0 & -2 & 2 & -1 & 1 & 0 \\
0 & 2 & 0 & 0 & 0 & 1
\end{array}
\right]
\Rightarrow
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
0 & -2 & 2 & -1 & 1 & 0 \\
0 & \color{red}{2} & 0 & 0 & 0 & 1
\end{array}
\right]
\Rightarrow
\] \[
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
\color{green}{0} & -2 & 2 & -1 & 1 & 0 \\
\color{green}{0} & \color{green}{0} & 2 & -1 & 1 & 1
\end{array}
\right]
\]

Step 2: transform to a diagonal matrix

The next step is to transform the upper diagonal matrix into a diagonal matrix, eliminating the upper triangular non-zero elements.

\[
\left[
\begin{array}{ccc|ccc}
2 & 2 & 0 & 1 & 0 & 0 \\
0 & -2 & \color{red}{2} & -1 & 1 & 0 \\
0 & 0 & 2 & -1 & 1 & 1
\end{array}
\right] \Rightarrow
\left[
\begin{array}{ccc|ccc}
2 & \color{red}{2} & 0 & 1 & 0 & 0 \\
0 & -2 & 0 & 0 & 0 & -1 \\
0 & 0 & 2 & -1 & 1 & 1
\end{array}
\right] \Rightarrow
\] \[
\left[
\begin{array}{ccc|ccc}
2 & \color{green}{0} & \color{green}{0} & 1 & 0 & -1 \\
0 & -2 & \color{green}{0} & 0 & 0 & -1 \\
0 & 0 & 2 & -1 & 1 & 1
\end{array}
\right]
\]

Step 3: scale to identity

Finally, we will scale the diagonal matrix into an identity matrix by going through the diagonal elements one by one.

\[
\left[
\begin{array}{ccc|ccc}
\color{green}{2} & 0 & 0 & 1 & 0 & -1 \\
0 & \color{green}{-2} & 0 & 0 & 0 & -1 \\
0 & 0 & \color{green}{2} & -1 & 1 & 1
\end{array}
\right] \Rightarrow
\left[
\begin{array}{ccc|ccc}
\color{green}{1} & 0 & 0 & {1 \over 2} & 0 & -{1 \over 2} \\
0 & \color{green}{1} & 0 & 0 & 0 & {1 \over 2} \\
0 & 0 & \color{green}{1} & -{1 \over 2} & {1 \over 2} & {1 \over 2}
\end{array}
\right]
\]

Now, what we have on the right side of the augmented matrix is the inverse of our example matrix.

Implementation in Python

In our implementation, we keep track of both sides of the augmented matrix separately. However, the steps that we take are exactly the same as in our example

Transform the given input matrix into an upper triangular, and do the same steps for the identity (right side).
Transform the upper triangular matrix into the diagonal matrix, and apply the same steps to the (right side).
Scale each row of the diagonal matrix into identity. Do the same scaling for the right side.

Here is the sample code for matrix inversion without NumPy or Scipy

import copy

def matrix_inverse_gausssian_elimination(A):
    # Initialize target inverse matrix to identity
    inverse_matrix = [[0] * N for n in range(N)]

    for n in range(N):
        inverse_matrix[n][n] = 1.0

    # Make copy of the test matrix for modifications
    test_matrix_copy = copy.deepcopy(A)

    # Make the test matrix upper triangular
    for i in range(N): # Row
        # The diagonal element
        coeff = test_matrix_copy[i][i] + 0.0 

        for j in range(i+1, N):
            # How should we add ith row to jth, to eliminate the ith element
            mult = -test_matrix_copy[j][i] / coeff

            for n in range(N):
                test_matrix_copy[j][n] += mult * test_matrix_copy[i][n] 
                inverse_matrix[j][n] += mult * inverse_matrix[i][n] 

    # Diagonalize (same as triangulation but in the opposite direction)
    for i in range(N-1, -1, -1): # Row
        coeff = test_matrix_copy[i][i] + 0.0

        for j in range(i-1, -1, -1):
            mult = -test_matrix_copy[j][i] / coeff

            for n in range(N):
                test_matrix_copy[j][n] += mult * test_matrix_copy[i][n] 
                inverse_matrix[j][n] += mult * inverse_matrix[i][n] 

    # Scale to identity
    for i in range(N): 
        # Take the scaling coefficient from the diagonal
        scaler = test_matrix_copy[i][i]

        # Scale to identity
        test_matrix_copy[i][i] /= scaler

        # Scale the ith row of inverse with the corresponding scaler
        for n in range(N):
            inverse_matrix[i][n] /= scaler

    return matrix_inverse

Python-only implementation of the matrix inversion via Gaussian elimination

Summary

You are now familiar with how to calculate matrix inversion using the essential Python libraries for linear algebra. Also, if you cannot rely on any external libraries, can implement the Gaussian elimination processing for matrix inversion to find the matrix inverse without any external dependencies outside core Python.

Inverse Matrix Fundamentals

Test data for the matrix inversion

Calculating matrix inverse using NumPy

Calculating matrix inversion using SciPy

Custom algorithm for matrix inversion

Gaussian elimination

Example

Implementation in Python

Summary

Further reading