Kernel creation for MultiVariateNormal

Hi there, I have 2D output Gaussian Processes Regression and I am using the similar approach suggested in Stan documentation for creating the Kernel over two separate dimensions that I have which is as follow (y is my observation):


data{
 int<lower=1> n_buckets;
 int<lower=1> K;
 matrix[K,n_buckets] y; // Outcome
 real<lower=0, upper=1.5> buck_rep[n_buckets];
}
model{
  matrix[K, n_buckets] latent_gp; 
  {

    matrix[n_buckets, n_buckets] L_K_x;
    matrix[n_buckets, n_buckets] L_x = cov_exp_quad(buck_rep, alpha, 0.1);
    for (n in 1:n_buckets){
      L_x[n,n] = L_x[n,n] + 1e-7;
    }
    L_K_x = cholesky_decompose(L_x);
    latent_gp = L_K_y * y_tilde * L_K_x';
    
  }
  L_K_y ~ lkj_corr_cholesky(2);
  to_vector(y_tilde) ~ normal(0, 1);
  sigma ~ normal(0, 1);
  alpha ~ normal(0, 1);
  to_vector(y) ~ normal(to_vector(latent_gp), sigma);  <== ???????????
}

I would like to directly use \text{multi_normal_cholesky}(0, L) instead of
\text{normal}(\text{to_vector}(\text{latent_gp}), sigma).

So, I need your guidance on how I should create L from L_K_y and L_K_x?

If you ask why I am doing this, well I believe with this replacement, I should see some run-time improvement.

Thanks for your help

latent_gp = L_K_y * y_tilde * L_K_x';

The way this model is written makes me think latent_gp is a non-centered parameterization for a certain GP. You could write that using multi_normal_cholesky and that’ll be a centered parameterization which, depending on the data, may sample more or less efficiently (with not much data the non-centered parameterization is usually expected to perform better).

What example in the manual is this based off of (and I’ll take a look)?

Thanks for looking at my quation @bbbales2. I got this model implemented based on the last example " Multiple-output Gaussian processes" here.

Ooof, all I got is a quick answer for you, so the centered version of the model you want is the stuff that appears after the text:

then our finite dimensional generative model for the above is:

The problem is there’s a matrix normal distribution in there, which isn’t one of the distributions in Stan. So to do this, you’ll either need to code up a matrix normal yourself, or you can do a bit of manipulation and fit this into a multivariate normal density.

The equation on the Wikipedia matrix normal page is this one.

I believe vec can be done with to_vector, but double check that, and you’ll have to implement the Kronecker product yourself (though this has been done before so if you search the forums you should be able to find code snippets).

The thing we’d be testing here is whether or not the centered or non-centered parameterizations are faster. This is an easy thing to test with a 1D output since the GPs are built in. My recommendation is try that first. If centered isn’t better for the 1D problem, then I wouldn’t worry about doing all the work to test it for the 2D problem.