# Could you please help to find out why my GP model got slower

**URL:** <https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579>\
**Category:** Modeling\
**Tags:** gaussian-process\
**Created:** [March 29, 2021, 9:21pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579 "2021-03-29T21:21:18Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [March 29, 2021, 9:21pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/1 "2021-03-29T21:21:18Z")

</div>

Dear community,

Sorry for writing so long. Maybe this is not a GP question but rather due to my poor Stan coding ability. I am playing with GP following this great case study from @betanalpha [https://betanalpha.github.io/assets/case\_studies/gaussian\_processes.html;](https://betanalpha.github.io/assets/case_studies/gaussian_processes.html;) stan manual [10.3 Fitting a Gaussian process | Stan User’s Guide](https://mc-stan.org/docs/2_26/stan-users-guide/fit-gp-section.html) and @avehtari’s HSGP paper [[2004.11408] Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming](https://arxiv.org/abs/2004.11408).

My practice is to estimate a 1D function f(x\_{i}), i = 1,2..N\_{f} which isconvolved (i.e., \*) with a vector event\_{j}, j = 1,2,...N\_{obs}, so that the observed data y\_{obs\_{j}} = event \*f + \sigma. Size N\_{f} for function f is just 100, so not computationly expensive to model it using a exponentiated quadratic covariance function. Vector event\_{j} has the same size as y\_{obs\_{j}}, N\_{obs} is around 2500. Values in event\_{j} is just binarized, so either 0 or 1. The apperaiance of 1 follows a known pattern (event[j==k]= 1; otherwise, event[j!=k] = 0, k is predefined). The idea is I simulated y\_{obs} using y\_{obs} = event \*f + \sigma, with known f, event and \sigma. Then I’d like to model f by GP and estimate it from y\_{obs} and event. Since I didn’t find a convolution fucntion that can be used inside Stan, I just toeplitzd vector event with an all 0 value vector that has the same size of f. This was done outside Stan code using Toeplitz function from the “pracma” package as following R code, the resultant matrix evt is just considered as data input to Stan model.

```
evt = t(Toeplitz(event, c(event[1], matrix(0, nrow = 1, ncol = length(f)-1))))

```

Then I can just do matrix multiplication of evt and f instead of convolution in the transformed parameter block. The model below (model #1) worked with good diagnostics and estimated well the shape of f. It took around 100s to finish the computation on my laptop. I can also mimic HSGP paper’s code (just for practicing to understand the method) to approximate the covariance function of f with m around 30 to speed up, although speed is not a big issue considering the size of f is just 100.

* * *

**model #1**

```stan
data{
    int<lower=1> N_obs;
    int<lower=1> N_f;
    int<lower=1> N_noise;
    vector[N_obs] y_obs;
    matrix[N_f,N_obs] evt;
    real x_f[N_f];
}

transformed data {
  real delta = 1e-10;
}

parameters {
  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> sigma_obs;
  vector[N_f] eta;
}

transformed parameters {
  vector[N_f] f;
  vector[N_obs] pred_obs; 
  {
    matrix[N_f, N_f] L_K;
    matrix[N_f, N_f] K = cov_exp_quad(x_f, alpha, rho);
    for (n in 1:N_f)
      K[n, n] = K[n, n] + delta;
    L_K = cholesky_decompose(K);
    f = L_K * eta ;
    pred_obs = evt' * f;
  }
}

model {
  //inv_gamma(9, 173); [10,50] , inv_gamma(4.6, 55.2); [5,50],  
  // inv_gamma(4.6, 22); [2,20], inv_gamma(7.3, 75);; [5,30], inv_gamma(4.6, 33.1); [3,30],
  rho ~ inv_gamma(4.6, 33.1);
  alpha ~ normal(0, 2);
  sigma_obs ~ std_normal();
  eta ~ std_normal();
  y_obs ~ normal(pred_obs, sigma_obs); 
}

```

**However** , for some reason, I tried to take this event as the input data and create this Toeplitz matrix evt inside the Stan code as shown in model #2, but the computation time got almost 10 times longer with the same results estimated by model #1. Toeplitz function was copied from the same function in the “pracma” package. The input and output should be exactly the same as I did outside Stan code `evt = t(Toeplitz(event, c(event[1], matrix(0, nrow = 1, ncol = length(f)-1))))`. It also doesn’t involve any parameter in Toeplitz function, but just predefined values.

**My question is** do you think the reason for slowing down is the multiplication ` pred_obs = generate_evt(event, N_f, N_obs) * f` called Toeplitz function so many times? I meant every draw for `f` will call this function once, although the function itself doesn’t involve any parameter to estimate. Is there any trick to only calculate evt once? I guess I should not put it in the transformed parameters block. However, the spectral density functions of different covariance functions are also called in transformed parameters in the birthday data [basis\_functions\_approach\_to\_GP/CaseII\_Birthday\_Approximate-GP.stan at master · gabriuma/basis\_functions\_approach\_to\_GP · GitHub](https://github.com/gabriuma/basis_functions_approach_to_GP/blob/master/Paper/Case-study_Birthday-data/CaseII_Birthday_Approximate-GP.stan).

The reason why I’d like to do this Toeplitz function inside Stan is that in the next step I would simulate data with a new event vector in which event[k] = alpha[k] instead of event[k] = 1, then estimate both f and alpha[k] as parameters.

* * *

**model #2**

```stan
functions {
  matrix Toeplitz(vector a, vector b){
  int n = rows(a);
  int m = rows(b);
  matrix[n, m] T;
  T[1,] = b';
  T[, 1] = a;
  for (i in 2:n) {
      T[i, 2:m] = T[i - 1, 1:(m - 1)];
  }
  return(T);
  }
  
  matrix generate_evt(vector event, int N1, int N2, int rec) {
  matrix[N2, N1] evt;
  evt = Toeplitz(event, append_row(event[1], rep_vector(0, N1-1)));
  return evt;
  }
}

data{
    int<lower=1> N_obs;
    int<lower=1> N_f;
    vector[N_obs] y_obs;
    vector[N_obs] event;
    real x_f[N_f];
}

transformed data {
  real delta = 1e-10;
}

parameters {
  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> sigma_obs;
  vector[N_f] eta;
  vector[N_noise] cof_noise;
}

transformed parameters {
  vector[N_f] f;
  vector[N_obs] pred_obs; 
  {
    matrix[N_f, N_f] L_K;
    matrix[N_f, N_f] K = cov_exp_quad(x_f, alpha, rho);

    // diagonal elements
    for (n in 1:N_f)
      K[n, n] = K[n, n] + delta;
    L_K = cholesky_decompose(K);
    f= L_K * eta ;
    pred_obs = generate_evt(event, N_f, N_obs) * f;
  }
}

model {
  //inv_gamma(9, 173); [10,50] , inv_gamma(4.6, 55.2); [5,50],  
  // inv_gamma(4.6, 22); [2,20], inv_gamma(7.3, 75);; [5,30], inv_gamma(4.6, 33.1);; [3,30],
  rho ~ inv_gamma(4.6, 33.1);
  alpha ~ normal(0, 2);
  sigma_obs ~ std_normal();
  eta ~ std_normal();
  cof_noise ~ normal(0, 10);
  y_obs ~ normal(pred_obs, sigma_obs); 
}

```

Thanks a lot,  
zcai

---

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [March 30, 2021, 1:50pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/2 "2021-03-30T13:50:53Z")

</div>

Indeed it was due to my poor Stan programming ability. I found the instruction here [8.6 Program block: transformed parameters | Stan Reference Manual](https://mc-stan.org/docs/2_26/reference-manual/program-block-transformed-parameters.html) saying

> Any variable that is defined wholly in terms of data or transformed data should be declared and defined in the transformed data block. Defining such quantities in the transformed parameters block is legal, but much less efficient than defining them as transformed data.

So I moved `generate_evt(event, N_f, N_obs)` from transformed parameters to transformed data and the speed is the same as if I would have given this consistent matrix as the input data.

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [March 30, 2021, 6:40pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/3 "2021-03-30T18:40:00Z")

</div>

Glad you got it sorted! Side note, we (Stan) should figure out how to do

```stan
diagonal(A) += delta;

```

to replace

```stan
    for (n in 1:N_f)
      K[n, n] = K[n, n] + delta;

```

Though we do allow

```stan
A += diag_matrix(rep_vector(N, 1.0));

```

though I think it will be a little slower than the first piece of code.

---

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [March 30, 2021, 6:50pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/4 "2021-03-30T18:50:17Z")

</div>

Thanks, Steve. [Robust Gaussian Process Modeling](https://betanalpha.github.io/assets/case_studies/gaussian_processes.html#3_Inferring_A_Gaussian_Process) was adding this k(x\_{i}, x\_{i}) as follows,

```
  matrix[N_obs, N_obs] cov = cov_exp_quad(x_obs, alpha, rho)
                             + diag_matrix(rep_vector(square(sigma), N_obs));

```

While I copied from Stan User’s Guide [10.3 Fitting a Gaussian process | Stan User’s Guide](https://mc-stan.org/docs/2_26/stan-users-guide/fit-gp-section.html) using

```
    // diagonal elements
    for (n1 in 1:N1)
      K[n1, n1] = K[n1, n1] + delta;

```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 30, 2021, 7:36pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/5 "2021-03-30T19:36:23Z")

</div>

> [@stevebronder](#):
>
> ```no-highlight
> for (n in 1:N_f)
> K[n, n] = K[n, n] + delta;
> 
> ```

I would use [add\_diag](https://mc-stan.org/docs/2_26/functions-reference/diagonal-matrix-functions.html)

```no-highlight
K = add_diag(K, delta);

```

> [@zcai](#):
>
> While I copied from Stan User’s Guide [10.3 Fitting a Gaussian process | Stan User’s Guide](https://mc-stan.org/docs/2_26/stan-users-guide/fit-gp-section.html) using

we should update that to use `add_diag`

---

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [March 30, 2021, 8:08pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/6 "2021-03-30T20:08:16Z")

</div>

Thanks a lot, @stevebronder and @avehtari. Just for your reference, I tested both options with the same data and seed using one chain.

Using `K = add_diag(K, delta);` I got

> Chain 1: Gradient evaluation took 0 seconds  
> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.  
> …  
> Chain 1: Elapsed Time: 258.85 seconds (Warm-up)  
> Chain 1: 297.712 seconds (Sampling)  
> Chain 1: 556.562 seconds (Total)

Using `for (n in 1:N_f) K[n, n] = K[n, n] + delta;;` I got

> Chain 1: Gradient evaluation took 0.002 seconds  
> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 20 seconds.  
> …  
> Chain 1: Elapsed Time: 285.46 seconds (Warm-up)  
> Chain 1: 301.017 seconds (Sampling)  
> Chain 1: 586.477 seconds (Total)

---

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [March 31, 2021, 12:53am UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/7 "2021-03-31T00:53:13Z")

</div>

Hi @avehtari, may I ask approximately the computation time of the birthday data? [basis\_functions\_approach\_to\_GP/CaseII\_Birthday\_Approximate-GP.stan at master · gabriuma/basis\_functions\_approach\_to\_GP · GitHub](https://github.com/gabriuma/basis_functions_approach_to_GP/blob/master/Paper/Case-study_Birthday-data/CaseII_Birthday_Approximate-GP.stan). In the paper, I just find the speed per iteration.

Another question if I may, do you think `optimizing()` would work for such a model with huge parameter number? I wonder if I can use it just to estimate the mean value of each parameter. Can I roughly say that the larger the parameter number and noiser the data are, the less accurate estimations `optimizing()` will provide?

---

<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:52am UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/8 "2021-03-31T08:52:26Z")

</div>

> [@zcai](#):
>
> I wonder if I can use it just to estimate the mean value of each parameter.

Optimizing will try to find the MAP parameter estimate, which is not the same as the posterior mean.

> [@zcai](#):
>
> Can I roughly say that the larger the parameter number and noiser the data are, the less accurate estimations `optimizing()` will provide?

Probably. The first problem is that the optimization problem is not convex, so you can’t be sure if the found local optimum is global. And I would think that when the number of parameters grows, it becomes harder and harder to find the global optimum (i.e. the actual MAP parameters).

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [April 2, 2021, 5:52pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/9 "2021-04-02T17:52:38Z")

</div>

One thing to keep in mind with convolution models is that if the convolution is too diffuse then the latent functional behavior can be extremely degenerate – in other words there will be _many_ functional behaviors compatible with any given observation, and the sampler will have to spend an inordinate amount of time exploring those configurations.

When investigating why a model is slow it helps to consider not only the total time (either for the entire model or just parts of the model using the new profiling tools) as well as the sampler behavior, in particular the total number of leapfrog steps and the distribution of the number of leapfrog steps per iteration. If there are a consistent, reasonably small (say less than 25) number of leapfrog steps per iteration then the slow down is probably due to the cost of a single gradient evaluation, but if the numerical trajectories are much longer then the number of gradient evaluations needed per iteration is probably more relevant. See for example some of the diagnostic investigations in [Identity Crisis](https://betanalpha.github.io/assets/case_studies/identifiability.html#5_Don't_Do_Me_Like_That).

---

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [April 2, 2021, 7:05pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/10 "2021-04-02T19:05:39Z")

</div>

> [@betanalpha](#):
>
> One thing to keep in mind with convolution models is that if the convolution is too diffuse then the latent functional behavior can be extremely degenerate – in other words there will be _many_ functional behaviors compatible with any given observation, and the sampler will have to spend an inordinate amount of time exploring those configurations.

Thanks a lot Michael for both replying and the nice case studies. I am not sure I get what you mean by “too diffuse”. If diffuse refers to how dense the 1s are distributed in vector event\_{j}, j =1,2...N\_{obs}, I agree with you that the latent function would capture many functional behaviours. In my simulation and the real applications, it is actually dense enough so that the distance between neighbours is shorter than the length of the latent function f. If these 1s are sparse, we would just extract the echoes and do an average to estimate f. The challenge is they overlay f along x-axis while convolution.

My observations from simulation results also proved your point at some level. Even if I configure 1s to be closer to each other and the distance is randomized to avoid the latent function capturing the periodic behaviour, the fitted GP will still capture some functional behaviour when I give the model a different event sequence that was not used to generate the observation. I understand GP will capture something anyhow, but I was expecting the posterior of marginal deviation (e.g. \alpha in your GP case study) in exponentiated quadratic function would reflect the uncertainty of the estimated latent function so that using the wrong event sequence would end up with larger \alpha. However, it is not the truth, \alpha is larger when giving the wrong sequence but not as large as I expect. The reason I guess so far is sigma\_obs in data likelihood “absorbed” the uncertainty. To solve this I tried to use the square of sigma\_obs in GP kernel instead of a nugget. This time the posterior of \alpha indeed reflected the uncertainty of f much better, but since I add sigma\_obs to the diagonal of the Gram matrix, the estimated function from GP is not smooth anymore.

---

<div class="post-metadata">

**Author:** ![betanalpha](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/betanalpha/32/4_2.png) [@betanalpha](https://discourse.mc-stan.org/u/betanalpha)\
**Post date:** [April 2, 2021, 7:30pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/11 "2021-04-02T19:30:37Z")

</div>

When the configuration of the latent phenomenology (here the gaussian process) and the measurement (here \sigma\_{\text{obs}}) are fit jointly then they too can be degenerate. Indeed these are often the most subtle yet damaging degeneracies. This is especially true when the model is known to be wrong and the likelihood function contorts to capture all of the least-worst model configurations.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [April 7, 2021, 7:21pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/12 "2021-04-07T19:21:44Z")

</div>

> [@zcai](#):
>
> Hi @avehtari, may I ask approximately the computation time of the birthday data?

In [this case study](https://avehtari.github.io/casestudies/Birthdays/birthdays.html)  
from seconds to minutes depending on the model (and it has faster code than the paper).

> [@zcai](#):
>
> Another question if I may, do you think `optimizing()` would work for such a model with huge parameter number?

Optimizing can work if the number of parameters is much smaller than the number of observations. In [the birthday case study](https://avehtari.github.io/casestudies/Birthdays/birthdays.html) optimizing gives quite good result for the simplest model, but not for the others.

---

<div class="post-metadata">

**Author:** ![zcai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/zcai/32/10242_2.png) [@zcai](https://discourse.mc-stan.org/u/zcai)\
**Post date:** [April 7, 2021, 7:45pm UTC](https://discourse.mc-stan.org/t/could-you-please-help-to-find-out-why-my-gp-model-got-slower/21579/13 "2021-04-07T19:45:39Z")

</div>

Thanks, @avehtari, didn’t notice this case study, I will go through it.