В 64-битных системах существует тип numpy.float128 d. (Я полагаю, что в 32-битных системах также существует тип float96 d). Хотя 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 впоследствии и отмасштабировать.
1026 * Е.Г. *
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. Однако это единственное решение, которое на самом деле выигрывает от дополнительной точности!