[Python] numpy matmul raises divide by zero warning

I'm implementing Gaussian Process Regression and am currently working on obtaining samples from the GP. For this, I'm using Equation 2.19 from Rasmussen/Williams:



However, I'm running into a weird problem when calculating the posterior mean where a matmul operation produces a divide by zero warning. I'm trying to understand why a matmul operation would even do a division in the first place. It also seems that the error only occurs when predicting for 500ish+ points, which leads me to believe there might be a memory connection. Can anyone explain why this error occurs and how I can stop it from happening?

Code:


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


Environment info/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",
]

Continue reading...