提升numpy.linalg运算精度：广义最小二乘拟合相关矩阵负对角元问题

最近在做高度相关数据的广义最小二乘（GLS）拟合，碰到了个头疼的问题——拟合出的最佳参数相关矩阵居然出现了负对角项，这明显不符合统计逻辑嘛。先给大家梳理下我的场景和当前实现：

变量说明 #

用到的核心变量定义很清晰：

• X：预测变量矩阵
• Y：因变量向量
• C：Y的协方差矩阵（规模较大，多数元素值接近1）
• BF：求解得到的最佳拟合参数
• BFCM：BF的协方差矩阵

当前代码实现 #

我用pylab写的核心计算代码如下：

from pylab import *
BF = inv(X.T @ inv(C) @ X) @ (X.T @ inv(C) @ Y)
BFCM = inv(X.T @ inv(C) @ X)

遇到的问题 #

计算完成后发现，BF对应的相关矩阵中出现了负对角项——正常来说相关矩阵的对角元素应该恒为1，这显然是数值计算精度不足导致的误差：由于C矩阵规模大且高度相关，直接求逆的过程中误差被不断放大，最终出现了不符合预期的结果。

优化方案 #

针对这个问题，我整理了几个提升运算精度、减少数值误差的可行办法：

1. 切换更高精度的数据类型 #

如果你的运行环境支持（比如numpy支持float128类型），将所有数组转换为更高精度的类型，能从根源上降低计算过程中的截断误差：

import numpy as np
X = X.astype(np.float128)
C = C.astype(np.float128)
Y = Y.astype(np.float128)

2. 避免直接求逆，用线性求解替代 #

直接调用inv()求逆很容易放大数值误差，改用numpy.linalg.solve求解线性方程组会更稳定。比如把BF的计算逻辑改写为：

# 先计算中间矩阵M
M = X.T @ inv(C) @ X
BF = np.linalg.solve(M, X.T @ inv(C) @ Y)
BFCM = np.linalg.inv(M)

更进一步，连inv(C)都可以用solve替代，彻底规避直接求逆的风险：

# 通过求解线性方程组得到 inv(C)@X 和 inv(C)@Y
invC_X = np.linalg.solve(C, X)
invC_Y = np.linalg.solve(C, Y)
M = X.T @ invC_X
BF = np.linalg.solve(M, X.T @ invC_Y)
BFCM = np.linalg.inv(M)

3. 加入正则化项提升矩阵稳定性 #

由于数据高度相关，中间矩阵M可能接近奇异（条件数极差），给M添加一个极小的单位矩阵正则项，能有效提升矩阵的可逆性，减少求逆时的数值波动：

eps = 1e-8  # 可根据数据实际情况调整大小
M_reg = M + eps * np.eye(M.shape[0])
BF = np.linalg.solve(M_reg, X.T @ invC_Y)
BFCM = np.linalg.inv(M_reg)

这些方法都能有效缓解数值精度问题，解决相关矩阵出现负对角项的异常情况。

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