LinAlgError: Matrix is not positive definite error when generating fake data in HSGP-Basic

[Armatron44]

Notebook title: HSGP-Basic
Notebook url: https://github.com/pymc-devs/pymc-examples/blob/main/examples/gaussian_processes/HSGP-Basic.ipynb

Issue description

When executing the code in fourth notebook cell (trying to plot some simulated data) a LinAlgError exception is raised with message "Matrix is not positive definite". This exception is thrown when executing f_true = pm.draw(gp_true, draws=1, random_seed=rng) in the simulated_1d function defined in the code of the third cell.

Tested on windows 10 with
pymc 5.25.1
pytensor 2.31.7
numpy 2.2.5
python 3.12.11
and also recreated on a google colab environment.

Expected output

Output should be matplotlib figure, as saved in output of hosted notebook.

Proposed solution

As this appears to be a numerical issue from decomposition for matrix inversion, I found passing the covariance function to pm.gp.util.stabilize stopped the exception being raised. Alternatively, one could reduce the size of x but that feels less robust given its a data-dependent solution.

I suggest changing the code in the third cell from

def simulate_1d(x, ell_true, eta_true, sigma_true):
...
cov_func = eta_true**2 * pm.gp.cov.Matern52(1, ell_true)
gp_true = pm.MvNormal.dist(mu=np.zeros(n), cov=cov_func(x[:, None]))
...

to

def simulate_1d(x, ell_true, eta_true, sigma_true):
...
cov_func = eta_true**2 * pm.gp.cov.Matern52(1, ell_true)
cov_mat_jitted = pm.gp.util.stabilize(cov_func(x[:, None]))
gp_true = pm.MvNormal.dist(mu=np.zeros(n), cov=cov_mat_jitted)
...