Compute matrix inverse with decimal object

Question

There's a related questions Matrix inverse with Decimal type NumPy 2015 a while ago which did not have a definite answer. There's a second question from me Is there a way for python to perform a matrix inversion at 500 decimal precision where hpaulj provided some updated suggestions.

Basically decimal is a standard python library capable of computing arbitrary precession value at arbitrary orders. It can be operated by most of the Numpy function such as poly to evaluate the polynomials

np.polyval([Decimal(1),Decimal(2)], Decimal(3.1) )

Decimal('5.100000000')

It can also be cased to a numpy array or being initiated as a dtype object array(Are Decimal 'dtypes' available in NumPy? 2011).

np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])

array([[Decimal('1'), Decimal('2')],
       [Decimal('3'), Decimal('4')]], dtype=object)

matrix_m=np.zeros((2,2) ,dtype=np.dtype)
for ix in range(0,2):
    for iy in range(0,2):
        matrix_m[ix,iy]=Decimal(ix)+Decimal(iy);

array([[Decimal('0'), Decimal('1')],
       [Decimal('1'), Decimal('2')]], dtype=object)

Some array operation from numpy also worked when Decimal was the element,

np.exp( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]]) )

array([[Decimal('2.718281828'), Decimal('7.389056099')],
       [Decimal('20.08553692'), Decimal('54.59815003')]], dtype=object)

np.sqrt( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]]) )

array([[Decimal('1'), Decimal('1.414213562')],
       [Decimal('1.732050808'), Decimal('2')]], dtype=object)

and, at single element, the numpy calculation agreed with decimal's native function

np.exp(Decimal(1))==Decimal(1).exp()
True

The useful constant was also provided

def pi():
    """Compute Pi to the current precision.
    #https://docs.python.org/3/library/decimal.html

    >>> print(pi())
    3.141592653589793238462643383

    """
    getcontext().prec += 2  # extra digits for intermediate steps
    three = Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n+na, na+8
        d, da = d+da, da+32
        t = (t * n) / d
        s += t
    getcontext().prec -= 2
    return +s               # unary plus applies the new precision

However, it turned out that both the determinate and the inverse of the matrix in the numpy

np.linalg.det(np.array([[Decimal(1),Decimal(2)],[Decimal(1),Decimal(3)]]))

File <__array_function__ internals>:180, in det(*args, **kwargs)
File ~\anaconda3\lib\site-packages\numpy\linalg\linalg.py:2154, in det(a)
   2152 t, result_t = _commonType(a)
   2153 signature = 'D->D' if isComplexType(t) else 'd->d'
-> 2154 r = _umath_linalg.det(a, signature=signature)
   2155 r = r.astype(result_t, copy=False)
   2156 return r

np.linalg.inv(np.array([[Decimal(1),Decimal(2)],[Decimal(1),Decimal(3)]]))

File <__array_function__ internals>:180, in inv(*args, **kwargs)
File ~\anaconda3\lib\site-packages\numpy\linalg\linalg.py:552, in inv(a)
    550 signature = 'D->D' if isComplexType(t) else 'd->d'
    551 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 552 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    553 return wrap(ainv.astype(result_t, copy=False))

returned the same error

UFuncTypeError: Cannot cast ufunc 'inv' input from dtype('O') to dtype('float64') with casting rule 'same_kind'

Which is not what was intended. It should just calculate the object according to the arithmetic and decimal should be able to compute the value itself. hpaulj's post provided an alternative method to cast the decimal object to mfp object of mpmath package

mp.matrix( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])) 

matrix(
[['1.0', '2.0'],
 ['3.0', '4.0']])

mp.matrix( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])) [0,0]

mpf('1.0')

and then perform the inverse in mpmath package.

mp.matrix( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])) **(-1)

matrix(
[['-2.0', '1.0'],
 ['1.5', '-0.5']])

This could work, however, it lost the nice functionally of decimal package and involved large amount casting elements from mpmath to numpy and decimal objects. The mpf() object's computational speed is also significantly slower than the calculation Decimal() object's.

Is there an easy way to write or improve the code from the numpy package directly so that a np.inverse() could be used on decimal array? Is there any way to compute the matrix inverse with decimal object?

The linked posts should make it clear that the `numpy` functions are compiled code, and only work with standard `c` types (float and double). They do not have a slow `object` dtype fall back option. `mpmath` has advanced math and linalg functions, hence my previous answer. `decimal` has none of that (little beyond `exp` and `log`). `sympy` can give you a horribly long analytical expression. — hpaulj, Mar 26 '23 at 16:35
@hpaulj It turned out with numpy.ndarray.tolist() and https://stackoverflow.com/a/61741074/11141816 the decimal package can work with numpy and compute most of the linear algebra functions, much faster and accurate than sympy or mpmath's calculation. The numerical stability was an issue though. It's curious why no one updated the packages to make the compatibility official. I tried it out in a complicated calculation, the decimal package combined with numpy was almost 1 million times faster than sympy with the same intermediate accuracy. — ShoutOutAndCalculate, Mar 27 '23 at 15:01
As long as the array shape is small (not much bigger than your sample (2,2)), this pure-python approach makes sense. I remember doing that kind of thing before-computers. `numpy` users and coders tend to think in terms of what's best for (1000,1000) sizes. — hpaulj, Mar 27 '23 at 15:57

Compute matrix inverse with decimal object

0 Answers0