谱分解重构矩阵与原矩阵不符的问题排查

原代码 #

import math
import numpy as np

def gauss_n(n):
    # Create tri-diagonal matrix
    A = np.zeros((n, n))
    for i in range(1, n):
        curr_val = i/math.sqrt((4*(i**2))-1)  # Compute values
        A[i][i-1] = curr_val  # Place value
        A[i-1][i] = curr_val  # Place value symmetrically-opposite

    eig_values, eig_vectors = np.linalg.eig(A)  # Compute eigenvalues and eigenvectors

    # Calculate: A = ZΛZ^T  - Spectral Decomposition (Z in this case is "eig_vectors")
    dia_eig = np.diag(eig_values)  # Create diagonal-matrix with eigenvalues (Λ)

    # Because our initial matrix was symmetrical, therefore it means that the set of eigenvectors are
    # all orthogonal to each other
    A2 = np.matmul(np.matmul(eig_vectors, dia_eig), np.linalg.inv(eig_vectors))  # The Z^T = Z^-1

疑问解答 #

你完全有理由期望A2和A相等，代码的问题出在两个关键点：

1. 正交矩阵的逆矩阵计算方式错误 #

实对称矩阵的特征向量矩阵是正交矩阵，满足转置等于逆矩阵（即Zᵀ = Z⁻¹），但直接调用np.linalg.inv()会引入额外的数值计算误差，正确的做法是直接使用特征向量矩阵的转置eig_vectors.T代替逆矩阵。

2. 选错了特征分解函数 #

np.linalg.eig是通用特征分解函数，对对称矩阵的处理没有针对性；而np.linalg.eigh是专门为Hermitian矩阵（实对称矩阵属于此类）设计的，它返回的特征向量是严格正交归一化的，数值稳定性更强，能大幅降低重构时的误差。

修正后的代码片段 #

把特征分解和矩阵重构部分替换为：

# 使用专门针对对称矩阵的eigh函数
eig_values, eig_vectors = np.linalg.eigh(A)
dia_eig = np.diag(eig_values)
# 用转置代替逆矩阵，简化矩阵乘法写法
A2 = eig_vectors @ dia_eig @ eig_vectors.T

补充说明 #

由于浮点数计算的固有精度限制，修正后的A和A2不会完全逐位相等，但两者的元素差异会极小（通常在1e-10量级），可以用np.allclose(A, A2)来验证两者是否近似相等。

内容的提问来源于stack exchange，提问作者Dylan Streicher