# GP regularization effect

**URL:** <https://discourse.pymc.io/t/gp-regularization-effect/6547>\
**Category:** Questions\
**Created:** [January 8, 2021, 1:17pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547 "2021-01-08T13:17:42Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![luisroque](https://avatars.discourse-cdn.com/v4/letter/l/7c8e57/32.png) [@luisroque](https://discourse.pymc.io/u/luisroque)\
**Post date:** [January 8, 2021, 1:17pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/1 "2021-01-08T13:17:42Z")

</div>

Hi,

In this case, I have 4 kernels, one for long term trend, one for short term trend, one for seasonal behavior and 1 for noise. I need to regularize the effect of the short term trend, otherwise, it overfits pretty quickly and the predictions go through the roof. The length scale parameters have prior distributions defined accordingly to the article from [Michael Betancourt](https://betanalpha.github.io/assets/case_studies/gaussian_processes.html), but he doesn’t touch much on the amplitudes.

When trying to regularize the hyperparameters of a GP I often get the following error:

> **trace**
>
> ```auto
> ---------------------------------------------------------------------------
> ValueError Traceback (most recent call last)
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/compile/function/types.py in __call__ (self, *args, **kwargs)
> 972 outputs = (
> --> 973 self.fn()
> 974 if output_subset is None
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gof/op.py in rval(p, i, o, n)
> 912 def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
> --> 913 r = p(n, [x[0] for x in i], o)
> 914 for o in node.outputs:
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/tensor/slinalg.py in perform(self, node, inputs, output_storage)
> 251 elif self.A_structure == "upper_triangular":
> --> 252 rval = scipy.linalg.solve_triangular(A, b, lower=False)
> 253 else:
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/linalg/basic.py in solve_triangular(a, b, trans, lower, unit_diagonal, overwrite_b, debug, check_finite)
> 333 a1 = _asarray_validated(a, check_finite=check_finite)
> --> 334 b1 = _asarray_validated(b, check_finite=check_finite)
> 335 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/_lib/_util.py in _asarray_validated(a, check_finite, sparse_ok, objects_ok, mask_ok, as_inexact)
> 271 toarray = np.asarray_chkfinite if check_finite else np.asarray
> --> 272 a = toarray(a)
> 273 if not objects_ok:
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
> 484 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
> --> 485 raise ValueError(
> 486 "array must not contain infs or NaNs")
> 
> ValueError: array must not contain infs or NaNs
> 
> During handling of the above exception, another exception occurred:
> 
> ValueError Traceback (most recent call last)
> <ipython-input-139-25fc0a8509e1> in <module>
> 40 mu=tt.exp(f),
> 41 observed=groups['train']['data'].reshape(groups['train']['s'], groups['train']['n']).T[:,i])
> ---> 42 mp = pm.find_MAP(maxeval=5000, progressbar = False)
> 43 
> 44 pred_samples_fit = pm.sample_posterior_predictive([mp], 
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/pymc3/tuning/starting.py in find_MAP(start, vars, method, return_raw, include_transformed, progressbar, maxeval, model, *args, **kwargs)
> 163 
> 164 try:
> --> 165 opt_result = minimize(
> 166 cost_func, x0, method=method, jac=compute_gradient, *args, **kwargs
> 167 )
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
> 615 **options)
> 616 elif meth == 'l-bfgs-b':
> --> 617 return _minimize_lbfgsb(fun, x0, args, jac, bounds,
> 618 callback=callback, **options)
> 619 elif meth == 'tnc':
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
> 358 # until the completion of the current minimization iteration.
> 359 # Overwrite f and g:
> --> 360 f, g = func_and_grad(x)
> 361 elif task_str.startswith(b'NEW_X'):
> 362 # new iteration
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in fun_and_grad(self, x)
> 198 if not np.array_equal(x, self.x):
> 199 self._update_x_impl(x)
> --> 200 self._update_fun()
> 201 self._update_grad()
> 202 return self.f, self.g
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in _update_fun(self)
> 164 def _update_fun(self):
> 165 if not self.f_updated:
> --> 166 self._update_fun_impl()
> 167 self.f_updated = True
> 168 
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in update_fun()
> 71 
> 72 def update_fun():
> ---> 73 self.f = fun_wrapped(self.x)
> 74 
> 75 self._update_fun_impl = update_fun
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in fun_wrapped(x)
> 68 def fun_wrapped(x):
> 69 self.nfev += 1
> ---> 70 return fun(x, *args)
> 71 
> 72 def update_fun():
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/optimize.py in __call__ (self, x, *args)
> 72 def __call__ (self, x, *args):
> 73 """ returns the the function value """
> ---> 74 self._compute_if_needed(x, *args)
> 75 return self._value
> 76 
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/optimize/optimize.py in _compute_if_needed(self, x, *args)
> 66 if not np.all(x == self.x) or self._value is None or self.jac is None:
> 67 self.x = np.asarray(x).copy()
> ---> 68 fg = self.fun(x, *args)
> 69 self.jac = fg[1]
> 70 self._value = fg[0]
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/pymc3/tuning/starting.py in __call__ (self, x)
> 227 value = -1.0 * nan_to_high(neg_value)
> 228 if self.use_gradient:
> --> 229 neg_grad = self.dlogp_func(pm.floatX(x))
> 230 if np.all(np.isfinite(neg_grad)):
> 231 self.previous_x = x
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/pymc3/blocking.py in __call__ (self, x)
> 270 
> 271 def __call__ (self, x):
> --> 272 return self.fa(self.fb(x))
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/pymc3/model.py in __call__ (self, state)
> 1537 
> 1538 def __call__ (self, state):
> -> 1539 return self.f(**state)
> 1540 
> 1541 
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/compile/function/types.py in __call__ (self, *args, **kwargs)
> 984 if hasattr(self.fn, "thunks"):
> 985 thunk = self.fn.thunks[self.fn.position_of_error]
> --> 986 gof.link.raise_with_op(
> 987 node=self.fn.nodes[self.fn.position_of_error],
> 988 thunk=thunk,
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gof/link.py in raise_with_op(node, thunk, exc_info, storage_map)
> 341 # Some exception need extra parameter in inputs. So forget the
> 342 # extra long error message in that case.
> --> 343 raise exc_value.with_traceback(exc_trace)
> 344 
> 345 
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/compile/function/types.py in __call__ (self, *args, **kwargs)
> 971 try:
> 972 outputs = (
> --> 973 self.fn()
> 974 if output_subset is None
> 975 else self.fn(output_subset=output_subset)
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gof/op.py in rval(p, i, o, n)
> 911 # default arguments are stored in the closure of `rval`
> 912 def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
> --> 913 r = p(n, [x[0] for x in i], o)
> 914 for o in node.outputs:
> 915 compute_map[o][0] = True
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/theano/tensor/slinalg.py in perform(self, node, inputs, output_storage)
> 250 rval = scipy.linalg.solve_triangular(A, b, lower=True)
> 251 elif self.A_structure == "upper_triangular":
> --> 252 rval = scipy.linalg.solve_triangular(A, b, lower=False)
> 253 else:
> 254 rval = scipy.linalg.solve(A, b)
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/linalg/basic.py in solve_triangular(a, b, trans, lower, unit_diagonal, overwrite_b, debug, check_finite)
> 332 
> 333 a1 = _asarray_validated(a, check_finite=check_finite)
> --> 334 b1 = _asarray_validated(b, check_finite=check_finite)
> 335 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
> 336 raise ValueError('expected square matrix')
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/scipy/_lib/_util.py in _asarray_validated(a, check_finite, sparse_ok, objects_ok, mask_ok, as_inexact)
> 270 raise ValueError('masked arrays are not supported')
> 271 toarray = np.asarray_chkfinite if check_finite else np.asarray
> --> 272 a = toarray(a)
> 273 if not objects_ok:
> 274 if a.dtype is np.dtype('O'):
> 
> ~/anaconda3/envs/phd/lib/python3.8/site-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
> 483 a = asarray(a, dtype=dtype, order=order)
> 484 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
> --> 485 raise ValueError(
> 486 "array must not contain infs or NaNs")
> 487 return a
> 
> ValueError: array must not contain infs or NaNs
> Apply node that caused the error: Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}(InplaceDimShuffle{1,0}.0, Elemwise{Composite{((i0 * i1) - i2)}}[(0, 0)].0)
> Toposort index: 72
> Inputs types: [TensorType(float64, matrix), TensorType(float64, matrix)]
> Inputs shapes: [(220, 220), (220, 220)]
> Inputs strides: [(1760, 8), (1760, 8)]
> Inputs values: ['not shown', 'not shown']
> Outputs clients: [[InplaceDimShuffle{1,0}(Solve{A_structure='upper_triangular', lower=False, overwrite_A=False, overwrite_b=False}.0)]]
> 
> Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1413, in access_grad_cache
> term = access_term_cache(node)[idx]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1064, in access_term_cache
> output_grads = [access_grad_cache(var) for var in node.outputs]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1064, in <listcomp>
> output_grads = [access_grad_cache(var) for var in node.outputs]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1413, in access_grad_cache
> term = access_term_cache(node)[idx]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1064, in access_term_cache
> output_grads = [access_grad_cache(var) for var in node.outputs]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1064, in <listcomp>
> output_grads = [access_grad_cache(var) for var in node.outputs]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1413, in access_grad_cache
> term = access_term_cache(node)[idx]
> File "/home/mach1ne/anaconda3/envs/phd/lib/python3.8/site-packages/theano/gradient.py", line 1223, in access_term_cache
> input_grads = node.op.L_op(inputs, node.outputs, new_output_grads)
> 
> HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.
> 
> ```

Here you can find some toy data:

> **data**
>
> ```auto
> y = np.asarray([936., 140., 118., 202., 122., 96., 124., 85., 8., 79., 82.,
> 282., 455., 219., 103., 124., 132., 191., 99., 78., 108., 113.,
> 140., 145., 574., 54., 52., 343., 66., 169., 174., 88., 121.,
> 140., 77., 251., 445., 220., 108., 343., 78., 54., 78., 80.,
> 68., 268., 185., 128., 248., 286., 85., 134., 84., 76., 67.,
> 72., 57., 158., 199., 238., 592., 175., 216., 298., 74., 68.,
> 148., 49., 82., 150., 125., 161., 579., 231., 79., 308., 164.,
> 68., 130., 104., 46., 93., 105., 116., 539., 100., 114., 129.,
> 50., 94., 56., 128., 93., 168., 92., 128., 792., 195., 89.,
> 246., 53., 147., 157., 86., 82., 213., 84., 145., 827., 208.,
> 114., 175., 62., 65., 118., 74., 92., 203., 138., 163., 715.,
> 166., 166., 133., 88., 65., 57., 123., 132., 121., 149., 171.,
> 726., 220., 71., 173., 24., 67., 127., 90., 82., 159., 113.,
> 187., 510., 131., 101., 189., 80., 82., 115., 99., 71., 120.,
> 121., 122., 615., 106., 88., 156., 91., 76., 140., 64., 84.,
> 147., 88., 165., 639., 204., 93., 103., 96., 82., 111., 119.,
> 109., 220., 94., 138., 514., 141., 91., 303., 103., 163., 116.,
> 77., 122., 87., 104., 205., 658., 94., 128., 126., 39., 75.,
> 71., 34., 66., 175., 89., 155., 831., 95., 210., 130., 31.,
> 95., 59., 95., 111., 136., 146., 184., 510., 101., 210., 138.])
> 
> ```

Model:

```auto
X = np.arange(220).reshape(-1,1)
with pm.Model() as model:

    l_t = pm.InverseGamma('l_t', 4, 220)
    l_ts = pm.Gamma('l_ts', 5, 1)
    l_p = pm.InverseGamma('l_p', 4, 220)
    η_trend = pm.HalfNormal('η_trend', 1)
    η_per = pm.HalfNormal('η_per', 1.5)
    η_ts = pm.HalfNormal('η_ts',0.15)
    σ = pm.HalfNormal("σ", sigma=0.05)
    
    a = pm.Normal('a', 
                  0.0,
                  5)

    # cov function for the GP 
    cov = (η_trend**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=l_t)
            + η_ts**2 * pm.gp.cov.ExpQuad(input_dim=1, ls=l_ts)
            + η_per**2 * pm.gp.cov.Periodic(1, period=12, ls=l_p)
            + pm.gp.cov.WhiteNoise(σ))

    gp = pm.gp.Latent(mean_func=pm.gp.mean.Constant(a), cov_func=cov)
    f = gp.prior('f', X=X, reparameterize=True)

    y_pred = pm.Poisson('y_pred', mu=tt.exp(f), observed=y)
    #prior_checks = pm.sample_prior_predictive(200)
    mp = pm.find_MAP(maxeval=5000, progressbar = True)

```

I understand that it is due to the fact that I am using a prior distribution with much mass near zero. Is there any other way to achieve this regularizing effect? What exactly is happening that triggers the error, is there any way to debug it a little deeper?

Tks,  
Luis

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [January 8, 2021, 7:50pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/2 "2021-01-08T19:50:09Z")

</div>

Yes it does look like you’re getting those errors because one or more of the etas is going to zero. Maybe you could try using Gammas or Lognormals for the eta parameters? Or other boundary avoiding priors?

---

<div class="post-metadata">

**Author:** ![Martin\_Ingram](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/martin_ingram/32/3535_2.png) [@Martin\_Ingram](https://discourse.pymc.io/u/Martin_Ingram)\
**Post date:** [January 9, 2021, 1:51am UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/3 "2021-01-09T01:51:10Z")

</div>

It looks like you’re trying to do MAP estimation here? Even with regularisation, that typically doesn’t work, unless you marginalise out the GP, because the posterior mode may well be at these pathological points. A similar thing can happen in hierarchical models where setting the prior variance of the random effects to zero gives a spike in the posterior distribution (there’s a bit of a discussion of that here, under “Why hierarchical models are Bayesian”, for example: [Why hierarchical models are awesome, tricky, and Bayesian — While My MCMC Gently Samples](https://twiecki.io/blog/2017/02/08/bayesian-hierchical-non-centered/)). I’ll try to dig up more references about this, but in short I’d suggest using the NUTS sampler or maybe trying variational inference.

---

<div class="post-metadata">

**Author:** ![luisroque](https://avatars.discourse-cdn.com/v4/letter/l/7c8e57/32.png) [@luisroque](https://discourse.pymc.io/u/luisroque)\
**Post date:** [January 9, 2021, 2:42pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/4 "2021-01-09T14:42:57Z")

</div>

Hi @bwengals,

My problem with changing to those priors is that the model quickly overfits the data (especially in the short term trend case). This results in estimations that go through the roof when the conditioning on new points (forecasting). I need some regularization effect here.

---

<div class="post-metadata">

**Author:** ![luisroque](https://avatars.discourse-cdn.com/v4/letter/l/7c8e57/32.png) [@luisroque](https://discourse.pymc.io/u/luisroque)\
**Post date:** [January 9, 2021, 3:04pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/5 "2021-01-09T15:04:28Z")

</div>

Hi @Martin_Ingram,

I completely agree and would do that if it was computationally doable. The problem is that I am using GPs as priors in a hierarchical model and running it with datasets with 500+ time series with 300+ points each.

Running NUTs with a simple model like this and just one time series takes several hours because of the GP complexity.

Any suggestion on how to improve performance and ensure that the posterior space is free of these pathological shapes?

---

<div class="post-metadata">

**Author:** ![Martin\_Ingram](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/martin_ingram/32/3535_2.png) [@Martin\_Ingram](https://discourse.pymc.io/u/Martin_Ingram)\
**Post date:** [January 9, 2021, 10:28pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/6 "2021-01-09T22:28:03Z")

</div>

I see, that does make things harder! I’d say there’s no perfect answer here, but some suggestions below.

The first thing I would probably check is whether it might be reasonable to make a Gaussian approximation to your Poisson likelihood. If your counts are large, this might be reasonable (see e.g. [Poisson normal approximation error](http://www.johndcook.com/blog/normal_approx_to_poisson/)). Then you could marginalise out the GP, at which point MAP estimation for the hyperparameters usually works fine. By the way, I just noticed your toy data is apparently continuous – it should be discrete for a Poisson! That might also be causing some numerical issues.

Another option could be to see if you could estimate the hyperparameters with NUTS on some smaller subset of data, and then turn them into a point estimate, e.g. by just using their mean, when fitting / predicting on your full set of time series. It might also be worth checking out the JAX-based samplers (see [here](https://pymc-devs.medium.com/the-future-of-pymc3-or-theano-is-dead-long-live-theano-d8005f8a0e9b)) and trying NUTS with those on a GPU – there can definitely be big speedups.

You could also try a sparse approximation to the GPs: [Sparse Approximations — PyMC3 3.10.0 documentation](https://docs.pymc.io/notebooks/GP-SparseApprox.html) . The idea there would roughly be that you approximate the GPs with 300+ points with ones with say 100, which can speed things up.

So no great answer from me here, sorry. Maybe someone else has a clear answer, but scaling GPs is often a bit tricky. You could do things with variational inference too, there are special techniques developed for GPs (e.g. [here](https://arxiv.org/abs/1411.2005)) that can work well, but you would likely have to do this outside pymc3.

---

<div class="post-metadata">

**Author:** ![luisroque](https://avatars.discourse-cdn.com/v4/letter/l/7c8e57/32.png) [@luisroque](https://discourse.pymc.io/u/luisroque)\
**Post date:** [January 11, 2021, 1:18pm UTC](https://discourse.pymc.io/t/gp-regularization-effect/6547/7 "2021-01-11T13:18:39Z")

</div>

Thank you for the detailed answer Martin!

> [@Martin\_Ingram](#):
>
> The first thing I would probably check is whether it might be reasonable to make a Gaussian approximation to your Poisson likelihood. If your counts are large, this might be reasonable (see e.g. [Poisson normal approximation error](http://www.johndcook.com/blog/normal_approx_to_poisson/)). Then you could marginalise out the GP, at which point MAP estimation for the hyperparameters usually works fine.

I started with this approach, but I had to change it pretty quickly as I was having negative estimations due to the flexibility of the GPs when fitting the data. I can’t really have that.

> [@Martin\_Ingram](#):
>
> By the way, I just noticed your toy data is apparently continuous – it should be discrete for a Poisson! That might also be causing some numerical issues.

This was a copy-paste error, I am indeed always using counts. Thank you for noticing.

> [@Martin\_Ingram](#):
>
> Another option could be to see if you could estimate the hyperparameters with NUTS on some smaller subset of data, and then turn them into a point estimate, e.g. by just using their mean, when fitting/predicting on your full set of time series.

Do you think that this is a better alternative? I would say that it isn’t because I am actually fitting my model to some subset of the data to recover hyperparameters, fitting it again to a different dataset (that includes the subset but that it is nonetheless significantly different) to recover parameters, and then conditioning on new data. It seems a bit tricky.

> [@Martin\_Ingram](#):
>
> You could also try a sparse approximation to the GPs: [Sparse Approximations — PyMC3 3.10.0 documentation](https://docs.pymc.io/notebooks/GP-SparseApprox.html) . The idea there would roughly be that you approximate the GPs with 300+ points with ones with say 100, which can speed things up.

I also did some digging here, but there is no current solution to do sparse approximations with latent GPs.

> [@Martin\_Ingram](#):
>
> You could do things with variational inference too, there are special techniques developed for GPs (e.g. [here](https://arxiv.org/abs/1411.2005)) that can work well, but you would likely have to do this outside pymc3.

I am gonna check this reference.

Thank you!