numpy matmul触发除零警告：高斯过程回归采样问题求助

我正在实现高斯过程回归（Gaussian Process Regression），需要从高斯过程（GP）中获取样本，采用Rasmussen/Williams著作中的公式2.19。但在计算后验均值时，遇到异常问题：matmul操作触发了divide by zero警告。我无法理解矩阵乘法操作为何会涉及除法，且该问题仅在预测500个以上点时出现，怀疑与内存有关。恳请解释该错误的成因并提供解决方法。

代码 #

import numpy as np
from numpy.random import RandomState

np.seterr(all='raise') # to error out

SEED = RandomState(42)
JITTER = 1e-6
observed_x = np.array([0.0, 1.0, 3.0])
observed_y = np.array([1., -1.0, 3.0])

noise = SEED.randn(10, 1000)

lengthscale = 1.
k_ = lambda xi, xj: np.exp(- 1 * (xi[:, None] - xj) ** 2 / (2 * lengthscale ** 2))
K = (
    k_(observed_x, observed_x)
)
K += JITTER * np.eye(len(observed_x))

x_star = np.linspace(-1, 4, 1000)
K_star = k_(x_star, observed_x)
L = np.linalg.cholesky(K)
alpha = np.linalg.solve(L.T, np.linalg.solve(L, observed_y))

mu_pred = K_star @ alpha

环境信息（pyproject.toml） #

[project]
name = "gp"
version = "0.1.0"
readme = "README.md"
requires-python = ">=3.10"
dependencies = [
    "ipykernel>=6.29.5",
    "ipympl>=0.9.7",
    "jupyter>=1.1.1",
    "matplotlib>=3.10.1",
    "numpy>=2.2.4",
    "scikit-learn>=1.6.1",
    "seaborn>=0.13.2",
    "notebook==6.5.6",
]

原因分析 #

1. 警告并非来自matmul本身：矩阵乘法本质是加法和乘法运算，不会直接触发除法错误。实际问题大概率出现在之前的Cholesky分解、线性求解步骤，但由于numpy的延迟计算机制，警告被延迟到matmul执行时才抛出。
2. 数值稳定性不足：尽管添加了JITTER保证协方差矩阵K正定，但1e-6的JITTER可能过小，当预测点数量增多时，数值误差被放大，导致Cholesky分解后出现接近零的对角元素，后续求解alpha时触发隐性除法（线性求解底层涉及除法操作）。
3. numpy版本特性：你使用的numpy 2.x版本在数值精度控制或警告机制上与1.x版本存在差异，可能将原本被忽略的数值问题以警告形式暴露。

解决方法 #

• 增大JITTER值：将JITTER从1e-6调整为1e-5或1e-4，增强协方差矩阵的正定性，避免Cholesky分解出现数值退化。
• 定位错误根源：在关键步骤后添加检查，确认中间结果是否异常：

  print("L对角元素:", np.diag(L))
  print("alpha取值:", alpha)
  

• 降级numpy版本：暂时切换到numpy 1.26.x版本，验证是否是2.x版本的兼容性问题。
• 替换线性求解方式：用矩阵逆运算替代Cholesky分解求解（仅用于调试，逆运算效率较低）：

  alpha = np.linalg.inv(K) @ observed_y
  

• 优化核函数实现：使用scikit-learn提供的RBF核函数，其数值稳定性优于手动实现的lambda函数。

内容的提问来源于stack exchange，提问作者till-m