Options for speeding up gaussian process regression

[Nick Menzies]

Hi guys, 

I am looking to get some advice on fitting a large(ish) gaussian process regression. The problem involves fitting a GP at ~6000 locations, with two predictors. 

From some early trials it is clear that to fit a GP to these data in Stan I will need to make some efficiency gains.  Options I have been looking at include...

(a) The reduced rank approach described by Maxwell Joseph here, and previously in a post on this forum.

(b) The kronecker inference approach described by Seth Flaxman here, and previously in a post on this forum (and thanks Seth for your earlier response).

The second of these approaches seems like it could work, as a separable covariance function should not be problematic (as far as I can tell), and I can organize my data to lie on a grid of 100 X 100.  I have been working through the appendix of Seth's paper, but the model takes a long time to evaluate, and longer than the speed results shown in his paper. 

I am copying my code below (for transparency, and in case I have made some obvious mistake), but my first question relates to the missingness pattern. While there is 10 to 20% missing in the data, extending the surface to create a grid produces a large area with no observations (see attached plot). Is it problematic that there are no data to fit the GP in this region?  I have no interest in making inference for this part of the surface, but I am wondering if this is causing problems with the sampler.

My more general question is how people proceed with this kind of problem (ie GPs that are too big to fit with the methods described the Stan manual) -- I have seen discussions on this list suggesting that it may be nice but not necessary to estimate the covariance function parameters simultaneously, and instead used fixed values for these. My overall analysis has similarities to the 'birthday problem' that has featured on Andrew's blog a number of times, in that there is a lot of structure in the data to deal with and the GP described is just step one. For that analysis they used GPstuff -- is that the way to go when the problem becomes sufficiently large?  Anyway, sorry some of these questions are general, but I didnt want to focus on optimizing the details if there are larger issues.

ps - thanks all for providing Stan, and for being so responsive on this list.

Nick

///////////////////////////////////////////////////////////

functions {

  // return (A * B) v where:

  // A is n1 x n1, B = n2 x n2, V = n2 x n1 = reshape(v,n2,n1)

  matrix kron_mvprod(matrix A, matrix B, matrix V) {

    return transpose(A * transpose(B * V));

    }

  // A is a length n1 vector, B is a length n2 vector.

  // Treating them as diagonal matrices, this calculates:

  // v = (A * B + sigma2)ˆ{-1}

  // and returns the n1 x n2 matrix V = reshape(v,n1,n2)

  matrix calculate_eigenvalues(vector A, vector B, int n1, int n2, real sigma2) {

    matrix[n1,n2] e;

    for(i in 1:n1) {

      for(j in 1:n2) {

        e[i,j] = (A[i] * B[j] + sigma2);

      }

    }

    return(e);

  }

}

///////////////////////////////////////////////////////////

data {

  int<lower=1> n1;          // number of locations on dimension 1

  int<lower=1> n2;          // number of locations on dimension 2

  int<lower=0> nmissing;    // number of missing obs

  vector[n1] x1;            // locations on dimension 1

  vector[n2] x2;            // locations on dimension 2

  matrix[n2,n1] y;          // matrix with response values for dimension 1 (cols) & 2 (rows)

  int x1missing[nmissing];  // dimension 1 locations for missing values

  int x2missing[nmissing];  // dimension 2 locations for missing values

}

///////////////////////////////////////////////////////////

transformed data {

  matrix[n1, n1] xd1;

  matrix[n2, n2] xd2;

  for(i in 1:n1) {

    xd1[i, i] = 0;

    for (j in (i+1):n1) {

      xd1[i, j] = -pow(x1[i]-x1[j],2);

      xd1[j, i] = xd1[i, j];

    }

  }

  for (i in 1:n2) {

    xd2[i, i] = 0;

    for (j in (i+1):n2) {

      xd2[i, j] = -pow(x2[i]-x2[j],2);

      xd2[j, i] = xd2[i, j];

    }

  }

}

