Matrix inversion without Numpy

Question

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

The reason is that I am using Numba to speed up the code, but numpy.linalg.inv is not supported, so I am wondering if I can invert a matrix with 'classic' Python code.

With numpy.linalg.inv an example code would look like that:

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

Probably not. There's no python "builtin" doing that for you and programming a matrix inversion yourself is anything but easy (see e.g. https://en.wikipedia.org/wiki/Invertible_matrix#Methods_of_matrix_inversion for a probably incomplete list of methods). I'm also not aware of any `numpy`-independent linear algebra package for python... — sebastian, Aug 20 '15 at 09:14
If you want to invert 3x3 matrices only, you can look up the formula [here](http://mathworld.wolfram.com/MatrixInverse.html). (You better specify the dimension and type of matrices you want to invert. In your example you use the most trivial identity matrix. Are they real? And regular?) — Falko, Aug 20 '15 at 09:17

stackPusher · Answer 1 · 2017-11-10T07:13:17.713

84

Here is a more elegant and scalable solution, imo. It'll work for any nxn matrix and you may find use for the other methods. Note that getMatrixInverse(m) takes in an array of arrays as input. 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

edited Nov 10 '17 at 07:13

answered Oct 05 '16 at 18:34

stackPusher

6,076
3
35
36

1

This works perfectly. According to the requirement, should be the accepted answer. The only minor change required is in `#base case for 2x2 matrix`. you need to explicitly convert to float. – Aju Antony Jun 09 '17 at 07:14
1

If the matrix is not square the transpose function will give an error, to find the transpose for a list simply we can do: zip(*theArray) Taken from: https://stackoverflow.com/questions/4937491/matrix-transpose-in-python – Mohanad Kaleia Nov 08 '17 at 18:59
2

@MohanadKaleia you're right, thanks. Although non square matrices don't have inverses, I do claim my answer is composed of reusable pieces so i've fixed the transpose function as per your suggestion. – stackPusher Nov 10 '17 at 07:19
1

@stackPusher I used part of your code with a non square matrix to find the projection matrix of it, that is why I had this issue. Thank you – Mohanad Kaleia Nov 11 '17 at 01:09
3

@stackPusher this is tremendous. I wish I could upvote more than once – fermi Jul 26 '18 at 22:23
2

@stackPusher I am getting this error on your code. Can you please see.. in getMatrixMinor(m, i, j) 3 4 def getMatrixMinor(m,i,j): ----> 5 return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])] 6 7 def getMatrixDeternminant(m): ValueError: operands could not be broadcast together with shapes (0,172877) (172876,172877) – shaifali Gupta Mar 23 '19 at 09:57
far from a modern method of computing the matrix inverse. – Faris Sbahi Apr 28 '19 at 06:35
what would be its time complexity ? I am getting it exceeds 5 sec for a 10*10 matrix – Subhang V May 07 '20 at 16:28
9

If you're using python3, then you need to define `transposeMatrix` as `list(map(list,zip(*m)))` instead of `map(list,zip(*m))` – TheChetan Aug 18 '20 at 05:44
1

@SubhangV, If you cache the `Deternminant` values for each matrix, then the algorithm speeds up tremendously – TheChetan Aug 19 '20 at 05:22
1

this solution only works for the 2x2 case when I try. otherwise i get errors like ValueError: operands could not be broadcast together with shapes (0,4) (3,4) – Ayan Mitra Apr 28 '22 at 04:08
This is a super slow way to calculate a matrix inverse. Even with memoization, it's O(n^4) time. Without memoization (like in the above code), it's O(n!^2) time. Good luck inverting anything larger than a 10x10 matrix! – Joseph Camacho Jun 08 '22 at 18:33
@shaifaliGupta you are likely passing in a numpy object which behaves differently to lists try passing in your matrix with .tolist() https://stackoverflow.com/q/47257427/8185373 https://numpy.org/doc/stable/reference/generated/numpy.ndarray.tolist.html – Joss Bird Oct 31 '22 at 16:19

