# How to Calculate Matrix Inverse in Python

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.

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

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

## 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.

### Further reading

- Wikipedia article on matrix inverse
- Wikipedia article on Gaussian elimination
- How to calculate singular value decomposition using Python
- How to calculate eigenvalue decomposition using Python
- How to calculate Cholesky decomposition using Python
- Reference manual for
`numpy.linalg.inv`

- Reference manual for
`scipy.linalg.inv`