Sparse GPs of different sizes

Modeling
1

Hi Stan folks,

I’m a big stuck on how to implement my type of model. The basic idea is that I have n independent blocks of data, each of different size. Each block can be modeled as a GP, and I’d like to specify either a covariance kernel for each block. One way to do this would be to use a block diagonal matrix, where the covariance between blocks is zero — but these is memory inefficient and would take a while to build each iteration. The better approach would be to have a covariance kernel for each block — but they’d have to be different sizes, and I’m not sure how this would work. As far as I can tell, it’s not possible to build a ragged array of covariance kernels.

I cannot find good instructions on how to go about setting this up. I’ve seen lots of (mostly older) discussions (e.g. here https://groups.google.com/g/stan-users/c/uxjpTq19AjY) but these seem inconclusive on what the best way to set this up would be. Any advice would be much appreciated!

thank you,
Vince

2

Can you just create n GPs?

3

How many dimensions do you have? Asking, because for 1D, 2D, and 3D you could probably get faster inference using basis function approximation and then handling different sizes of blocks would be easy

4

Thank you for the suggestions! I realized my problem was that I was thinking that the covariance matrices would have to be created in transformed_parameters and stored somehow, not created as needed in the model block. Ultimately, my model mixed very poorly and I think might be too hard to fit with GPs. But, for completeness and future folks that stumble on this with similar problems, below is a rough sketch of what did end up working. For some context, the model was a multivariate normal approximation to a binomial, where the dependency between observations in a block was modeled by the GP.

data {
  int nb; // num blcoks
  int nx;  // total number of obs
  int block_lens[nb]; // for block lengths
  // raw data:
  vector[nx] nd;
  vector[nx] nt;
  // the block start/end indices
  int start[nx];
  int end[nx];
}

// log_x comes from transformed parameters and other data

model {
    // ... 
    int l;
    int u;
    int cl;
    for (i in 1:nb) {
      l = start[i];
      u = end[i];
      bl = block_lens[i];

      matrix[bl,bl] Sigma = gp_exponential_cov(pos[l:u], sigma, lscale);

      vector[bl] tau;
      tau = sqrt(nt[l:u] .* exp(log_x[l:u]) .* exp(log1m_exp(log_x[l:u])));

      nd[l:u] ~ multi_normal(nt[l:u].*exp(log_x[l:u]), quad_form_diag(Sigma, tau));
    }
}