Wednesday, June 3, 2015

On Matrix Inversion in Python

In the introductory programming course that I lecture, the big programming project that I set for the students involved writing a number of matrix operations in pure Python (that is, using built-in functions only - I specifically disallowed the use of any imported modules like NumPy). Transposition and matrix multiplication are fairly simple to do, but the challenge is in writing a function to invert a matrix.

Of course, Gauss-Jordan elimination (which should be covered in any first year algebra course) is one of the easier ways to do this. I wrote a fairly long but simple-to-understand function to give out to the students as a model solution at the end of the project, but I was wondering, how short could a working matrix inversion function be?

I came up with the function below, complete with partial pivoting, which seems to work correctly. It could probably be even shorter, but I think I achieved what I set out to do:

The input argument A is simply a nested list, such as
    A = [[1,2,3],[1,5,9],[7,5,2]]

Line 3 simply augments the matrix with the identity matrix. The for loop on line 4 loops through each row, and sets the element on the main diagonal to 1, and every element above and below that element to zero. This happens in three steps. The first step, on line 5, is a simple pivoting step. It sorts all of the remaining rows such that the current row will be divided by the number with the largest possible magnitude (to avoid division by zero, and improve the accuracy if the matrix contains very small numbers). The second step (line 6) is to divide each element in the row by the element on the main diagonal, leaving a 1 on the main diagonal. The third step (line 7) is a little more complicated. For each row that is not the current row, we subtract the current row multiplied by the element in the column index equal to the index of the current row. This will set each element in the current column to zero, except for the element in the current row, which remains 1.

The net result is that the left hand side of the augmented matrix becomes the identity matrix, while the right hand side becomes the inverse. All that's left is to separate the two, and return the inverse. This happens on line 8.

I don't know who would find this post useful, because it's better to just use numpy, but it is certainly a fun exercise.

If you enjoyed this post, then don't forget to like, tweet, +1, or upvote on reddit. If you have any questions, comments or complaints, post them using the form below.
. . . . . . . . . . . . . . . . . . . . . . . .


Unknown said...

Thanks for your hard work, I used your solution in a stats hacking challenge to maintain a personal pure-python goal. Once you import numpy for one thing, it's hard to keep from using all of it.

Ryan Connor said...

Awesome. Was using another pure python implementation I found, but it started failing. This did the trick!