///////////////////////////////////////////////////////////

parameters {

  real<lower=0> var1; 

  real<lower=0> inv_rho1; 

  real<lower=0> inv_rho2;

  real<lower=0.0001> sigma1;

  vector[nmissing] ymissing;

}

///////////////////////////////////////////////////////////

transformed parameters {

  matrix[n2,n1] ystar;

  ystar = y;

  for(i in 1:nmissing) {

    ystar[x2missing[i],x1missing[i]] = ymissing[i];

  }

}

///////////////////////////////////////////////////////////

model {

  matrix[n1, n1] Sigma1;

  matrix[n2, n2] Sigma2;

  matrix[n1, n1] Q1;

  vector[n1] R1;

  matrix[n2, n2] Q2;

  vector[n2] R2;

  matrix[n2,n1] eigenvalues;

  Sigma1 = var1 * exp(xd1 * inv_rho1);

  for(i in 1:n1)

    Sigma1[i,i] = Sigma1[i,i] + 0.0001;

  Sigma2 = exp(xd2 * inv_rho2);

  for(i in 1:n2)

    Sigma2[i,i] = Sigma2[i,i] + 0.0001;

  Q1 = eigenvectors_sym(Sigma1);

  R1 = eigenvalues_sym(Sigma1);

  Q2 = eigenvectors_sym(Sigma2);

  R2 = eigenvalues_sym(Sigma2);

  eigenvalues = calculate_eigenvalues(R2,R1,n2,n1,sigma1);

  var1     ~ lognormal(0,1);

  inv_rho1 ~ normal(0, 1);

  inv_rho2 ~ normal(0, 1);

  sigma1   ~ lognormal(0,1);

  target +=  -0.5 * sum(ystar .* kron_mvprod(Q1,Q2,

              kron_mvprod(transpose(Q1),transpose(Q2),ystar) ./ eigenvalues)) // calculates -0.5 * y’ (K1 * K2) y

             -0.5 * sum(log(eigenvalues)); // calculates logdet(K1 * K2)

}

[Andrew Gelman]

Hi, I'm hoping some of the actual GP experts chime in here, but my quick thought, in addition to trying out lower dimensional models, is that maybe those 0.0001's could be causing some instability?  I have no idea but they kinda look funny to me.

Andrew

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<Figure_MissingnessPattern.pdf>

[Bob Carpenter]

Which ones? The ones lower-bounding scales or the ones
adding jitter to the covariance matrix?

The lower-bounding of scales shouldn't be necessary in
Stan, or I should say, if it is, you have bigger problems
and need something more like a zero-avoiding prior.

Are there recommendations for how much jitter to add based
on the scale of the covariance matrix? It does seem
odd to just add this small number.

- Bob

[Aki Vehtari]

Hi,

Use of the built-in covariance function cov_exp_quad should make this faster.

Rob Trangucci has an example code for Kronecker product GP using the builtin covariance function, and I guess he'll post the link here when he notices this thread.

Rob was able to get good efficiency with 10000 observations.

Missing observations will slow the sampling, but it can vary how much depending on the covariance matrix.

Instead of regular grid in x space, you could use  the spectral representation of Gaussian processes discussed, e.g. here https://arxiv.org/abs/1611.06740. This approach allows irregularily placed x's, and then you would not need to have those missing values. There are also other basis function approaximations which could be useful in this kind of 2 dimensional case (they usually don't scale well to higher dimensions).

Don't fix the covariance function parameters (Andrew is only one who has suggested that, and we will eventually convince him, too, that it's not a good idea)

GPstuff (Matlab/Octave), GPy (Python) or GPflow (Python+Tensorflow) are specialized software packages for GPs and thus can handle large number of observations, but on the other hand they have limited flexibility how you can use GPs as part of a bigger hierarchical model without loosing speed. If you only have spatial part in your model as shown in code example, then these specialized software will be much faster than Stan. But if you need something which is not yet implemented in them, it might be easier to write the model in Stan.

