Unimodular Matrix Exact Inverse Python

Question

numpy has a method to invert matrices. However i have a matrix with huge integer entries and the matrix is unimodular, so the inverse is also an integer matrix. The error is too big just to use numpy to calculate it with float and then round the values. Is there a way to calculate the entries exact with integers or any other library?

Answer 1

You can use sympy.

For example, here's a numpy integer array:

In [148]: m = np.array([[10, 11, 0, 0], [9, 10, 0, 0], [100, 0, 20, 3], [-10000, 200, 133, 20]])

In [149]: m
Out[149]: 
array([[    10,     11,      0,      0],
       [     9,     10,      0,      0],
       [   100,      0,     20,      3],
       [-10000,    200,    133,     20]])

The exact determinant is 1, but both np.linalg.det() and np.linalg.inv() introduce floating point errors:

In [150]: np.linalg.det(m)
Out[150]: 0.99999999999887024

In [151]: np.linalg.inv(m)
Out[151]: 
array([[  1.00000000e+01,  -1.10000000e+01,  -6.69432763e-16,
          8.92627614e-17],
       [ -9.00000000e+00,   1.00000000e+01,   6.34413157e-16,
         -8.43354404e-17],
       [ -3.25400000e+05,   3.58000000e+05,   2.00000000e+01,
         -3.00000000e+00],
       [  2.16900000e+06,  -2.38630000e+06,  -1.33000000e+02,
          2.00000000e+01]])

Create a sympy Matrix object from the numpy array:

In [156]: import sympy

In [157]: M = sympy.Matrix(m)

In [158]: M
Out[158]: 
Matrix([
[    10,  11,   0,  0],
[     9,  10,   0,  0],
[   100,   0,  20,  3],
[-10000, 200, 133, 20]])

The sympy calculations are exact:

In [159]: M.det()
Out[159]: 1

In [160]: M.inv()
Out[160]: 
Matrix([
[     10,      -11,    0,  0],
[     -9,       10,    0,  0],
[-325400,   358000,   20, -3],
[2169000, -2386300, -133, 20]])

To convert the inverse back to a numpy array, you can do:

In [185]: Minv = M.inv()

In [186]: minv = np.asarray(Minv).astype(int)

In [187]: minv
Out[187]: 
array([[      10,      -11,        0,        0],
       [      -9,       10,        0,        0],
       [ -325400,   358000,       20,       -3],
       [ 2169000, -2386300,     -133,       20]])

It was necessary to use the astype(int) method, because the numpy array has dtype object without it:

In [188]: np.asarray(Minv)
Out[188]: 
array([[10, -11, 0, 0],
       [-9, 10, 0, 0],
       [-325400, 358000, 20, -3],
       [2169000, -2386300, -133, 20]], dtype=object)

That result is a numpy array of sympy integers.

Be careful with the conversion back to numpy: the integers in the inverse might be larger than can be represented with a 64 bit integer.