Python Inverse of a Matrix

Question

How do I get the inverse of a matrix in python? I've implemented it myself, but it's pure python, and I suspect there are faster modules out there to do it.

Solution 1

You should have a look at numpy if you do matrix manipulation. This is a module mainly written in C, which will be much faster than programming in pure python. Here is an example of how to invert a matrix, and do other matrix manipulation.

from numpy import matrix
from numpy import linalg
A = matrix( [[1,2,3],[11,12,13],[21,22,23]]) # Creates a matrix.
x = matrix( [[1],[2],[3]] )                  # Creates a matrix (like a column vector).
y = matrix( [[1,2,3]] )                      # Creates a matrix (like a row vector).
print A.T                                    # Transpose of A.
print A*x                                    # Matrix multiplication of A and x.
print A.I                                    # Inverse of A.
print linalg.solve(A, x)     # Solve the linear equation system.

You can also have a look at the array module, which is a much more efficient implementation of lists when you have to deal with only one data type.

Solution 2

Make sure you really need to invert the matrix. This is often unnecessary and can be numerically unstable. When most people ask how to invert a matrix, they really want to know how to solve Ax = b where A is a matrix and x and b are vectors. It's more efficient and more accurate to use code that solves the equation Ax = b for x directly than to calculate A inverse then multiply the inverse by B. Even if you need to solve Ax = b for many b values, it's not a good idea to invert A. If you have to solve the system for multiple b values, save the Cholesky factorization of A, but don't invert it.

See Don't invert that matrix.

Solution 3

It is a pity that the chosen matrix, repeated here again, is either singular or badly conditioned:

A = matrix( [[1,2,3],[11,12,13],[21,22,23]])

By definition, the inverse of A when multiplied by the matrix A itself must give a unit matrix. The A chosen in the much praised explanation does not do that. In fact just looking at the inverse gives a clue that the inversion did not work correctly. Look at the magnitude of the individual terms - they are very, very big compared with the terms of the original A matrix...

It is remarkable that the humans when picking an example of a matrix so often manage to pick a singular matrix!

I did have a problem with the solution, so looked into it further. On the ubuntu-kubuntu platform, the debian package numpy does not have the matrix and the linalg sub-packages, so in addition to import of numpy, scipy needs to be imported also.

If the diagonal terms of A are multiplied by a large enough factor, say 2, the matrix will most likely cease to be singular or near singular. So

A = matrix( [[2,2,3],[11,24,13],[21,22,46]])

becomes neither singular nor nearly singular and the example gives meaningful results... When dealing with floating numbers one must be watchful for the effects of inavoidable round off errors.

Solution 4

For those like me, who were looking for a pure Python solution without pandas or numpy involved, check out the following GitHub project: https://github.com/ThomIves/MatrixInverse.

It generously provides a very good explanation of how the process looks like "behind the scenes". The author has nicely described the step-by-step approach and presented some practical examples, all easy to follow.

This is just a little code snippet from there to illustrate the approach very briefly (AM is the source matrix, IM is the identity matrix of the same size):

def invert_matrix(AM, IM):
    for fd in range(len(AM)):
        fdScaler = 1.0 / AM[fd][fd]
        for j in range(len(AM)):
            AM[fd][j] *= fdScaler
            IM[fd][j] *= fdScaler
        for i in list(range(len(AM)))[0:fd] + list(range(len(AM)))[fd+1:]:
            crScaler = AM[i][fd]
            for j in range(len(AM)):
                AM[i][j] = AM[i][j] - crScaler * AM[fd][j]
                IM[i][j] = IM[i][j] - crScaler * IM[fd][j]
    return IM

But please do follow the entire thing, you'll learn a lot more than just copy-pasting this code! There's a Jupyter notebook as well, btw.

Hope that helps someone, I personally found it extremely useful for my very particular task (Absorbing Markov Chain) where I wasn't able to use any non-standard packages.

Solution 5

Numpy will be suitable for most people, but you can also do matrices in Sympy

Try running these commands at http://live.sympy.org/

M = Matrix([[1, 3], [-2, 3]])
M
M**-1

For fun, try M**(1/2)

View more solutions

Comments

[futurecat]

numpy is also featured in the book "Beautiful Code". :-)

[Sakie]

Note here also, that there's no inversion happening, and that the system is solved directly, as per John D. Cook's answer.

[Kragen Javier Sitaker]

What if my matrix members are exact rationals? It seems like that avoid the accuracy problem, although of course at the cost of making the performance problem a lot worse.

[Juh_]

Never used R, but why would an external program and its python binder be better than the most well known scientific package of python?

[georg]

one may also check A==A.I.I in order to verifiy the result

[asmeurer]

The problem is that humans pick matrices at "random" by entering simple arithmetic progressions in the rows, like 1, 2, 3 or 11, 12, 13. The problem is that if you have at least three rows like this they are always linearly dependent.

[primo]

I found that numpy.linalg was giving inexact results for matrices containing large integers, whereas the results from sympy are exact. +1

[Deniz Kaplan]

With an approximate precision, Sympy is a good and live terminal. I checked with command (M**-1)*M and it gave unit matrix (not exactly but very close one)

[Praveen]

Consider using numpy arrays instead of matrices. Along with numpy.linalg, you can get pretty much everything you want, without using the unwieldy numpy.matrix.