Numpy矩阵逆计算精度异常：求解M⁻¹K特征值遇难题

问题背景 #

计算矩阵乘积M⁻¹K的特征值时，通过np.linalg.inv(M)得到的逆矩阵M_inv与M相乘后，结果偏离单位矩阵，导致无法得到有效特征值。矩阵M和K的初始化代码如下：

import numpy as np

def mij(i, j, h):
    if i==j:
        return 2.0 * h / 3.0
    else:
        return h / 6.0

def kij(i, j, h):
    if i==j:
        return 2.0 / h
    else:
        return -1 / h

n = 500
h = 1 / n

M = np.zeros((n, n))
K = np.zeros((n, n))

for i in range(0, n):
    M[i,i] = mij(i, i, h)
    if i+1 < n:
        M[i,i+1] = mij(i, i+1, h)
    if i-1 >= 0:
        M[i,i-1] = mij(i, i-1, h)

    K[i,i] = kij(i, i, h)
    if i+1 < n:
        K[i,i+1] = kij(i, i+1, h)
    if i-1 >= 0:
        K[i,i-1] = kij(i, i-1, h)

解决方案 #

1. 直接求解广义特征值（最优方案） #

无需显式计算M的逆矩阵，直接求解广义特征值问题Kx = λMx（等价于M⁻¹Kx = λx），使用np.linalg.eigvals的广义特征值模式：

# 求解M⁻¹K的特征值，等价于广义特征值问题Kx=λMx
eigenvalues = np.linalg.eigvals(K, M)

该方法避免了矩阵求逆带来的精度损失，计算效率和稳定性远优于显式求逆后再计算特征值。

2. 使用更高精度的数据类型 #

将矩阵M和K转换为更高精度的浮点类型（如np.float128）后再求逆，减少精度损失：

# 转换为高精度类型
M_high = M.astype(np.float128)
K_high = K.astype(np.float128)

# 计算逆矩阵
M_inv_high = np.linalg.inv(M_high)

# 验证逆矩阵
with np.printoptions(threshold=np.inf, precision=10, suppress=True, linewidth=20000):
    print(np.matrix(M_inv_high @ M_high))  # 使用@运算符做矩阵乘法，比*更可靠

# 计算特征值
eigenvalues_high = np.linalg.eigvals(M_inv_high @ K_high)

3. 利用稀疏矩阵特性 #

M和K都是三对角稀疏矩阵，使用scipy.sparse库处理稀疏矩阵，既能节省内存，也能提升计算精度和速度：

from scipy.sparse import diags
from scipy.sparse.linalg import eigs

# 构造稀疏矩阵M
m_diag = np.full(n, 2.0 * h / 3.0)
m_off_diag = np.full(n-1, h / 6.0)
M_sparse = diags([m_off_diag, m_diag, m_off_diag], [-1, 0, 1], shape=(n, n))

# 构造稀疏矩阵K
k_diag = np.full(n, 2.0 / h)
k_off_diag = np.full(n-1, -1 / h)
K_sparse = diags([k_off_diag, k_diag, k_off_diag], [-1, 0, 1], shape=(n, n))

# 求解广义特征值（可指定需要的特征值数量，比如前10个）
eigenvalues_sparse, _ = eigs(K_sparse, M=M_sparse, k=10, which='LM')

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