Strange result in math.inv() function

[isma4444]

Somehow I get different results using the inverse of matrices with python and with math.js

Original matrix:

A =
6.123233995736766e-17 1 0 788 
-1 6.123233995736766e-17 0 692   
0 0 1 0   
0 0 0 1

The result by using math.js

math.inv(A)=
0 -1 0 0 
1 6.123233995736766e-17 0 -788  
0 0 1 0
0 0 0 1

Using python:

[[ 6.123234e-17 -1.000000e+00 -0.000000e+00  6.920000e+02]
 [ 1.000000e+00  6.123234e-17  0.000000e+00 -7.880000e+02]
 [ 0.000000e+00  0.000000e+00  1.000000e+00  0.000000e+00]
 [ 0.000000e+00  0.000000e+00  0.000000e+00  1.000000e+00]]

The element (0,3) is different and by a big margin.

[harrysarson]

For reference in MATLAB I get the same as python:

      6.12323399573677e-17                        -1                         0                       692
                         1      6.12323399573677e-17                         0                      -788
                         0                         0                         1                         0
                         0                         0                         0                         1

The matrix A is very nearly singular with is causing these problems. I dunno if these errors can be avoided or not.

@isma4444 what does Ainv = math.inv(A); math.inv(Ai) give?

[ericman314]

Changing element (0,0) to exactly 0 results in the inv(A) giving the correct inverse. I suspect this is not an issue with a singular matrix but an actual bug related to some conditional check equating that element to 0 during gaussian elimination or something like that. I'm intrigued and interested in finding the bug but I won't be able to get to it for a day or two, likely.

[ericman314]

I managed to implement a Gauss-Jordan elimination with full pivoting in inv.js which solves this issue. Unfortunately, I shamelessly copied the algorithm from Numerical Recipes in C so we probably can't include it in math.js. I will need to study the algorithm and try to reproduce it from memory.

[josdejong]

Thanks for looking into this Eric. Maybe there is an open source implementation of the Gauss-Jordan elimination that you could port, like from Octave or so?

[harrysarson]

numpy must have something aswell, might be worth a look

[ericman314]

Thanks, those are indeed good resources. But it ended up being simple enough to add pivoting to the current implementation; see PR #1114.

[josdejong]

This issue is fixed now in v4.4.1. Thanks again Eric!