Aki

[Aki Vehtari]

On Tuesday, February 28, 2017 at 7:02:01 AM UTC+2, Bob Carpenter wrote:

Are there recommendations for how much jitter to add based
on the scale of the covariance matrix?  It does seem
odd to just add this small number.

There are different recommendations. In case of MCMC and non-Gaussian likelihoods Radford Neal recommends quite a lot. In some specific cases and carefully implemented numerically stable computations, it's possible to manage with zero jitter. In some cases a good justification based on the floating point accuracy etc. is to have jitter slightly larger than sqrt(eps). Depending on these case I've used something between 1e-4 and 1e-8 (usually smaller).

Aki

[Rob Trangucci]

Hi guys,

Here's a link to an R script https://github.com/rtrangucci/multi_output_gps/blob/master/gp_2dim_array_dgp.R and a Stan program https://github.com/rtrangucci/multi_output_gps/blob/master/gp_2dim_ncp.stan that runs a 100 X 100 GP in about 45 minutes on my laptop:  I haven't run it on real data yet, but this would be a good test! You could definitely modify the stan code to handle missing observations. I'm not sure how much it would slow down the evaluation though. The way that I would do it is to pass in a vector of y's and a vector of the row indices and the column indices in the 100 X 100 grid that correspond to each observation. I'd change the program like so to allow for missing observations:

functions {

  matrix gp_exp_quad_chol(real[] x, real alpha, real len, real jitter) {

    int dim_x = size(x);

    matrix[dim_x, dim_x] L_K_x;

    matrix[dim_x, dim_x] K_x = cov_exp_quad(x, alpha, len);

    for (n in 1:dim_x)

      K_x[n,n] = K_x[n,n] + jitter;

    L_K_x = cholesky_decompose(K_x);

    return L_K_x;

  }

}

data {

  // Model is for outcomes that have been

  // observed on a 2 dimensional grid

  int<lower=1> dim_x_1; // size of grid dimension 1

  int<lower=1> dim_x_2; // size of grid dimension 2

  int<lower=1> N; // Number of observations

  vector[N] y; // Outcome

  int row_i[N]; // row position of each observation

  int col_j[N]; // col position of each observation

  real x_1[dim_x_1]; // locations - 1st dimension

  real x_2[dim_x_2]; // locations - 2nd dimension

}

parameters {

  real<lower=0> len_scale_x_1; // length-scale - 1st dimension

  real<lower=0> len_scale_x_2; // length-scale - 2nd dimension

  real<lower=0> alpha; // Scale of outcomes

  real<lower=0> sigma; // Scale of noise

  // Standardized latent GP

  matrix[dim_x_1, dim_x_2] y_tilde;

}

model {

  matrix[dim_x_1, dim_x_2] latent_gp; 

  {

    matrix[dim_x_1, dim_x_1] L_K_x_1 = gp_exp_quad_chol(x_1, 1.0, len_scale_x_1, 1e-12);

    matrix[dim_x_2, dim_x_2] L_K_x_2 = gp_exp_quad_chol(x_2, alpha, len_scale_x_2, 1e-12);

    // latent_gp is matrix-normal with among-column covariance K_x_1

    // among-row covariance K_x_2

    

    latent_gp = L_K_x_1 * y_tilde * L_K_x_2';

  }

  // priors

  sigma ~ normal(0, 1);

  len_scale_x_1 ~ gamma(8, 2);

  len_scale_x_2 ~ gamma(8, 2);

  alpha ~ normal(0, 1);

  to_vector(y_tilde) ~ normal(0, 1);

  

  // likelihood

  y ~ normal(latent_gp[row_i, col_j], sigma);

}

Rob

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Nick Menzies]

Great, thanks for the quick feedback (Rob, Aki, Bob, Andrew)!

I will give this a go and report back. 45 minutes is certainly reasonable so if I can get somewhere close to that it will be great. 

