Matrix Inverse Without Numpy

Question

I want to invert the matrix without using numpy.linalg.inv .

The reason is that I use Numba to speed up the code, but numpy.linalg.inv is not supported, so I wonder if I can invert the matrix with the "classic" Python code.

With numpy.linalg.inv, a sample code would look like this:

import numpy as np M = np.array([[1,0,0],[0,1,0],[0,0,1]]) Minv = np.linalg.inv(M) 

Answer 1

I used the formula http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html to write a function that performs an inversion of the 4x4 matrix:

 import numpy as np def myInverse(A): detA = np.linalg.det(A) b00 = A[1,1]*A[2,2]*A[3,3] + A[1,2]*A[2,3]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] - A[1,3]*A[2,2]*A[3,1] b01 = A[0,1]*A[2,3]*A[3,2] + A[0,2]*A[2,1]*A[3,3] + A[0,3]*A[2,2]*A[3,1] - A[0,1]*A[2,2]*A[3,3] - A[0,2]*A[2,3]*A[3,1] - A[0,3]*A[2,1]*A[3,2] b02 = A[0,1]*A[1,2]*A[3,3] + A[0,2]*A[1,3]*A[3,1] + A[0,3]*A[1,1]*A[3,2] - A[0,1]*A[1,3]*A[3,2] - A[0,2]*A[1,1]*A[3,3] - A[0,3]*A[1,2]*A[3,1] b03 = A[0,1]*A[1,3]*A[2,2] + A[0,2]*A[1,1]*A[2,3] + A[0,3]*A[1,2]*A[2,1] - A[0,1]*A[1,2]*A[2,3] - A[0,2]*A[1,3]*A[2,1] - A[0,3]*A[1,1]*A[2,2] b10 = A[1,0]*A[2,3]*A[3,2] + A[1,2]*A[2,0]*A[3,3] + A[1,3]*A[2,2]*A[3,0] - A[1,0]*A[2,2]*A[3,3] - A[1,2]*A[2,3]*A[3,0] - A[1,3]*A[2,0]*A[3,2] b11 = A[0,0]*A[2,2]*A[3,3] + A[0,2]*A[2,3]*A[3,0] + A[0,3]*A[2,0]*A[3,2] - A[0,0]*A[2,3]*A[3,2] - A[0,2]*A[2,0]*A[3,3] - A[0,3]*A[2,2]*A[3,0] b12 = A[0,0]*A[1,3]*A[3,2] + A[0,2]*A[1,0]*A[3,3] + A[0,3]*A[1,2]*A[3,0] - A[0,0]*A[1,2]*A[3,3] - A[0,2]*A[1,3]*A[3,0] - A[0,3]*A[1,0]*A[3,2] b13 = A[0,0]*A[1,2]*A[2,3] + A[0,2]*A[1,3]*A[2,0] + A[0,3]*A[1,0]*A[2,2] - A[0,0]*A[1,3]*A[2,2] - A[0,2]*A[1,0]*A[2,3] - A[0,3]*A[1,2]*A[2,0] b20 = A[1,0]*A[2,1]*A[3,3] + A[1,1]*A[2,3]*A[3,0] + A[1,3]*A[2,0]*A[3,1] - A[1,0]*A[2,3]*A[3,1] - A[1,1]*A[2,0]*A[3,3] - A[1,3]*A[2,1]*A[3,0] b21 = A[0,0]*A[2,3]*A[3,1] + A[0,1]*A[2,0]*A[3,3] + A[0,3]*A[2,1]*A[3,0] - A[0,0]*A[2,1]*A[3,3] - A[0,1]*A[2,3]*A[3,0] - A[0,3]*A[2,0]*A[3,1] b22 = A[0,0]*A[1,1]*A[3,3] + A[0,1]*A[1,3]*A[3,0] + A[0,3]*A[1,0]*A[3,1] - A[0,0]*A[1,3]*A[3,1] - A[0,1]*A[1,0]*A[3,3] - A[0,3]*A[1,1]*A[3,0] b23 = A[0,0]*A[1,3]*A[2,1] + A[0,1]*A[1,0]*A[2,3] + A[0,3]*A[1,1]*A[2,0] - A[0,0]*A[1,1]*A[2,3] - A[0,1]*A[1,3]*A[2,0] - A[0,3]*A[1,0]*A[2,1] b30 = A[1,0]*A[2,2]*A[3,1] + A[1,1]*A[2,0]*A[3,2] + A[1,2]*A[2,1]*A[3,0] - A[1,0]*A[2,1]*A[3,2] - A[1,1]*A[2,2]*A[3,0] - A[1,2]*A[2,0]*A[3,1] b31 = A[0,0]*A[2,1]*A[3,2] + A[0,1]*A[2,2]*A[3,0] + A[0,2]*A[2,0]*A[3,1] - A[0,0]*A[2,2]*A[3,1] - A[0,1]*A[2,0]*A[3,2] - A[0,2]*A[2,1]*A[3,0] b32 = A[0,0]*A[1,2]*A[3,1] + A[0,1]*A[1,0]*A[3,2] + A[0,2]*A[1,1]*A[3,0] - A[0,0]*A[1,1]*A[3,2] - A[0,1]*A[1,2]*A[3,0] - A[0,2]*A[1,0]*A[3,1] b33 = A[0,0]*A[1,1]*A[2,2] + A[0,1]*A[1,2]*A[2,0] + A[0,2]*A[1,0]*A[2,1] - A[0,0]*A[1,2]*A[2,1] - A[0,1]*A[1,0]*A[2,2] - A[0,2]*A[1,1]*A[2,0] Ainv = np.array([[b00, b01, b02, b03], [b10, b11, b12, b13], [b20, b21, b22, b23], [b30, b31, b32, b33]]) / detA return Ainv 

