Background

I wrote this code back on the second semester of my freshman year in university for an advanced level assignment. At the time I was still very much inexperienced with programming and it took a lot of pains to convert ideas to code. We can easily find the inverse of a matrix using libraries such as numpy or scipy, but I thought doing so would miss the intent of the question and would not net me the extra grade. To this day I do not know if I did end up getting it or not.

I had to follow video lectures on Youtube to learn about LU (lower-upper) decomposition. You can find the videos at the bottom of this post.

LU decomposition is computationally cheaper over Gauss-Jordan elimination (also linked a video below). Having said that, the python code I created will always be slower compared to the equivalent scipy or numpy functions which are implemented in C++ using LAPACK. No beating that :)

The Code

Identity Matrix Function

This function is used to create an identity matrix. It will be used in the inverse matrix function to create an identity matrix that is the same size as the input.

1
2
3
4
5
6
7

def identity(m):
    result = []
    for i in range(m):
        row = [0]*m
        row[i] = 1
        result.append(row)
    return result

Inverse Matrix Function

Now the interesting part. As the name implies, we will first decompose our input nxn matrix into a lower matrix and an upper matrix. We do that using the code below. Note that a matrix in python is represented as a list of lists.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

def matrix_inverse(a):
    # accepts any square, non-singular matrix of dimension 2x2 or bigger
    #LU decomposition
    L = identity(len(a))
    for j in range(0, len(a[0])):
        # checks the lower triangular for non-zeros
        for i in range(j + 1, len(a)):
            if a[i][j] != 0:
                # divide non zero aij elements by ajj i.e. the diagonal element of that column
                mc = a[i][j] / a[j][j] 
                # create an eliminator list using mc to cancel out the non-zero
                eliminator = [a[j][x] * mc for x in range(len(a[j]))]
                trouble_row = [a[i][y] - eliminator[y] for y in range(len(a[i]))]
                a[i] = trouble_row
                L[i][j] = mc

    U = a

The first step of LU decomposition is Forward Elimination, which turns the lower triangular elements to zeros by creating, as I like to call it, eliminator rows using a multiplier calculated from the non-zero element divided by the diagonal matrix element directly above it. Subtract the row where the non-zero lower triangular element resides and cancel it out to 0. Repeat until we have our upper triangular matrix. This process of checking non-zeros goes column by column.

The multipliers used to create the subtraction rows are actually the elements of our lower triangular matrix. Their positions in the matrix are the same as the positions of the non-zero elements they helped cancel out. The diagonal elements of the lower triangular matrix are 1s and the upper triangle elements are 0s.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

    #forward-substitution
    idm = identity(len(a))
    # create a zero matrix with th same size as a
    Z = [[0 for y in range(len(U))] for x in range(len(L))]
    for f in range(len(L[0])):
        for r in range(len(L)):
            Z[r][f] += idm[r][f] - sum([L[r][x] * Z[x][f] for x in range(len(L[0]))])

    #back-substitution
    X = [[0 for y in range(len(U))] for x in range(len(L))]
    for e in range(len(X[0])):
        for w in range(-1, -(1 + len(U)), -1):
            for q in range(-1, -(1 + len(U[0])), -1):
                if U[w][q] == U[w][w]:
                    X[w][e] += Z[w][e] / U[w][q]
                else:
                    Z[w][e] -= U[w][q] * X[q][e]

    return X

work in progress…