实对称矩阵用eigvalsh、eigvals、eigh计算特征值结果不一致

问题场景 #

计算实对称矩阵A的特征值时，np.linalg.eigvalsh、np.linalg.eigvals、np.linalg.eigh返回结果不一致，且矩阵行列式与特征值乘积不匹配。

代码实现 #

import numpy as np

def convert_Matrix(Var, dim):
    var_i = np.zeros((dim, dim))
    indices = np.tril_indices(dim)
    var_i[indices] = Var
    var = np.dot(var_i, var_i.T)
    return var

Var = np.array([4.04640469e+05, 9.98970587e+07, 1.75522655e+02, 1.85808145e+08,
            2.88277248e+08, 7.50545040e+05])
dim = 3

A= convert_Matrix(Var, dim)

det_A = np.linalg.det(A)
eigvalsh_values = np.linalg.eigvalsh(A)
eigvals_values = np.linalg.eigvals(A)
eigh_values = np.linalg.eigh(A)[0]  

print("A:\n", A)
print("Eigenvalues from eigvalsh:", eigvalsh_values)
print("Eigenvalues from eigvals:", eigvals_values)
print("Eigenvalues from eigh:", eigh_values)
print('det:', det_A)

输出结果 #

A:[[1.63733909e+11 4.04223927e+13 7.51854949e+13]
   [4.04223927e+13 9.97942234e+15 1.85617378e+16]
   [7.51854949e+13 1.85617378e+16 1.17629002e+17]]

Eigenvalues from eigvalsh:[2.43186146e-04 6.86886614e+15 1.20739722e+17]
Eigenvalues from eigvals: [1.20739722e+17 -5.60191504e-01  6.86886614e+15]
Eigenvalues from eigh:    [4.41546949e-04 6.86886614e+15 1.20739722e+17]

det: -1.4510608257526929e+28

差异原因 #

1. 函数算法的适配性差异 #

• eigvals：通用特征值求解函数，针对任意矩阵设计，不会利用实对称矩阵的特性。由于浮点数计算误差，会将本应接近0的极小特征值算出非零甚至负数（如输出中的-5.60191504e-01），属于数值误差导致的伪结果。
• eigvalsh和eigh：专门为实对称/厄米特矩阵优化的函数，利用矩阵对称性减少计算误差，返回的特征值理论上均非负（因为矩阵A是var_i * var_i.T，属于半正定矩阵）。两者的微小差异源于内部算法细节（如分治策略、精度控制），属于正常的数值精度范围波动。

2. 行列式与特征值乘积不匹配的根源 #

矩阵A的特征值存在极端数量级差异：最大特征值约1e17，极小特征值约1e-4，跨度达21个数量级。

• np.linalg.det基于LU分解计算，面对极端条件数的矩阵时，数值稳定性极差，会出现精度丢失甚至符号错误（半正定矩阵行列式应非负，但输出为负）。
• 直接计算特征值乘积时，极端数值的乘法会超出浮点数有效精度范围，导致结果偏差。

解决方案 #

1. 优先使用对称矩阵专用函数 #

对于实对称矩阵，完全放弃eigvals，选择eigvalsh或eigh：

• eigvalsh直接返回升序排列的特征值，适合仅需特征值的场景。
• eigh返回特征值与对应特征向量，适合需要特征向量的场景。

若需进一步统一两者结果，可通过指定driver参数对齐算法：

# 指定分治算法，与eigvalsh默认算法一致
eigh_values = np.linalg.eigh(A, driver='evd')[0]

2. 基于特征值正确计算行列式 #

避免直接调用np.linalg.det，改用对称矩阵特征值计算行列式：

# 利用eigvalsh结果计算，保证半正定矩阵乘积非负
det_from_eig = np.prod(eigvalsh_values)
print('det from eigvalsh:', det_from_eig)

该结果会符合半正定矩阵的行列式特性，且与特征值乘积一致。

3. 优化极端条件数矩阵的计算（可选） #

针对特征值跨度极大的矩阵，可先缩放矩阵元素，减少数值精度损失：

scale_factor = np.max(A)
A_scaled = A / scale_factor
# 计算缩放后矩阵的特征值，再还原为原矩阵特征值
eigval_scaled = np.linalg.eigvalsh(A_scaled)
eigval_original = eigval_scaled * scale_factor

内容的提问来源于stack exchange，提问作者skye W