Answer 2

Here is a more elegant and scalable solution, imo. It will work for any nxn matrix, and you can use it for other methods. Note that getMatrixInverse (m) takes an array of arrays as an array. Please feel free to ask any questions.

 def transposeMatrix(m): return map(list,zip(*m)) def getMatrixMinor(m,i,j): return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])] def getMatrixDeternminant(m): #base case for 2x2 matrix if len(m) == 2: return m[0][0]*m[1][1]-m[0][1]*m[1][0] determinant = 0 for c in range(len(m)): determinant += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c)) return determinant def getMatrixInverse(m): determinant = getMatrixDeternminant(m) #special case for 2x2 matrix: if len(m) == 2: return [[m[1][1]/determinant, -1*m[0][1]/determinant], [-1*m[1][0]/determinant, m[0][0]/determinant]] #find matrix of cofactors cofactors = [] for r in range(len(m)): cofactorRow = [] for c in range(len(m)): minor = getMatrixMinor(m,r,c) cofactorRow.append(((-1)**(r+c)) * getMatrixDeternminant(minor)) cofactors.append(cofactorRow) cofactors = transposeMatrix(cofactors) for r in range(len(cofactors)): for c in range(len(cofactors)): cofactors[r][c] = cofactors[r][c]/determinant return cofactors 

Answer 3

For a 4x4 matrix, it’s probably just OK to use the mathematical formula that you can find using the Google formula for the 4 by 4 matrix. For example, here (I cannot vouch for its accuracy):

http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html

In general, inverting a common matrix is ​​not for the faint of heart. You must know all mathematically complex cases and know why they will not be applied to your use, and catch them when you receive mathematically pathological inputs (which either return results with low accuracy or numerical garbage in the knowledge that it does not matter in your use case, if you don’t actually divide by zero or overflow MAXFLOAT ... which you can catch with an exception handler and present it as "Error: the matrix is ​​the only or very close to it").

As a rule, it is better for a programmer to use library code written by specialists in mathematical mathematics if you do not want to spend time understanding the physical and mathematical nature of the specific problem you are addressing and becoming your own mathematics expert in your own specialist field.

Answer 4

At least July 16, 2018, Numba has a fast inverse matrix. (You can see how they overload the standard NumPy inverse and other operations here .)

Here are the results of my benchmarking:

 import numpy as np from scipy import linalg as sla from scipy import linalg as nla import numba def gen_ex(d0): x = np.random.randn(d0,d0) return xT + x @numba.jit def inv_nla_jit(A): return np.linalg.inv(A) @numba.jit def inv_sla_jit(A): return sla.inv(A) 

For small matrices, this is especially fast:

 ex1 = gen_ex(4) %timeit inv_nla_jit(ex1) # NumPy + Numba %timeit inv_sla_jit(ex1) # SciPy + Numba %timeit nla.inv(ex1) # NumPy %timeit sla.inv(ex1) # SciPy 

[Of]

 2.54 µs ± 467 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each) 67.3 µs ± 9.18 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 63.5 µs ± 7.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 56.6 µs ± 5.03 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

Note that acceleration only works for inverting NumPy, not SciPy (as expected).

Slightly larger matrix:

 ex2 = gen_ex(40) %timeit inv_nla_jit(ex2) # NumPy + Numba %timeit inv_sla_jit(ex2) # SciPy + Numba %timeit nla.inv(ex2) # NumPy %timeit sla.inv(ex2) # SciPy 

[Of]

 131 µs ± 12.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 278 µs ± 26.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) 231 µs ± 24.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) 189 µs ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) 

So there is still acceleration, but SciPy is catching up.