numpy: maintain decimal class through inverse calculation

Question

I'm trying to find the inverse of a matrix made up of a specific class (decimal.Decimal) and keep the values as Decimal objects throughout the process (to preserve exactness throughout the calculation).

My problem is numpy.linalg.inverse always returns the matrix values as floats. I've figured out a work around by changing the type from floats to Decimal objects after the inverse is calculated but I'd prefer to maintain the class of the original matrix throughout (I'm worried I may be losing accuracy when the numbers are converted to floats)?

So I guess I have a few questions: (1) am I losing accuracy/exactness when the values of the matrix are converted to float types (I'm dealing with an 11 by 11 matrix); if so, (2) is there anyway to keep the values as decimal.Decimal throughout the calculation using numpy; if not, (3) is there another module / method I should consider for this type of calculation?

Here's an example of what my code will look like:

import numpy as np
from decimal import Decimal as D

a = np.array( [ [D('1'), D('2'), D('3')],
                [D('4'), D('5'), D('6')],
                [D('7'), D('8'), D('9')] ] )

print type(a[0,0])
# <class 'decimal.Decimal'>

inverse_a = np.linalg.inv(a)

print type(inverse_a[0,0])
# <type 'numpy.float64'>

inverse_a_Decimal_flat = [D(str(i)) for i in inverse_a.flat]  # change all values to decimal.Decimal
inverse_a_Decimal = np.reshape(inverse_a_Decimal_flat, [3, 3])  # reshape to 3x3

print type(inverse_a_Decimal[0,0]), d.shape
# <class 'decimal.Decimal'> (3,3)

1. Yes. 2. Not that I'm aware of. 3. Have you tried writing it yourself? — jonrsharpe, Dec 08 '16 at 18:20
So if it's exactness I'm after it looks like I may need to use the `np.longdouble` type [link](http://stackoverflow.com/questions/25181642/how-set-numpy-floating-point-accuracy?noredirect=1&lq=1). How does this differ from `decimal` in terms of exactness? And it sounds like it will take longer to compute so I'm wondering if it's even worth the trouble. `float64` may be good enough for one step in the calculation (I'm calculating the chemical composition of a mixture). — Johnny Metz, Dec 08 '16 at 18:41
If you want exactness, then `Decimal` is not enough - you need `fractions.Fraction`, else you'll end up rounding to some (negative) power of 10. You might want to look into using [`sympy`](http://docs.sympy.org/0.7.2/modules/matrices/matrices.html) for infinite-precision matrix algebra — Eric, Dec 08 '16 at 23:03
For a 3x3 matrix, there is an analytical expression for the inverse. It's something like the 9 2x2 determinants divided by the full determinant. Check a symbolic math program. — hpaulj, Dec 09 '16 at 02:56

score 2 · Answer 1 · answered Dec 08 '16 at 23:11

Is using sympy.Matrix acceptable?

If we use that:

>>> import sympy as sy
>>> a = sy.Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> a.inv()
ValueError: Matrix det == 0; not invertible.

We realize that inverting this matrix doesn't even make sense!

Trying again swapping the 9 with a 10, gives

Matrix([
[-2/3, -4/3,  1],
[-2/3, 11/3, -2],
[   1,   -2,  1]])

Which is the desired infinite-precision result.

hpaulj · Accepted Answer · 2016-12-08T19:12:51.523

Your a is a object array:

In [255]: a
Out[255]: 
array([[Decimal('1'), Decimal('2'), Decimal('3')],
       [Decimal('4'), Decimal('5'), Decimal('6')],
       [Decimal('7'), Decimal('8'), Decimal('9')]], dtype=object)

In my newly updated numpy I get an error when trying to do an inverse.

In [257]: np.linalg.inv(a)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-257-5c2063aa5f16> in <module>()
----> 1 np.linalg.inv(a)

/usr/local/lib/python3.5/dist-packages/numpy/linalg/linalg.py in inv(a)
    524     signature = 'D->D' if isComplexType(t) else 'd->d'
    525     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 526     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    527     return wrap(ainv.astype(result_t, copy=False))
    528 

TypeError: No loop matching the specified signature and casting
was found for ufunc inv
In [258]: np.__version__
Out[258]: '1.11.2'

inv is probably using compiled code (_umath_linalg.inv), which will be coded to work with standard c floats (float or double).

You asked about longdouble:

In [266]: a1=np.arange(1,10,dtype=np.longdouble).reshape(3,3)
In [267]: a1
Out[267]: 
array([[ 1.0,  2.0,  3.0],
       [ 4.0,  5.0,  6.0],
       [ 7.0,  8.0,  9.0]], dtype=float96)
In [268]: np.linalg.inv(a1)
...
TypeError: array type float96 is unsupported in linalg

Again, it's a question of how linalg was compiled.

In [269]: a1=np.arange(1,10,dtype=np.double).reshape(3,3)
In [270]: np.linalg.inv(a1)
Out[270]: 
array([[  3.15221191e+15,  -6.30442381e+15,   3.15221191e+15],
       [ -6.30442381e+15,   1.26088476e+16,  -6.30442381e+15],
       [  3.15221191e+15,  -6.30442381e+15,   3.15221191e+15]])

I tried to dot this and didn't get the expected np.eye array. But this is a nearly singular array (determinant is 0, or nearly so).

In [274]: np.linalg.det(a1)
Out[274]: -9.5171266700777698e-16

So it's a bad example.

Shuffling values to the array is no longer singular

In [288]: a1=np.array([[3,1,2],[4,6,5],[8,9,7]],np.float)
In [289]: a1.dot(np.linalg.inv(a1))
Out[289]: 
array([[  1.00000000e+00,   0.00000000e+00,   2.22044605e-16],
       [ -4.44089210e-16,   1.00000000e+00,   4.44089210e-16],
       [  0.00000000e+00,   0.00000000e+00,   1.00000000e+00]])

I get the same thing if a1 is np.int, since the np.int is upcast to float64. I can test with np.float32 to see the effect of lower precision. As noted, I can't go to longfloat.

inv often is calculated as a I/a1 linear equation solution:

In [297]: np.linalg.solve(a1,np.eye(3))
Out[297]: 
array([[ 0.14285714, -0.52380952,  0.33333333],
       [-0.57142857, -0.23809524,  0.33333333],
       [ 0.57142857,  0.9047619 , -0.66666667]])
In [298]: np.linalg.inv(a1)
Out[298]: 
array([[ 0.14285714, -0.52380952,  0.33333333],
       [-0.57142857, -0.23809524,  0.33333333],
       [ 0.57142857,  0.9047619 , -0.66666667]])

I'm using version 1.9.3 -- I'll update and try again. So if this is the case would you recommend converting to `float64` or a different data type for exactness (maybe `np.longdouble`)? — Johnny Metz, Dec 08 '16 at 18:44
Any idea what the exactness difference is between `float64` and `np.longdouble` (or `float96` in this case - mine evaluates to `float128`)? Like how many points after the decimal are accurate in each case — Johnny Metz, Dec 08 '16 at 18:54
Why are you worried about precision? Why is important in this problem? Often taking a standalong `inv` is discouraged. Instead we are told to use `np.linalg.solve`. In fact `inv` might actually be implemented as `solve(a,I)`. — hpaulj, Dec 08 '16 at 19:10
This is one of the steps in calculating the chemical composition of a gaseous solution. I would like to be as precise as possible. Sounds like `float64` should be good enough for this step though (and using `decimal` the rest of the way). `np.longdouble` looks like a headache, may not be supported on some machines and is easy to lose ([Extended Precision at bottom](https://docs.scipy.org/doc/numpy/user/basics.types.html)) — Johnny Metz, Dec 08 '16 at 19:19
@hpaulj: In particular, `inv` has failed on the original input because the matrix is singular, but `solve` may not — Eric, Dec 08 '16 at 23:11
The singularity of the original array creates these outlandish `e15` values in the inverse. It's a bad example to test precision on. — hpaulj, Dec 08 '16 at 23:41

numpy: maintain decimal class through inverse calculation

2 Answers2