Nick

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/kwt4XFh55Ek/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Rob Trangucci]

I should have noted that I got around 500 effective samples for the 45 minutes of runtime (timing is for 4 chains in parallel), which is only about 0.18 effective samples per second. If you need more effective samples you'll need to run the sampler for longer. You'll likely get lower effective samples per second on real data.

Rob

[Andrew Gelman]

But still a good idea to give strong priors on these, I think.

[Michael Betancourt]

Yes, especially on the lengthscale! See, for example, the
informative priors used in the Kronekcer code Rob sent.

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Andrew Gelman]

Happy to be in full agreement with Mike for once!

[andre.p...@googlemail.com]

Dear Nick,

did you consider GP by convolution by David Higdon? See reference here:

http://bobby.gramacy.com/teaching/rsm/lect6.html#112

Andre

[Nick Menzies]

Hi Andre -- thanks for this, and no, no yet. Am trying to understand the slide (I am sure the slide is clear, more my ability to understand it). Is this a case of identifying knots positioned throughout the space (I see the spline motivation), or other reduced rank approach?

 Will look at it in more detail after I have worked through the earlier suggestions.  N

Nick Menzies
nick.m...@gmail.com
404 217 1076

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/kwt4XFh55Ek/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Nick Menzies]

Hi All -- I wanted to provide an update on this effort. I implemented Rob's code (thanks Rob!) and it has mostly been good, the one major challenge has been that I couldn't achieve good sampling -- Rhat always high, n-eff very low. I suspect that this is linked to the large region in the grid where I have no observations. The fix I have come to so far is to add fake observations to fill this part of the grid, each of these set to the mean of the other observations. The attached shows a heatmap of the surface I am fitting, and the additional points in red. I have tried to keep these additional points away from the real observations so they don't bias the estimation, but obviously this is a hack.

Anyway, I wanted to see whether this raises any obvious flags, or if there are other better solutions I might try.

As a side note, are there any approaches to GPs based on taxi-cab distance rather than euclidian distance? My reason for asking is that this would greatly reduce the number of discrete distances involved, which I assume could speed the matrix calculations. 

Thanks again for the help -- Nick

[Aki Vehtari]

Hi,

On Monday, March 13, 2017 at 3:11:50 AM UTC+2, Nick Menzies wrote:

Hi All -- I wanted to provide an update on this effort. I implemented Rob's code (thanks Rob!) and it has mostly been good, the one major challenge has been that I couldn't achieve good sampling -- Rhat always high, n-eff very low. I suspect that this is linked to the large region in the grid where I have no observations.

I expected this could happen. You should probably use inducing point approach (in spatial or spectral space). These require a bit more equations, but not impossible to implement (I can send some links later, but don't have time now).

 

As a side note, are there any approaches to GPs based on taxi-cab distance rather than euclidian distance? My reason for asking is that this would greatly reduce the number of discrete distances involved, which I assume could speed the matrix calculations. 

Yes, but you need to be careful as not all covariance functions combined with taxi-cab distance produce positive-semi-definite covariance matrix. And I don't think you would get a speed advantage.

Aki

 

Thanks again for the help -- Nick

On Friday, 3 March 2017 11:32:00 UTC-5, Nick Menzies wrote:

Hi Andre -- thanks for this, and no, no yet. Am trying to understand the slide (I am sure the slide is clear, more my ability to understand it). Is this a case of identifying knots positioned throughout the space (I see the spline motivation), or other reduced rank approach?

 Will look at it in more detail after I have worked through the earlier suggestions.  N

Nick Menzies
nick.m...@gmail.com
404 217 1076

On Fri, Mar 3, 2017 at 10:12 AM, andre.pfeuffer via Stan users mailing list <stan-...@googlegroups.com> wrote:

Dear Nick,

did you consider GP by convolution by David Higdon? See reference here:

http://bobby.gramacy.com/teaching/rsm/lect6.html#112

Andre

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/kwt4XFh55Ek/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.