I implemented a Gaussian process model in Stan for modelling the difference of nonlinear data from two groups, for which only data from one group is observed at each time point.

The model should ideally run on data set of 500~1000 data points. Currently, the sampling of 500 HMC samples takes more than a day, so I am hoping to improve the speed…

I did some benchmarking on simulated data (number of data points J = 100). Surprisingly, Nystr"om approximation was *slower* than the original implementation…

The Nystr"om method approximates the kernel matrix K \in R^{J times J} using a subset of the rows and columns of K:

where C \in R^{J \times M} contains a subset of the columns M < J from K, and G \in R^{M \times M} is a subset of M rows and M columns from K.

Time complexity of the original algorithm for evaluating a Gaussian process sample f is O(J^3) due to Cholesky decomposition.

Nystr"om approximation should have time complexity O(J M^2), which is linear in the number of data points.

Of course, sampling complicates matter…

Implementation is available at:

The original implementation (triggered when M >= J) obtains 500 HMC samples in 30 seconds.

The Nystr"om implementation (M = 0.2 J) obtains 500 HMC samples in 3 minutes.

I ran 4 chains each, at least 4 times…

Does subsetting regressors for Gaussian process regression not play well with Hamiltonian Monte Carlo?

Or did I screw up the implementation of Nystr"om reproduced below?

```
// time complexity is O(J M^2) due to matrix multiplication for B
vector approximate_f(matrix K, vector f_eta, vector u, int J, int M) {
// subset index
int idx[M];
// inverse of subset kernel matrix
matrix[M, M] W;
idx = sort_indices_asc(u)[1:M];
// approximate K as \tilde{K} = C W C^T
// where
// C \in R^{J \times M} is subset of columns
// G \in R^{M \times M} is subset of columns and rows
// W = G^{-1}
W = inverse_spd(K[idx, idx]);
// define B \in R^{J \times M} s.t. \tilde{K} = B B^T
// B = C * cholesky_decompose(W)
// then
// f = B * f_eta
return K[, idx] * (cholesky_decompose(W) * f_eta);
}
```

Environmental information:

```
R version 3.3.1
rstan_2.14.2
```