Gaussian processes
Bob Carpenter
Gaussian processes. I start with a simple Stan model that
samples from the prior, discuss a Cholesky transform for
efficient multivariate Gaussian generation, extend it to obvious joint model
that can compute the posterior, then to hyperparameter estimation,
then I unfold it with the analytical posterior. I also do logistic
regression. So it's quite a few models. I tested them all
on simulated data.
The examples are in src/models/misc/gaussian-process
and I included R code to generate the data and also run
the examples in Rstan.
It all works surprisingly well, at least up to a few
hundred dimensions. And it's pretty easy to code in Stan.
(Seeing Yuanjun's code for Andrew's class really helped!)
One issue I was wondering about is how to code the linear
algebra most effectively. The sources I read (Neal, Rasmussen
and Williams, and MacKay) all do it differently from each other
and from the way I did it.
K: M x N
C: N x N, pos def
y: N x 1
W: M x M
posterior mean: K' * inv(C) * y
posterior covar: W - K' * inv(C) * K
One approach (what you see from Radford), is to compute
inv(C) * y = inv(LL') * y
= L' \ L \ y
I did this the other way around, computing
K' * inv(C) = K / L / L'
so that I could reuse the computation.
Should I be doing it the other way? I wouldn't think
it would really matter much because the Cholesky
factorization should dominate on any sizable problem
and the auto-dif on anything smaller. (These will go
lots faster when I add the partials for matrix division.)
We could add a gaussian_process_sq_exp_log() function
that takes a parameterization of a squared exponential.
Something like this for the prior:
real gaussian_process_sq_exp_log(vector y, vector x,
real eta_sq, real rho_sq, real sigma_sq);
and this for the posterior:
real gaussian_process_sq_exp_log(vector y_tilde, vector x_tilde,
vector y, vector x,
real eta_sq, real rho_sq, real sigma_sq);
And we could do the same for the logistic case, though there's
no good analytic optimization for the posterior. (There are on
the other hand, a whole lot of approximate methods.)
- Bob
Marcus Brubaker
We could add a gaussian_process_sq_exp_log() function
that takes a parameterization of a squared exponential.
Something like this for the prior:
real gaussian_process_sq_exp_log(vector y, vector x,
real eta_sq, real rho_sq, real sigma_sq);
and this for the posterior:
real gaussian_process_sq_exp_log(vector y_tilde, vector x_tilde,
vector y, vector x,
real eta_sq, real rho_sq, real sigma_sq);
Bob Carpenter
On 1/16/13 3:24 PM, Marcus Brubaker wrote:
> On Wed, Jan 16, 2013 at 7:50 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
...
> algebra most effectively. The sources I read (Neal, Rasmussen
> and Williams, and MacKay) all do it differently from each other
> and from the way I did it.
>
> K: M x N
> C: N x N, pos def
> y: N x 1
> W: M x M
>
> posterior mean: K' * inv(C) * y
>
> posterior covar: W - K' * inv(C) * K
>
> One approach (what you see from Radford), is to compute
>
> inv(C) * y = inv(LL') * y
>
> = L' \ L \ y
>
> I did this the other way around, computing
>
> K' * inv(C) = K / L / L'
>
> so that I could reuse the computation.
>
> Should I be doing it the other way? I wouldn't think
> it would really matter much because the Cholesky
> factorization should dominate on any sizable problem
> and the auto-dif on anything smaller. (These will go
> lots faster when I add the partials for matrix division.)
>
>
>
> It shouldn't matter too much. Since K' * inv(C) is needed if you want the covar as well your way is probably faster.
> Sometimes people only want the posterior mean, in which case the other way is significantly faster.
it differently.
> A big optimization from the AD perspective will be implementing the derivatives of log(det(A)) directly.
I was just thinking about the multi-normal and the
basic matrix ops.
>
> We could add a gaussian_process_sq_exp_log() function
> that takes a parameterization of a squared exponential.
> Something like this for the prior:
>
>
> and this for the posterior:
>
> real eta_sq, real rho_sq, real sigma_sq);
>
>
> I'd suggest splitting the kernel away from the posterior computation. E.G.,
>
> matrix sq_exp_kernel(vector x, ...);
>
> And
>
> so that people can easily create new kernels and mix and match them.
which is how I coded it already. Though I didn't add a sq_exp_kernel function.
The point would be to do the derivates through the kernel function.
Do you think this kernel isn't general/popular enough to be worth bothering with?
I'm still at the textbook stage of understanding these things.
- Bob
Ben Goodrich
I just added a 10-page chapter to the manual on coding
Gaussian processes. I start with a simple Stan model that
samples from the prior, discuss a Cholesky transform for
efficient multivariate Gaussian generation, extend it to obvious joint model
that can compute the posterior, then to hyperparameter estimation,
then I unfold it with the analytical posterior. I also do logistic
regression. So it's quite a few models. I tested them all
on simulated data.
One thing I would add is that in the case where the covariance matrix is fixed, then calculating the Cholesky factorization once and using multi_normal_cholesky(mu, L) only gets you part way to optimum speed because although L is a constant, the inverse of L is calculated repeatedly inside multi_normal_cholesky (except not literally due to the triangular solving). To get optimum speed, you would want to invert L once, transpose that, and do something like this in the model block
lp__ <- lp__ - 0.5 * dot_self(tL_inv * (y - mu));
But when the covariance matrix is a function of unknown hyperparameters, it doesn't make any difference whether you do multi_normal() or multi_normal_cholesky(). Supporting --- or perhaps even providing convenience functions that implement --- various covariance or correlation functions would be nice, although some of them involve special functions that are not implemented in Stan yet and may or may not be in Boost yet.
Ben
Marcus Brubaker
- Bob
--
Bob Carpenter
> The squared exponential is, by far, the most popular kernel so making it fast is a good thing. That said, discrete data
> and people working with non-stationary data will require different kernels or will use combinations of kernels. Plus,
> some will want to use ARD, others won't. There's a lot of variety out there in how people use GP kernels.
spelling this out -- I just didn't know whether to
trust them :-)
> That said, do you think there's a lot to be gained in terms of speed by having such a specialized function?
along the same lines as the products that you (Marcus)
unfolded.
After that, the gain will depend on how much we can
lean on the algebra of the whole problem and save
subcomputations. I don't think things like vectorization
are going to matter because those matrix ops are a little
more expensive than virtual function calls and a little
pointer chasing.
- Bob
Bob Carpenter
On 1/16/13 5:57 PM, Ben Goodrich wrote:
> On Wednesday, January 16, 2013 2:50:55 PM UTC-5, Bob Carpenter wrote:
>
> I just added a 10-page chapter to the manual on coding
> Gaussian processes. I start with a simple Stan model that
> samples from the prior, discuss a Cholesky transform for
> efficient multivariate Gaussian generation, extend it to obvious joint model
> that can compute the posterior, then to hyperparameter estimation,
> then I unfold it with the analytical posterior. I also do logistic
> regression. So it's quite a few models. I tested them all
> on simulated data.
>
>
> One thing I would add is that in the case where the covariance matrix is fixed, then calculating the Cholesky
> factorization once and using multi_normal_cholesky(mu, L) only gets you part way to optimum speed because although L is
> a constant, the inverse of L is calculated repeatedly inside multi_normal_cholesky (except not literally due to the
> triangular solving).
out. What I meant by the Cholesky transform is that
we generate
z ~ normal(0,1); // vectorized
then in the derived quantities block, set the posterior
sample:
y <- mu + L * z;
where L is just computed once from the fixed covariance
matrix. At least that's the simple example I used.
> To get optimum speed, you would want to invert L once, transpose that, and do something like this
> in the model block
>
> |
> lp__ <- lp__ - 0.5 * dot_self(tL_inv * (y - mu));
> |
it, but I should mention that.
> But when the covariance matrix is a function of unknown hyperparameters, it doesn't make any difference whether you do
> multi_normal() or multi_normal_cholesky(). Supporting --- or perhaps even providing convenience functions that implement
> --- various covariance or correlation functions would be nice, although some of them involve special functions that are
> not implemented in Stan yet and may or may not be in Boost yet.
it. The simple example just uses straightup multi_normal,
which I think is easiest to understand.
Michael B's been adding special functions that he needed
for the cumulative distribtuion functions and their
gradients. There's a lot of opportunity there for
optimization, too, because when we have a truncated
y ~ Foo(theta) [L,U];
there's often a lot of shared work among the density
Foo, and the two cumulative distribution evals at L and U.
Those optimizations are on the future to-do list.
- Bob
Nigel/essence
I spend all my time using gaussian processes with MCMC.
It is quite surprising where the time is spent. Constructing K and C can be a bottleneck in some cases.
Realising what is dependent and independent on x (the current position) is vital, and storing intermediate vectors of things which are independent of x.
I don;t have time to understand Bob's notation, but my guess is that K is dependent on x and inv(C) * y is independent. So most of the work will be spent in repeatedly contsructing K, little work is done in cholesky factorisation or back substitution, at least for the mean. The variance does need repeated back substitutions.
Of course, to calculate
K' * inv(C) * K
one should calculate
P = L \ K
and then construct
P' * P
There are as many varieties of gaussian processes as there are names. Surrogate, proxy, linear bayes, estimator, radial basis function, etc etc.
There are many flavours of the algorithms. I come more from a linear bayes perspective. I also pay significant attention to a linear regreassion model as a starting point, with second order terms. I use sophisticated techniques to determine the efficient terms in the polynomial.
There are issues such as whether one should include the variance of the regression beta coefficients on the variance.
I actually use Matern correlation functions rather than gaussian. They are more flexible but more time consuming, unless (as i do) you use approximate methods. Taking the derivatives of matern adds extra complexity...
There is another range of choices in the hyperparameters, you can go from one global correlation length to adapting every term in the correlation matrix - but rasmussen discusses these alternatives quite well.
As a point of note, I caculate all the derivatives and use NUTS very successfully.
A lot (all?) of the computation can be easily vectorised. This is where R is excellent, at things like matrix * matrix calculations. I don;t know exactly what underlying libraries they use, but some version of BLAS must be there. A C implementation should use similar libraries.
Nigel/essence
More commonly I have 50 dimensions and 400 data points. Hence constructing K is a bottleneck.
I also use not just one gaussian process, but up to 500 independent gaussian processes, all of which contribute to a likelihood function which is used in MCMC.
I repeat this MCMC calculation maybe 50 times, it takes a few days for the whole process!
Bob Carpenter
On 1/17/13 5:28 PM, Nigel/essence wrote:
> ps the problems I tackle can be up to 1000 dimensions, have up to 5000 data points.
Stan to do in terms of matrix ops and gradients.
A naive implementation would take ~ 1K * 5K * 5K) ops to construct
the covar matrix, then ~ (5K)^3 for the Cholesky factorization.
> More commonly I have 50 dimensions and 400 data points. Hence constructing K is a bottleneck.
were saying -- that whether it's the gradient calc or the
matrix or solving will depend on the exact dimensions of the problem.
- Bob
Bob Carpenter
On 1/17/13 5:21 PM, Nigel/essence wrote:
...
> at things like matrix * matrix calculations. I don;t know exactly what underlying libraries they use, but some
version of BLAS must be there. A C implementation should use similar libraries.
loops because Stan loops are compiled to C++ loops.
Having said that, we're using the Eigen C++ lib, which
indeed does efficient matrix ops and linear algebra.
The big deal for us is vectorizing things that involve
the gradient (that is, involve estimated parameters),
because they can dramatically cut down on virtual function
calls, memory fiddling (we manage our own in an arena-based
arrangement), and the size of the expression graph.
- Bob
Nigel/essence
Maybe we are talking with different ssumptions, but....
I am talking about sampling a single Gaussian process using MCMC.
So if we have nS samples, for the mean we do cholesky factorisation 1 time, we construct the k matrix nS times, and we multiply the K matrix by the stored constant vector nS times.
The K matrix takes nData x dimension operations.
So in total the K matrix contruction dominates, it take O(S x dimension x data points) operations.
That is for the mean. The variance is different, of course, you have to do a lot of back substitutions, but again factorisation is only done once for all samples.
Maybe you are talking about optimising the hyper parameters, where of course you have to factorise each time they change.
It is late and i am tired so I may be wrong, I have not checked against my code, but more likely we are talking about different things. If in doubt I can send some profiling results.
Bob Carpenter
it's easy to see in Stan where the time's going
to go.
There are several different models with different
costs and factorizations. There's the plain
old simulation, hyperparameter estimation,
and joint (or marginal) predictive inference. And
then the same things again with a logistic model for
binary data (or generalized to K-way classification
data).
For Stan, a big cost comes in auto-diff for the
hyperparameter estimation, which we need for HMC/NUTS.
- Bob
Nigel/essence
Hyperarameter estimation is a subject in itself, there are papers giving the complcated formulae for the derivative of the likelihood function (if doing maximum likelihood evaluation), and other papers taking a bayesian approach to estimation.
Take a look at the www.mucm.ac.uk web site. Technical papers 09/05 and 11/03 may be of interest.