我有一个 numpy 2d 数组 [中/大型 - 比如说 500x500]。我想找到它的元素指数的特征值。问题是某些值是非常负的（-800、-1000 等），并且它们的指数下溢（意味着它们非常接近于零，因此 numpy 将它们视为零）。无论如何在numpy中使用任意精度？
我梦想的方式：
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)
我用 gmpy 和 mpmath 寻找解决方案无济于事。任何想法都会受到欢迎。
SymPy 可以计算任意精度：
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
输出：
-2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807463648841185335e-261
2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807466621962539464e-261
[[-0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999994391176386872]
[ 1]]
[[1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000560882361313]
[ 1]]
您可以将 N(x, 100) 中的 100 更改为其他精度，但是，当我尝试 1000 时，特征向量的计算失败。
在 64 位系统上，有一个numpy.float128dtype。float96（我相信32 位系统上也有dtype）虽然numpy.linalg.eig不支持 128 位浮点数，但scipy.linalg.eig（有点）支持。
但是，从长远来看，这一切都无关紧要。任何针对特征值问题的通用求解器都将是迭代的，而不是精确的，因此通过保持额外的精度你不会获得任何东西！ np.linalg.eig适用于任何形状，但从不返回精确解。
如果您总是在求解 2x2 矩阵，那么编写自己的求解器就很简单了，它应该更精确。我将在最后展示一个例子......
无论如何，向前推进毫无意义的精确内存容器：
import numpy as np
import scipy as sp
import scipy.linalg
a = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)
ex = np.exp(a)
print ex
eigvals, eigvecs = sp.linalg.eig(ex)
# And to test...
check1 = ex.dot(eigvecs[:,0])
check2 = eigvals[0] * eigvecs[:,0]
print 'Checking accuracy..'
print check1, check2
print (check1 - check2).dot(check1 - check2), '<-- Should be zero'
然而，你会注意到你得到的和做的一样np.linalg.eig(ex.astype(np.float64)。事实上，我很确定这就是scipy正在做的事情，同时numpy引发错误而不是默默地做。我可能完全错了，虽然...
如果您不想使用 scipy，一种解决方法是在求幂之后但在求解特征值之前重新缩放事物，将它们转换为“正常”浮点数，求解特征值，然后将事物重新转换为 float128 并重新缩放。
例如
import numpy as np
a = np.array([[-800.21,-600.00],[-600.00,-1000.48]], dtype=np.float128)
ex = np.exp(a)
factor = 1e300
ex_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, check2
print (check1 - check2).dot(check1 - check2), '<-- Should be zero'
最后，如果您只求解 2x2 或 3x3 矩阵，您可以编写自己的求解器，它将返回这些矩阵形状的精确值。
import numpy as np
def 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, root2
def 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 vals
def 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 vecs
def eig_2x2(matrix_2x2):
vals = eigvals(matrix_2x2)
vecs = eigvecs(matrix_2x2, vals)
return vals, vecs
a = 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, check2
print (check1 - check2).dot(check1 - check2), '<-- Should be zero'
这个返回一个真正精确的解决方案，但仅适用于 2x2 矩阵。然而，这是唯一真正受益于额外精度的解决方案！
据我所知，numpy 不支持高于双精度（float64），如果没有指定，这是默认的。
尝试使用这个：http ://code.google.com/p/mpmath/
功能列表（除其他外）
算术：
所以你需要 350 位的精度。使用 IEEE 浮点数（numpy 使用的是什么）你不会得到它。您可以使用 bc 程序获得它：
$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.
scale=350
e(-800)
.<hundreds of zeros>00366
我特别没有使用 numpy 的经验，但是添加一个带有一定数量的零的小数点可能会有所帮助。例如，使用 1.0000 而不是 1。在我遇到此问题的普通 python 脚本中，这会有所帮助，因此除非您的问题是由 numpy 中的怪异引起的并且与 python 无关，否则这应该会有所帮助。
祝你好运！