我目前正在尝试寻找条件差的实对称矩阵的特征值，其元素可能非常小/很大（从1e-8到1e15），并且可能为负。矩阵虽然很小（4x4），所以执行速度不是主要问题。
对角化之前的所有操作都使用“ logsumexp”技巧的变体完成（适应于此线程numerically stable way to multiply log probability matrices in numpy的矩阵-矩阵/矩阵矢量乘法，并进一步适应负系数），所以我最终得到两个矩阵（ sgnM，logM），分别包含M中元素的绝对值的符号和对数。这部分效果很好。
但是，我没有找到任何等效于eigvalsh的文档来将其作为输入，并一直使用这些数字技巧直到它返回特征值。现在，我只使用
scipy.linalg.eigvalsh(sgnM*np.exp(logM))
仍然会遇到精度问题（当我更改参数时，我感兴趣的特征值从1e-4变为1e-9左右，我可以看到1e-9周围的特征值比第一个嘈杂。直到从那里得到的结果不再有意义）。
是否有某种功能可以使我寻找的东西隐藏在现有的线性代数引擎中？如果不是，我是否仍然可以依靠Lapack / MKL / Blas例程来实现它，还是需要从头开始？
我设法实现了它，很不幸，它对我的情况没有太大的改善。万一有人想尝试另一个问题（基于https://github.com/mateuv/MetodosNumericos/blob/master/python/NumericalMethodsInEngineeringWithPython/jacobi.py中的代码，我仍将其保留在此处）
import numpy as np
from scipy.special import logsumexp
def eigsh_jacobi_stable(log_a, sgn_a, tol = 1.0e-14): # Jacobi method
def maxElem(log_a): # Find largest off-diag. element a[k,l]
n = len(log_a)
log_aMax = np.log(1e-36)
for i in range(n-1):
for j in range(i+1,n):
if log_a[i,j] >= log_aMax:
log_aMax = log_a[i,j]
k = i; l = j
return log_aMax,k,l
def rotate(log_a, sgn_a, log_p, sgn_p, k, l): # Rotate to make a[k,l] = 0
n = len(log_a)
log_aDiff, sgn_aDiff = logsumexp(a=[log_a[l,l], log_a[k,k]],
b=[sgn_a[l,l], -sgn_a[k,k]],
return_sign=True)
if log_a[k,l] - log_aDiff < np.log(1.0e-36):
log_t = log_a[k,l] - log_aDiff
sgn_t = sgn_a[k,l] * sgn_aDiff
else:
log_phi = log_aDiff - np.log(2) - log_a[k,l]
sgn_phi = sgn_aDiff * sgn_a[k,l]
log_t = -logsumexp(a=[log_phi, .5*logsumexp(a=[2*log_phi, 0])])
sgn_t = 1.
if sgn_phi < 0:
sgn_t = -1
log_c = - .5 * logsumexp(a=[2*log_t, 0])
sgn_c = 1.
log_s = log_t + log_c
sgn_s = sgn_t * sgn_c
log_tau = log_s - logsumexp(a=[log_c, 0])
sgn_tau = sgn_s
log_temp, sgn_temp = log_a[k,l], sgn_a[k,l]
log_a[k,l] = np.log(1.0e-36)
sgn_a[k, l] = 0
log_a[k,k], sgn_a[k,k] = logsumexp(a=[log_a[k,k], log_t + log_temp],
b=[sgn_a[k,k], -sgn_t * sgn_temp],
return_sign=True)
log_a[l,l], sgn_a[l,l] = logsumexp(a=[log_a[l,l], log_t + log_temp],
b=[sgn_a[l,l], sgn_t * sgn_temp],
return_sign=True)
for i in range(k): # Case of i < k
log_temp, sgn_temp = log_a[i,k], sgn_a[i,k]
log_plop, sgn_plop = logsumexp(a=[log_a[i, l], log_tau + log_temp],
b=[sgn_a[i, l], sgn_tau * sgn_temp],
return_sign=True)
log_a[i,k], sgn_a[i,k] = logsumexp(a=[log_temp, log_s + log_plop],
b=[sgn_temp, -sgn_plop * sgn_s],
return_sign=True)
log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_a[i, l]],
b=[sgn_temp, -sgn_tau * sgn_a[i, l]],
return_sign=True)
log_a[i,l], sgn_a[i,l] = logsumexp(a=[log_a[i, l], log_s + log_plop],
b=[sgn_a[i, l], sgn_plop * sgn_s],
return_sign=True)
for i in range(k+1,l): # Case of k < i < l
log_temp, sgn_temp = log_a[k, i], sgn_a[k,i]
log_plop, sgn_plop = logsumexp(a=[log_a[i, l], log_tau + log_a[k,i]],
b=[sgn_a[i, l], sgn_tau * sgn_a[k,i]],
return_sign=True)
log_a[k, i], sgn_a[k, i] = logsumexp(a=[log_temp, log_s + log_plop],
b=[sgn_temp, -sgn_plop * sgn_s],
return_sign=True)
log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_a[i, l]],
b=[sgn_temp, -sgn_tau * sgn_a[i, l]],
return_sign=True)
log_a[i,l], sgn_a[i,l] = logsumexp(a=[log_a[i, l], log_s + log_plop],
b=[sgn_a[i, l], sgn_plop * sgn_s],
return_sign=True)
for i in range(l+1,n): # Case of i > l
log_temp, sgn_temp = log_a[k, i], sgn_a[k,i]
log_plop, sgn_plop = logsumexp(a=[log_a[l,i], log_tau + log_temp],
b=[sgn_a[l,i], sgn_tau * sgn_temp],
return_sign=True)
log_a[k, i], sgn_a[k, i] = logsumexp(a=[log_temp, log_s + log_plop],
b=[sgn_temp, -sgn_plop * sgn_s],
return_sign=True)
log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_a[l,i]],
b=[sgn_temp, -sgn_tau * sgn_a[l,i]],
return_sign=True)
log_a[l,i], sgn_a[l,i] = logsumexp(a=[log_a[l,i], log_s + log_plop],
b=[sgn_a[l, i], sgn_plop * sgn_s],
return_sign=True)
for i in range(n): # Update transformation matrix
log_temp, sgn_temp = log_p[i,k], sgn_p[i,k]
log_plop, sgn_plop = logsumexp(a=[log_p[i,l], log_tau + log_p[i,k]],
b=[sgn_p[i,l], sgn_tau * sgn_p[i,k]],
return_sign=True)
log_p[k, i], sgn_p[k, i] = logsumexp(a=[log_temp, log_s + log_plop],
b=[sgn_temp, -sgn_plop * sgn_s],
return_sign=True)
log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_p[i, l]],
b=[sgn_temp, -sgn_tau * sgn_p[i, l]],
return_sign=True)
log_p[i,l], sgn_p[i,l] = logsumexp(a=[log_p[i, l], log_s + log_plop],
b=[sgn_p[i, l], sgn_plop * sgn_s],
return_sign=True)
n = len(log_a)
maxRot = 5*(n**2) # Set limit on number of rotations
log_p, sgn_p = np.log(1e-36) * np.ones_like(log_a), np.zeros_like(log_a)
np.fill_diagonal(log_p, 0)
np.fill_diagonal(sgn_p, 1)
for i in range(maxRot): # Jacobi rotation loop
log_aMax,k,l = maxElem(log_a)
# print(log_aMax)
if np.exp(log_aMax) < tol:
return np.sort(np.exp(np.diag(log_a)) * np.diag(sgn_a))
rotate(log_a, sgn_a, log_p, sgn_p ,k,l)
raise RuntimeError('Jacobi method did not converge')
if __name__ == '__main__':
from scipy.linalg import eigvalsh
for _ in range(100):
test = np.random.randn(10, 10)
test += test.transpose()
log_test, sgn_test = np.log(np.abs(test)), np.sign(test)
eigs_jacobi = np.sort(eigsh_jacobi_stable(log_test, sgn_test)) # this sort is unnecessary, already done internally
eigs_scipy = np.sort(eigvalsh(test))
print(eigs_jacobi-eigs_scipy)