numpy 任意精度线性代数

Question

我有一个大小为500x500的numpy二维数组。我想找出这个数组每个元素的指数值（也就是e的幂）。问题是，有些值非常负，比如-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来解决这个问题，但没有成功。欢迎任何建议。

Answer 1

据我所知，numpy不支持比双精度（float64）更高的精度，如果不特别说明，默认就是双精度。

你可以试试这个：http://code.google.com/p/mpmath/

这个工具的功能列表（还有其他功能）

算术运算：

• 可以处理任意精度的实数和复数
• 指数的大小没有限制

Answer 2

在64位系统上，有一种叫做 numpy.float128 的数据类型。（我记得在32位系统上也有 float96 的数据类型。）虽然 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的矩阵。不过，它是唯一一个真正能从额外精度中受益的解法！

Answer 3

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，结果计算特征向量的时候就失败了。