GP-MeansAndCovs.ipynb throws LinAlgError

[MarcoGorelli]

Cell 6 contains:

lengthscale = 0.2
eta = 2.0
cov = eta ** 2 * pm.gp.cov.ExpQuad(1, lengthscale)

X = np.linspace(0, 2, 200)[:, None]
K = cov(X).eval()

plt.figure(figsize=(14, 4))

plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K).random(size=3).T)
plt.title("Samples from the GP prior")
plt.ylabel("y")
plt.xlabel("X");

which throws

ValueError: input operand has more dimensions than allowed by the axis remapping

As per #11 , I can fix this by specifying the shape of the MvNormal:

lengthscale = 0.2
eta = 2.0
cov = eta ** 2 * pm.gp.cov.ExpQuad(1, lengthscale)

X = np.linspace(0, 2, 200)[:, None]
K = cov(X).eval()

plt.figure(figsize=(14, 4))

plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K, shape=K.shape[0]).random(size=3).T)
plt.title("Samples from the GP prior")
plt.ylabel("y")
plt.xlabel("X");

but then I get

~/pymc3-dev/pymc3/distributions/multivariate.py in random(self, point, size)
    277 
    278         if self._cov_type == "cov":
--> 279             chol = np.linalg.cholesky(param)
    280         elif self._cov_type == "chol":
    281             chol = param

<__array_function__ internals> in cholesky(*args, **kwargs)

~/miniconda3/envs/pymc3-dev-py38/lib/python3.8/site-packages/numpy/linalg/linalg.py in cholesky(a)
    762     t, result_t = _commonType(a)
    763     signature = 'D->D' if isComplexType(t) else 'd->d'
--> 764     r = gufunc(a, signature=signature, extobj=extobj)
    765     return wrap(r.astype(result_t, copy=False))
    766 

~/miniconda3/envs/pymc3-dev-py38/lib/python3.8/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_nonposdef(err, flag)
     89 
     90 def _raise_linalgerror_nonposdef(err, flag):
---> 91     raise LinAlgError("Matrix is not positive definite")
     92 
     93 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):

LinAlgError: Matrix is not positive definite

cc @Sayam753 any suggestions?

[Sayam753]

So, the covariance matrix is singular.

That's why Cholesky decomposition fails. I would add a small WhiteNoise to stabilize the diagonal.

Before MvNormal.random was fixed for shape errors in pymc-devs/pymc#4207, it supported cov being singular using scipy.stats. (see before 3.10.0)

I can lookup the implementation in scipy this week and add to PyMC3 to also support singular matrices.

[MarcoGorelli]

cool, thanks for your explanation!

Do you to change it to, e.g.

K = cov(X).eval()
K += np.random.randn(*K.shape)*1e-1

?

If so, np.linalg.det(K) becomes non-zero, but I still get the same linalg error

[Sayam753]

Something like

cov = eta ** 2 * pm.gp.cov.ExpQuad(1, lengthscale) + pm.gp.cov.WhiteNoise(1e-6)

[MarcoGorelli]

Ah, makes sense, thanks, that works!