I’m pretty sure I’ve seen expressions involving White Noise in the literature where some additive white noise is included in the kernel. I think this should be possible:

with pymc.Model() as model:
...
k_se = pymc.gp.cov.ExpQuad(10, ls=[1.0]*10)
k_wn = pymc.gp.cov.WhiteNoise(sigma=1)
k = k_se + k_wn
gp = pymc.gp.Latent(cov_func=k)
f = gp.prior("f", X)
...

This raises:

line 189, in <setcomp>
input_dims = {factor.input_dim for factor in factor_list if isinstance(factor, Covariance)}
AttributeError: 'WhiteNoise' object has no attribute 'input_dim'

Not sure why this is happening, these evaluate to the same shape, so elementwise addition should be possible, i.e.

X = np.random.rand(442,10)
k_se.full(X).shape.eval() # 425 x 425
k_wn.full(X).shape.eval() # 425 x 425

Often people add a small amount of “white noise” in order to make the covariance matrix numerically positive semidefinite for the Cholesky decomposition. Most of the GP classes have an argument called jitter that does this. It defaults to 1e-6, but you can increase or decrease that. I think it’s helpful to differentiate whether you’re actually adding noise, as in your stochastic process has a white noise component, or you’re just adding a bit to make the Cholesky decomposition work – even though mathematically the same thing happens.

The cases I’ve seen use the white noise kernel to model stochastic noise in the underlying physical process - I wasn’t even aware of the numerical stability case. Good to know though

pm.gp.cov.WhiteNoise won’t work with HSGP. Adding a white noise kernel with an unknown sigma parameter is equivalent to adding a Gaussian random effect to your model.

For HSGP you don’t need to add a bit to the diagonal like you do for GPs that use the Cholesky because it doesn’t calculate the covariance matrix.