Asad-ullah Khan · Answer 2 · 2020-12-30T18:02:32.793

11

Here is another way, using gaussian elimination instead:

def eliminate(r1, r2, col, target=0):
    fac = (r2[col]-target) / r1[col]
    for i in range(len(r2)):
        r2[i] -= fac * r1[i]

def gauss(a):
    for i in range(len(a)):
        if a[i][i] == 0:
            for j in range(i+1, len(a)):
                if a[i][j] != 0:
                    a[i], a[j] = a[j], a[i]
                    break
            else:
                raise ValueError("Matrix is not invertible")
        for j in range(i+1, len(a)):
            eliminate(a[i], a[j], i)
    for i in range(len(a)-1, -1, -1):
        for j in range(i-1, -1, -1):
            eliminate(a[i], a[j], i)
    for i in range(len(a)):
        eliminate(a[i], a[i], i, target=1)
    return a

def inverse(a):
    tmp = [[] for _ in a]
    for i,row in enumerate(a):
        assert len(row) == len(a)
        tmp[i].extend(row + [0]*i + [1] + [0]*(len(a)-i-1))
    gauss(tmp)
    ret = []
    for i in range(len(tmp)):
        ret.append(tmp[i][len(tmp[i])//2:])
    return ret

edited Dec 30 '20 at 18:02

answered May 11 '20 at 23:18

Asad-ullah Khan

1,573
18
22

7

I required this technique to solve a Markov chain. – quartarian Dec 27 '20 at 12:07
1

Ha! That was the reason I made this as well – Asad-ullah Khan Dec 28 '20 at 15:32
7

foobar challenge? – luizfls Mar 12 '21 at 14:52
yup you got it! – Asad-ullah Khan Mar 13 '21 at 17:38
This seems more efficient than stackPusher's answer, right? – Mercury Feb 14 '22 at 21:34
I think Big-O is the same, maybe my constant is lower? – Asad-ullah Khan Feb 24 '22 at 20:04
More foobar here. :) – rrrrrrrrrrrrrrrr Sep 15 '22 at 15:07
The function `gauss()` return an unused value (`tmp`) – Yogev Neumann Nov 26 '22 at 23:34
@YogevNeumann gauss modifies `tmp` in place and also returns it – Asad-ullah Khan Nov 28 '22 at 00:57
@rrrrrrrrrrrrrrrr Did you get a callback from Google after finishing the challenge? Or is this just an easter egg? – desert_ranger Jul 03 '23 at 23:33
@desert_ranger Haven't gotten around to finishing the challenge just yet. :) – rrrrrrrrrrrrrrrr Jul 19 '23 at 08:00
Great answer, I rank this higher than the accepted answer. – Thomas Kimber Aug 30 '23 at 20:28

score 9 · Answer 3 · edited Aug 23 '21 at 17:26

As of at least July 16, 2018 Numba has a fast matrix inverse. (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 x.T + 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 it is particularly 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

[Out]

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)

Notice that the speedup only works for NumPy inverse, 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

[Out]

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's still a speedup here but SciPy is catching up.

score 2 · Answer 4 · answered Aug 20 '15 at 10:23

For a 4 x 4 matrix it's probably just about OK to use the mathematical formula, which you can find using Googling "formula for 4 by 4 matrix inverse". For example here (I can't vouch for its accuracy):

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

In general inverting a general matrix is not for the faint-hearted. You have to be aware of all the mathematically difficult cases and know why they won't apply to your usage, and catch them when you are supplied with mathematically pathological inputs (that, or return results of low accuracy or numerical garbage in the knowledge that it won't matter in your usage case provided you don't actually end up dividing by zero or overflowing MAXFLOAT ... which you might catch with an exception handler and present as "Error: matrix is singular or very close thereto").

It's generally better as a programmer to use library code written by numerical mathematics experts, unless you are willing to spend time understanding the physical and mathematical nature of the particular problem that you are addressing and become your own mathematics expert in your own specialist field.

