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 があります。( float9632 ビット システムにも dtype があると思います) numpy.linalg.eig128 ビットの float をサポートしていませんが、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 を使用したくない場合の回避策の 1 つは、べき乗の後、固有値を解く前に再スケーリングし、それらを「通常の」浮動小数点としてキャストし、固有値を解決してから、後で 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 の代わりに 1.0000 を使用します。この問題が発生した通常の python スクリプトでは、これが役に立ちました。問題が numpy の異常によって引き起こされ、python とは関係がない場合を除き、これは役立つはずです。
幸運を！