Inverting matrix in python slightly off

Question

I'm trying to use this http://www.irma-international.org/viewtitle/41011/ algorithm to invert a nxn matrix.

I ran the function on this matrix

[[1.0, -0.5],

 [-0.4444444444444444, 1.0]]

and got the output

[[ 1.36734694,  0.64285714]

 [ 0.57142857,  1.28571429]]

the correct output is meant to be

[[ 1.28571429,  0.64285714]

 [ 0.57142857,  1.28571429]]

My function:

def inverse(m):
    n = len(m)
    P = -1
    D =  1
    mI = m
    while True:
        P += 1
        if m[P][P] == 0:
           raise Exception("Not Invertible")
        else:
            D = D * m[P][P]
            for j in range(n):
                if j != P:
                    mI[P][j] =  m[P][j] / m[P][P]
            for i in range(n):
                if i != P:
                    mI[i][P] = -m[i][P] / m[P][P]
            for i in range(n):
                for j in range(n):
                    if i != P and j != P:
                        mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
            mI[P][P] = 1 / m[P][P]
            if P == n - 1: # All elements have been looped through
                break

    return mI

Where am I making my mistake?

You may be interested in this: https://stackoverflow.com/a/39881366/4909087 — cs95, Feb 16 '18 at 04:13
@cᴏʟᴅsᴘᴇᴇᴅ Thanks for the link. I actually saw that one but am trying to get this algorithm to work. Theres something short and simple about this one that I like. — huddie96, Feb 16 '18 at 04:18
Do you have a link to the entire paper, It almost seems like the leema at the end of the page is setting up for a special case in 2x2 which is fairly common. — gdahlm, Feb 16 '18 at 04:49
@gdahlm I don't sorry. My mistake was pointed out beautifully below — huddie96, Feb 16 '18 at 04:59
It looks to me that this is (related to) the SWEEP method, used in statistics. Note though that it will not invert all invertible matrices, because of the requirement that M(0,0) be non-zero; it will fail to invert rotation by 90 degrees for example. — dmuir, Feb 16 '18 at 08:55

TTT · Accepted Answer · 2018-02-16T04:53:05.440

https://repl.it/repls/PowerfulOriginalOpensoundsystem

Output

inverse: [[Decimal('1.285714285714285693893862813'), Decimal('0.6428571428571428469469314065')], [Decimal('0.5714285714285713877877256260'), Decimal('1.285714285714285693893862813')]] numpy: [[ 1.28571429 0.64285714] [ 0.57142857 1.28571429]]

from decimal import Decimal
import numpy as np

def inverse(m):
    m = [[Decimal(n) for n in a] for a in m]
    n = len(m)
    P = -1
    D =  Decimal(1)
    mI = [[Decimal(0) for n in a] for a in m]
    while True:
        P += 1
        if m[P][P] == 0:
           raise Exception("Not Invertible")
        else:
            D = D * m[P][P]
            for j in range(n):
                if j != P:
                    mI[P][j] =  m[P][j] / m[P][P]
            for i in range(n):
                if i != P:
                    mI[i][P] = -m[i][P] / m[P][P]
            for i in range(n):
                for j in range(n):
                    if i != P and j != P:
                      mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
            mI[P][P] = 1 / m[P][P]
            m = [[Decimal(n) for n in a] for a in mI]
            mI = [[Decimal(0) for n in a] for a in m]
            if P == n - 1: # All elements have been looped through
                break

    return m

m = [[1.0, -0.5],

 [-0.4444444444444444, 1.0]]

print(inverse(m))
print(np.linalg.inv(np.array(m)))

My thought process:

At first, I thought you might have lurking floating point roundoff errors. This turned out to not be true. That's what the Decimal jazz is for.

Your bug is here

mI = m # this just creates a pointer that points to the SAME list as m

and here

for i in range(n):
                for j in range(n):
                    if i != P and j != P:
                        mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
            mI[P][P] = 1 / m[P][P] 
            # you are not copying mI to m for the next iteration
            # you are also not zeroing mI 
            if P == n - 1: # All elements have been looped through
                break

    return mI

In adherence to the algorithm, every iteration creates a NEW a' matrix, it does not continue to modify the same old a. I inferred this to mean that in the loop invariant, a becomes a'. Works for your test case, turns out to be true.

Be careful of making reference copies of containers in python! — TTT, Feb 16 '18 at 04:59

Inverting matrix in python slightly off

1 Answers1