Moved to http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche0023.html — Csaba Toth, Feb 10 '22 at 00:37

score 2 · Answer 5 · answered Jan 07 '22 at 02:51

I found that Gaussian Jordan Elimination Algorithm helped a lot when attempting this. If you're going to use a given matrix (any size, i.e 5x5) where the hardcore formula for it is 49 pages long. It's best to use this. To inverse a matrix place it as a 2D array and then run the Inverse function

# Python test Guassion Jordan Elimination
# Inputs are 2D array not matrix

Test_Array = [[3,3,2,1,1],[2,1,3,2,3],[1,3,3,2,2],[2,3,3,1,1], 
[3,1,2,1,2]]

# Creating storage & initalizing for augmented matrix
# this is the same as the np.zeros((n,2*n)) function
def nx2n(n_Rows, n_Columns):
    Zeros = []
    for i in range(n_Rows):
        Zeros.append([])
        for j in range(n_Columns*2):
            Zeros[i].append(0)
    return Zeros

# Applying matrix coefficients
def update(inputs, n_Rows, n_Columns, Zero):
    for i in range(n_Rows):
        for j in range(n_Columns):
            Zero[i][j] = inputs[i][j]
    return Zero

# Augmenting Identity Matrix of Order n
def identity(n_Rows, n_Columns, Matrix):
    for i in range(n_Rows):
        for j in range(n_Columns):
            if i == j:
                Matrix[i][j+n_Columns] = 1
    return Matrix

# Applying & implementing the GJE algorithm
def Gussain_Jordan_Elimination(n_Rows, n_Columns, Matrix):
    for i in range(n_Rows):
        if Matrix[i][i] == 0:
            print('error cannot divide by "0"')
    
        for j in range(n_Columns):
            if i != j:
                ratio = Matrix[j][i]/Matrix[i][i]

                for k in range(2*n_Columns):
                    Matrix[j][k] = Matrix[j][k] - ratio * Matrix[i][k]
    return Matrix

# Row Operation to make Principal Diagonal Element to '1'
def row_op(n_Rows, n_Columns, Matrix):
    for i in range(n_Rows):
        divide = Matrix[i][i]
        for j in range(2*n_Columns):
            Matrix[i][j] = Matrix[i][j]/divide
    return Matrix

# Display Inversed Matix
def Inverse(Matrix):
    returnable = []
    number_Rows = int(len(Matrix))
    number_Columns = int(len(Matrix[0]))
    Inversed_Matrix = (row_op(number_Rows, number_Columns, 
        Gussain_Jordan_Elimination(number_Rows, number_Columns, 
        identity(number_Rows, number_Columns, 
        update(Matrix, number_Rows, number_Columns, 
        nx2n(number_Rows, number_Columns))))))

    for i in range(number_Rows):
        returnable.append([])
        for j in range(number_Columns, 2*number_Columns):
            returnable[i].append(Inversed_Matrix[i][j])
    return returnable

print(Inverse(Test_Array))

score 0 · Answer 6 · answered Nov 04 '21 at 19:23

Simply add all methods

import math

def getMinorIndex(matrixLocal, x, y):
  minor = []
  for i in range(3):
    minorRow = []
    if i == x:
      continue
    for j in range(3):
      if j == y:
        continue
      minorRow.append(matrixLocal[i][j])
    minor.append(minorRow)
  return minor

def getDeterminant2By2(matrixLocal):
  determinant = matrixLocal[0][0] * matrixLocal[1][1] - matrixLocal[0][1] * matrixLocal[1][0]
  return determinant

def getDeterminant(matrixLocal):
  determinant = 0
  for x in range(3):
    t = getDeterminant2By2(getMinorIndex(matrixLocal, 0, x))
    e = matrixLocal[0][x]
    determinant += (t * e * math.pow(-1, x))
  return determinant

