# Improve NUT sampling rate in Gaussian Process when using custom Covariance

**URL:** <https://discourse.pymc.io/t/improve-nut-sampling-rate-in-gaussian-process-when-using-custom-covariance/12073>\
**Category:** v5\
**Tags:** gaussian\_process\
**Created:** [May 7, 2023, 7:07am UTC](https://discourse.pymc.io/t/improve-nut-sampling-rate-in-gaussian-process-when-using-custom-covariance/12073 "2023-05-07T07:07:28Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![TPK](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tpk/32/6355_2.png) [@TPK](https://discourse.pymc.io/u/TPK)\
**Post date:** [May 7, 2023, 7:07am UTC](https://discourse.pymc.io/t/improve-nut-sampling-rate-in-gaussian-process-when-using-custom-covariance/12073/1 "2023-05-07T07:07:28Z")

</div>

I am trying to implement Hierarchical Gaussian Process based on GPy implementation [here](https://discourse.pymc.io/t/how-to-implement-hierarchical-gaussian-process-in-pymc/12056/2). The implementation works, but fitting to the data (NUTS Sampling) is a hit-or-miss i.e. Some data I can sample within a few minutes, but others it took hours, even days to finish. How can I improve the performance

---

<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:** [May 8, 2023, 5:33am UTC](https://discourse.pymc.io/t/improve-nut-sampling-rate-in-gaussian-process-when-using-custom-covariance/12073/2 "2023-05-08T05:33:05Z")

</div>

I think this looks like a perfect situation to try out the new HSGP approximation. I made a gist [here](https://gist.github.com/bwengals/90e79d64fe9e2d496a048f0b2c6d33da) implementing it for a similar data set.

There are a few differences between this and the GPy implementation, but the biggest one is there’s no need to form a big blocked covariance matrix, we can implement the hierarchical GP directly. Sampling should be much much faster too. If you have a lot of GP replicates (`f`), there’s a better implementation now that there is a sparse matrix dot product (credit @ricardoV94 and @tcapretto). I’ll add that to the gist soon.

Just to make sure, how many GP replicates do you have, and how big is each GP? And are your GPs 1 or 2 dimensional?

---

<div class="post-metadata">

**Author:** ![TPK](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tpk/32/6355_2.png) [@TPK](https://discourse.pymc.io/u/TPK)\
**Post date:** [May 8, 2023, 6:02am UTC](https://discourse.pymc.io/t/improve-nut-sampling-rate-in-gaussian-process-when-using-custom-covariance/12073/3 "2023-05-08T06:02:14Z")

</div>

Thanks! for the reply!  
I’m not familiar with HSGP, and I think I have to read up on them first. The data I am trying to model don’t require too many replicates (from `3` to `10`). And each one will have about `100` to `200` datapoints