numpy arbitrary precision linear algebra

SymPy can calculate arbitrary precision:

from sympy import exp, N, Sfrom sympy.matrices import Matrixdata = [[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 valueprint vecs[1][0] # eigen valueprint vecs[0][2] # eigen vectprint vecs[1][2] # eigen vect

output:

-2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807463648841185335e-2612.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.

On 64-bit systems, there's a numpy.float128 dtype. (I believe there's a float96 dtype on 32-bit systems, as well) While numpy.linalg.eig doesn't support 128-bit floats, scipy.linalg.eig (sort of) does.

However, none of this is going to matter, in the long run. Any general solver for an eigenvalue problem is going to be iterative, rather than exact, so you're not gaining anything by keeping the extra precision! np.linalg.eig works for any shape, but never returns an exact solution.

If you're always solving 2x2 matrices, it's trivial to write your own solver that should be more exact. I'll show an example of this at the end...

Regardless, forging ahead into pointlessly precise memory containers:

import numpy as npimport scipy as spimport scipy.linalga = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)ex = np.exp(a)print exeigvals, eigvecs = sp.linalg.eig(ex)# And to test...check1 = ex.dot(eigvecs[:,0])check2 = eigvals[0] * eigvecs[:,0]print 'Checking accuracy..'print check1, check2print (check1 - check2).dot(check1 - check2), '<-- Should be zero'

However, you'll notice that what you get is identical to just doing np.linalg.eig(ex.astype(np.float64). In fact, I'm fairly sure that's what scipy is doing, while numpy raises an error rather than doing it silently. I could be quite wrong, though...

If you don't want to use scipy, one workaround is to rescale things after the exponentiation but before solving for the eigenvalues, cast them as "normal" floats, solve for the eigenvalues, and then recast things as float128's afterwards and rescale.

E.g.

import numpy as npa = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)ex = np.exp(a)factor = 1e300ex_rescaled = (ex * factor).astype(np.float64)eigvals, eigvecs = np.linalg.eig(ex_rescaled)eigvals = eigvals.astype(np.float128) / factor# And to test...check1 = ex.dot(eigvecs[:,0])check2 = eigvals[0] * eigvecs[:,0]print 'Checking accuracy..'print check1, check2print (check1 - check2).dot(check1 - check2), '<-- Should be zero'

Finally, if you're only solving 2x2 or 3x3 matrices, you can write your own solver, which will return an exact value for those shapes of matrices.

import numpy as npdef quadratic(a,b,c):    sqrt_part = np.lib.scimath.sqrt(b**2 - 4*a*c)    root1 = (-b + sqrt_part) / (2 * a)    root2 = (-b - sqrt_part) / (2 * a)    return root1, root2def eigvals(matrix_2x2):    vals = np.zeros(2, dtype=matrix_2x2.dtype)    a,b,c,d = matrix_2x2.flatten()    vals[:] = quadratic(1.0, -(a+d), (a*d-b*c))    return valsdef eigvecs(matrix_2x2, vals):    a,b,c,d = matrix_2x2.flatten()    vecs = np.zeros_like(matrix_2x2)    if (b == 0.0) and (c == 0.0):        vecs[0,0], vecs[1,1] = 1.0, 1.0    elif c != 0.0:        vecs[0,:] = vals - d        vecs[1,:] = c    elif b != 0:        vecs[0,:] = b        vecs[1,:] = vals - a    return vecsdef eig_2x2(matrix_2x2):    vals = eigvals(matrix_2x2)    vecs = eigvecs(matrix_2x2, vals)    return vals, vecsa = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)ex = np.exp(a)eigvals, eigvecs =  eig_2x2(ex) # And to test...check1 = ex.dot(eigvecs[:,0])check2 = eigvals[0] * eigvecs[:,0]print 'Checking accuracy..'print check1, check2print (check1 - check2).dot(check1 - check2), '<-- Should be zero'

This one returns a truly exact solution, but will only work for 2x2 matrices. It's the only solution that actually benefits from the extra precision, however!

As far as I know, numpy does not support higher than double precision (float64), which is the default if not specified.

Try using this: http://code.google.com/p/mpmath/

List of features (among others)

Arithmetic:

• Real and complex numbers with arbitrary precision
• Unlimited exponent sizes / magnitudes