Search code examples
pythonmathnumpylinear-algebraarbitrary-precision

numpy arbitrary precision linear algebra


I have a numpy 2d array [medium/large sized - say 500x500]. I want to find the eigenvalues of the element-wise exponent of it. The problem is that some of the values are quite negative (-800,-1000, etc), and their exponents underflow (meaning they are so close to zero, so that numpy treats them as zero). Is there anyway to use arbitrary precision in numpy?

The way I dream it:

import numpy as np

np.set_precision('arbitrary') # <--- Missing part
a = np.array([[-800.21,-600.00],[-600.00,-1000.48]])
ex = np.exp(a)  ## Currently warns about underflow
eigvals, eigvecs = np.linalg.eig(ex)

I have searched for a solution with gmpy and mpmath to no avail. Any idea will be welcome.


Solution

  • SymPy can calculate arbitrary precision:

    from sympy import exp, N, S
    from sympy.matrices import Matrix
    
    data = [[S("-800.21"),S("-600.00")],[S("-600.00"),S("-1000.48")]]
    m = Matrix(data)
    ex = m.applyfunc(exp).applyfunc(lambda x:N(x, 100))
    vecs = ex.eigenvects()
    print vecs[0][0] # eigen value
    print vecs[1][0] # eigen value
    print vecs[0][2] # eigen vect
    print vecs[1][2] # eigen vect
    

    output:

    -2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807463648841185335e-261
    2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807466621962539464e-261
    [[-0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999994391176386872]
    [                                                                                                      1]]
    [[1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000560882361313]
    [                                                                                                    1]]
    

    you can change 100 in N(x, 100) to other precision, but, as I tried 1000, the calculation of eigen vect failed.