# Gaussian process regression

**URL:** <https://discourse.mc-stan.org/t/gaussian-process-regression/21526>\
**Category:** General\
**Created:** [March 26, 2021, 10:09pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526 "2021-03-26T22:09:28Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![linas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/linas/32/5342_2.png) [@linas](https://discourse.mc-stan.org/u/linas)\
**Post date:** [March 26, 2021, 10:09pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/1 "2021-03-26T22:09:28Z")

</div>

Hello,

I was trying to run Gaussian process regression in Stan with 5000 observations and 14 factors. It runs very very slowly. Are there any packages that run such problems? The Stan code is attached. The code does well for 100 observations

Thanks for any advice.  
[GP.stan](https://discourse.mc-stan.org/uploads/short-url/3qeAI0M3QvdZhOM1kHlk1LSym6k.stan) (3.2 KB)

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [March 29, 2021, 5:55pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/2 "2021-03-29T17:55:48Z")

</div>

I don’t see any obvious routes to speedups. Presumably you’ve already tried on a GPU?

One thing that’s unlikely to help much but worth a try if you’re curious is to precompute the set of unique differences in `X`. This saves a small amount of compute during sampling but if the set of unique differences is actually much smaller than the total number of differences, you might save more substantial compute. In my tests long ago I found the speedup from this didn’t match simply using cov\_exp\_quad, but possibly since you’re not using cov\_exp\_quad yourself it might come in handy.

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [March 29, 2021, 5:58pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/3 "2021-03-29T17:58:21Z")

</div>

Link to code from ages ago: [Comparing cov\_exp\_quad to alternative gp optimizations · GitHub](https://gist.github.com/mike-lawrence/b56ac2d9450233837d426dd17fe20557)

---

<div class="post-metadata">

**Author:** ![linas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/linas/32/5342_2.png) [@linas](https://discourse.mc-stan.org/u/linas)\
**Post date:** [March 30, 2021, 10:43pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/4 "2021-03-30T22:43:03Z")

</div>

Thanks a lot. In my case I have 14 factors. This example is one factor.  
Also I wonder if there are any R packages which do Gaussian process regression? I am specifically interested in the package by Aki Vehtari.

---

<div class="post-metadata">

**Author:** ![jtimonen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jtimonen/32/9034_2.png) [@jtimonen](https://discourse.mc-stan.org/u/jtimonen)\
**Post date:** [March 30, 2021, 11:03pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/5 "2021-03-30T23:03:34Z")

</div>

you can do additive GP regression with [Longitudinal Gaussian Process Regression • lgpr](https://jtimonen.github.io/lgpr-usage/) but it won’t be any faster than your code

---

<div class="post-metadata">

**Author:** ![jtimonen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jtimonen/32/9034_2.png) [@jtimonen](https://discourse.mc-stan.org/u/jtimonen)\
**Post date:** [March 30, 2021, 11:06pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/6 "2021-03-30T23:06:02Z")

</div>

also brms: [gp: Set up Gaussian process terms in 'brms' in brms: Bayesian Regression Models using 'Stan'](https://rdrr.io/cran/brms/man/gp.html)

---

<div class="post-metadata">

**Author:** ![jtimonen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jtimonen/32/9034_2.png) [@jtimonen](https://discourse.mc-stan.org/u/jtimonen)\
**Post date:** [March 31, 2021, 2:09am UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/7 "2021-03-31T02:09:35Z")

</div>

Those use Stan to sample kernel parameters. If you want to just optimize hyperparameters, there is for example gplite: [gplite Quickstart](https://cran.r-project.org/web/packages/gplite/vignettes/quickstart.html)

---

<div class="post-metadata">

**Author:** ![jbaranowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jbaranowski/32/5734_2.png) [@jbaranowski](https://discourse.mc-stan.org/u/jbaranowski)\
**Post date:** [March 31, 2021, 5:59am UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/8 "2021-03-31T05:59:40Z")

</div>

I was not reading it very carefully, but I’ve noticed that your calcP function requires you to do a nested for loop, by the number of observations on each level. It is certainly a delaying factor on as you need O(25mln) multiplications per every step. Maybe you can vectorize it?

---

<div class="post-metadata">

**Author:** ![jtimonen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jtimonen/32/9034_2.png) [@jtimonen](https://discourse.mc-stan.org/u/jtimonen)\
**Post date:** [March 31, 2021, 8:37am UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/9 "2021-03-31T08:37:28Z")

</div>

The Cholesky decomposition is still the computationally most demanding part there. I am certainly interested if there is some way to get less autodiff variables and speed up computation that way, but I am surprised if there is any way to get that implementation running faster than several days for \> 2000 observations.

---

<div class="post-metadata">

**Author:** ![linas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/linas/32/5342_2.png) [@linas](https://discourse.mc-stan.org/u/linas)\
**Post date:** [March 31, 2021, 8:06pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/10 "2021-03-31T20:06:45Z")

</div>

It would be nice. Do you have any suggestions?

---

<div class="post-metadata">

**Author:** ![jbaranowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jbaranowski/32/5734_2.png) [@jbaranowski](https://discourse.mc-stan.org/u/jbaranowski)\
**Post date:** [April 1, 2021, 6:38am UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/11 "2021-04-01T06:38:31Z")

</div>

Maybe use already implemented covariance function? [5.13 Covariance functions | Stan Functions Reference](https://mc-stan.org/docs/2_26/functions-reference/covariance.html)  
It will certainly be more efficient than an explicit loop

---

<div class="post-metadata">

**Author:** ![linas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/linas/32/5342_2.png) [@linas](https://discourse.mc-stan.org/u/linas)\
**Post date:** [April 1, 2021, 4:46pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/12 "2021-04-01T16:46:08Z")

</div>

Thank you. The problem is that I need length scale parameter to be a vector. The current implementation is real.

---

<div class="post-metadata">

**Author:** ![jbaranowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jbaranowski/32/5734_2.png) [@jbaranowski](https://discourse.mc-stan.org/u/jbaranowski)\
**Post date:** [April 2, 2021, 2:06pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/13 "2021-04-02T14:06:32Z")

</div>

Oh, I have not noticed that. But vectorization will be easy.  
You need to generate two matrices of repeated X’s one as rows and one as columns. Subtract them from each other, multiply by a diagonal matrix of rho’s from the appropriate side, call element wise square, divide by two and exponentiate.  
Unless I’ve missed something it should work.

---

<div class="post-metadata">

**Author:** ![jtimonen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jtimonen/32/9034_2.png) [@jtimonen](https://discourse.mc-stan.org/u/jtimonen)\
**Post date:** [April 2, 2021, 5:53pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/14 "2021-04-02T17:53:48Z")

</div>

That is how you could do it if you had a one-dimensional x and you want the cov\_exp\_quad kernel, but here the ARD kernel was used and each x is a vector with length Ncol (=14?).

---

<div class="post-metadata">

**Author:** ![jtimonen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jtimonen/32/9034_2.png) [@jtimonen](https://discourse.mc-stan.org/u/jtimonen)\
**Post date:** [April 2, 2021, 5:57pm UTC](https://discourse.mc-stan.org/t/gaussian-process-regression/21526/15 "2021-04-02T17:57:48Z")

</div>

But I just have the feeling that even if you can avoid a loop in kernel matrix computation, it’s not going to help a lot here, because you still have to do cholesky for the 5000 x 5000 matrix and have to sample the 5000 eta parameters