def getCofactorMatrix(matrixLocal):
  cofactorMatrix = []
  for i in range(3):
    row = []
    for j in range(3):
      e = matrixLocal[i][j]
      t = getDeterminant2By2(getMinorIndex(matrixLocal, i, j))
      row.append(t * math.pow(-1, i + j))
    cofactorMatrix.append(row)
  return cofactorMatrix

def transpose(matrixLocal):
  transposeMatrix = []
  for i in range(3):
    row = []
    for j in range(3):
      e = matrixLocal[j][i]
      row.append(e)
    transposeMatrix.append(row)
  return transposeMatrix

def divideMatrix(matrixLocal, divisor):
  ansMatrix = []
  for i in range(3):
    row = []
    for j in range(3):
      e = matrixLocal[i][j]/divisor
      row.append(e)
    ansMatrix.append(row)
  return ansMatrix

cofactor = getCofactorMatrix(matrix)
adjoint = transpose(cofactor)
det = getDeterminant(matrix)
inverse = divideMatrix(adjoint, det)
inverse

score -1 · Answer 7 · answered May 28 '21 at 10:21

Inverse matrix of 3x3 without numpy [python3]

import pprint


def inverse_3X3_matrix():
    I_Q_list = [[0, 1, 1],
                [2, 3, -1],
                [-1, 2, 1]]
    det_ = I_Q_list[0][0] * (
            (I_Q_list[1][1] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][1])) - \
           I_Q_list[0][1] * (
                   (I_Q_list[1][0] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][0])) + \
           I_Q_list[0][2] * (
                   (I_Q_list[1][0] * I_Q_list[2][1]) - (I_Q_list[1][1] * I_Q_list[2][0]))
    co_fctr_1 = [(I_Q_list[1][1] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][1]),
                 -((I_Q_list[1][0] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][0])),
                 (I_Q_list[1][0] * I_Q_list[2][1]) - (I_Q_list[1][1] * I_Q_list[2][0])]

    co_fctr_2 = [-((I_Q_list[0][1] * I_Q_list[2][2]) - (I_Q_list[0][2] * I_Q_list[2][1])),
                 (I_Q_list[0][0] * I_Q_list[2][2]) - (I_Q_list[0][2] * I_Q_list[2][0]),
                 -((I_Q_list[0][0] * I_Q_list[2][1]) - (I_Q_list[0][1] * I_Q_list[2][0]))]

    co_fctr_3 = [(I_Q_list[0][1] * I_Q_list[1][2]) - (I_Q_list[0][2] * I_Q_list[1][1]),
                 -((I_Q_list[0][0] * I_Q_list[1][2]) - (I_Q_list[0][2] * I_Q_list[1][0])),
                 (I_Q_list[0][0] * I_Q_list[1][1]) - (I_Q_list[0][1] * I_Q_list[1][0])]

    inv_list = [[1 / det_ * (co_fctr_1[0]), 1 / det_ * (co_fctr_2[0]), 1 / det_ * (co_fctr_3[0])],
                [1 / det_ * (co_fctr_1[1]), 1 / det_ * (co_fctr_2[1]), 1 / det_ * (co_fctr_3[1])],
                [1 / det_ * (co_fctr_1[2]), 1 / det_ * (co_fctr_2[2]), 1 / det_ * (co_fctr_3[2])]]

    pprint.pprint(inv_list)


inverse_3X3_matrix()

score -8 · Accepted Answer · answered Aug 20 '15 at 10:37

I used the formula from http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html to write the function that does the inversion of a 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

You don't want to use `np.linalg.inv` but `np.linalg.det` is fine? That's a really awkward requirement... — sebastian, Aug 20 '15 at 11:15
Of course one needs to write another 'brute force' implementation for the determinant calculation as well. Or just calculate the det outside the Numba function and pass it as an argument — Alessandro Vianello, Aug 20 '15 at 12:27
Moved to http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche0023.html — Csaba Toth, Feb 10 '22 at 00:36

Matrix inversion without Numpy

8 Answers8

Linked

Related