Photo by Markus Spiske / Unsplash

Intuitive method of finding the inverse of a matrix

Python Feb 11, 2022

Finding a matrix inverse is very useful not only in Computer science but when solving equations by hand. However, using the formula Adj(A) / det(A) can be quite tedious, especially if you don't know how exactly it works. In this article I will show how to intuitively and efficiently find the inverse of a matrix.

Let us begin with the following equation:

A * I = A

By a series of changes to matrices `I` and `A`, we want to reach the following equation:

A * A^-1 = I

The changes we can perform are:

  • multiplying a row by a number
  • adding one row multiplied by a number to another

Each change must be made to both matrices for the equation to still be true.

It is easier to think about getting from `A` to `I`, than it is about getting from `I` to `A^-1` because we already know what `A` and `I` look like, so that's where we'll start.

We'll take the following matrix A for example:

    [ 2 3 2 ]
A = [ 3 2 3 ]
    [ 2 3 3 ]
    
    [ 1 0 0 ]
I = [ 0 1 0 ]
    [ 0 0 1 ]

Let us multiply the first row with the inverse of the first element of the first row

     [ 1  1.5  1 ]
A' = [ 3   2   3 ]
     [ 2   3   3 ]
     
     [ 0.5  0  0 ]
I' = [  0   1  0 ]
     [  0   0  1 ]

Notice that now we can substract the first row times 3 from the second row to get 0 in the position (2, 1)

     [ 1  1.5  1 ]
A' = [ 0 -2.5  0 ]
     [ 2   3   3 ]

     [ 0.5  0  0 ]
I' = [-1.5  1  0 ]
     [  0   0  1 ]

We can do the same for the third row

     [ 1  1.5  1 ]
A' = [ 0 -2.5  0 ]
     [ 0   0   1 ]
     
     [ 0.5  0  0 ]
I' = [-1.5  1  0 ]
     [ -1   0  1 ]

Now the first column of A' is the first column from the identity matrix. In the next iteration we repeat the same process, but instead we set the position (2, 2) to 1 and all others to 0.

     [ 1  1.5  1 ]
A' = [ 0   1   0 ]
     [ 0   0   1 ]
     
     [ 0.5   0   0 ]
I' = [ 0.6 -0.4  0 ]
     [ -1    0   1 ]
     [ 1  0  1 ]
A' = [ 0  1  0 ]
     [ 0  0  1 ]
     
     [-0.4  0.6  0 ]
I' = [ 0.6 -0.4  0 ]
     [ -1    0   1 ]

Lastly, you only need to substract the last row from the first row to transform A' to I.

     [ 1  0  0 ]
A' = [ 0  1  0 ] = I
     [ 0  0  1 ]
     
     [ 0.6  0.6 -1 ]
I' = [ 0.6 -0.4  0 ] = A^-1
     [ -1    0   1 ]

As you can easily check by matrix multiplication, we have found the inverse.

The method works because once you solve a column, it can no longer be affected by future changes (because all fields are zero except for one, which will not be added or substracted from anything anymore) and since all changes are made to both matrices, once A' becomes I, I' will become A^-1.

The algorithm is:

Iterate over columns
    let i denote index of the column
    divide row i by value at position (i, i)
    iterate over all rows except for row i
    	let j denote index of the row
        substract row i times value at position (j, i) from row j

Below is an implementation of the algorithm in Python.

# Implementation of the algorithm
# in Python 3

# multiply matrices
def multiply(a, b):
    res = []
    for i in range(len(a)):
        rez.append([])
        for j in range(len(a)):
            res[i].append(0)
            for k in range(len(a)):
                res[i][j] += a[i][k] * b[k][j]
    return res


# divide a row with a number
def divide(row, num):
    k = 1 / num
    for i in range(len(row)):
        row[i] *= k

# add one row times a factor to another
def addrow(row, row2, f):
    for i in range(len(row)):
        row[i] += row2[i] * f

# print out a matrix
def printmatrix(m):
    for i in range(len(m)):
        for j in range(len(m)):
            print('\t{}'.format(round(m[i][j], 1)), end=' ')
        print()
    print()

# calculate inverse of a matrix
def inv(m2):
    # create a copy of the matrix
    m = [[m2[j][i] for i in range(len(m2))] for j in range(len(m2))]
    # identity matrix
    im = [[0 for j in range(len(m))] for i in range(len(m))]
    for i in range(len(m)):
        im[i][i] = 1
    # solve column by column
    for i in range(len(m)):
        # set m[i][i] to 1
        divide(im[i], m[i][i])
        divide(m[i], m[i][i])
        for j in range(len(m)):
            if j == i: continue
            # set m[j][i] to 0
            addrow(im[j], im[i], -m[j][i])
            addrow(m[j], m[i], -m[j][i])
    return im

## TEST ##

ma = [
    [2, 3, 2],
    [3, 2, 3],
    [2, 3, 3]
]

iv = inv(ma)

print("MATRIX: ")

printmatrix(ma)

print('INVERSE: ')

printmatrix(iv)

Tags