-1

i'm using pyomo and need a function that gives me the inverse matrix of a squared matrix without using numpy. I already found a function online but can't find the issue that results in the following error message:

operands could not be broadcast together with shapes (0,4) (3,4)

def transposeMatrix(m):
    return list(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


a = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [10, 4, 6, 2]])
getMatrixInverse(a)
Paul Rogge
  • 17
  • 5
  • 2
    You are asking "without `np`", but are using `np`... Which is it? – Reblochon Masque Jul 08 '21 at 16:08
  • 1
    Now's a good time to learn [how to debug small programs](//ericlippert.com/2014/03/05/how-to-debug-small-programs/) and [to use a debugger](//stackoverflow.com/q/25385173/843953) Step through your code and observe what each line of code does. Identify where your program differs from your expectations by comparing these intermediate results with expected results. Work backwards from there to narrow down the cause of the problem. Then ask a _specific_ question if you are still confused by your code's behavior. It's not okay to dump your code and expect other people to debug it for you. – Pranav Hosangadi Jul 08 '21 at 16:54
  • 1
    How big might your matrix be? The cofactor method is O(N!) (yes, N factorial) for NxN matrices and so is impractical for any but the smallest matrices; 4x4 should be ok. Do you know anthing about the matrices? For example is it a covariance (and so positive definite)? – dmuir Jul 08 '21 at 17:00

1 Answers1

1

If we get rid of the NumPy reference and just use

a = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [10, 4, 6, 2]]

then the code works. However, your example has a zero determinant so you will get a divide by zero error when calculating the inverse. If instead you use a matrix like

a = [[ 0., 2.0, 3.0, 4.0 ], [ 5.0, 6.0, 7.0, 8. ],[ 9.0, 10.0, 11.0, 12.0 ], [ 13.,14.,15., 0. ]]

which has a determinant of -64, then you get the correct inverse

[[-1.0, 2.0, -1.0, 0.0], [2.0, -6.8125, 3.875, -0.0625], [-1.0, 4.625, -2.75, 0.125], [0.0, -0.0625, 0.125, -0.0625]]
Salix alba
  • 7,536
  • 2
  • 32
  • 38