Gaussian Processes in Stan

[Marcus Brubaker]

Hey everyone,

I've been implementing a couple GP models in Stan recently and I've
come up with two issues I wanted to run past the list.

First, there is a need for a "GP distribution" in order to more
efficiently support multivariate GP outputs. Basically, if you have a
vector valued function, people typically treat each output as
independent but with related hyperparameters (and thus, often the same
kernel matrix). Right now in Stan this works out looking something
like:

vector[N] y_trans[d];
vector[d] w;
matrix[N,N] K_sqExp;
matrix[N,N] Sigma;
...
// Output likelihood
for (i in 1:yDim) {
// Compute the kernel matrix for this output
Sigma <- w[i] * K_sqExp;
y_trans[i] ~ multi_normal(mu,Sigma);
}

where y_trans[i] is a vector containing the ith dimension of output
for each training point. This is one of the common cases, where each
output dimension is given a different scale (w[i]) but the underlying
kernel matrix (K_sqExp) is shared.

I'm almost done implementing a multigp_log(y,Sigma,w) function which
would replace the above with:

matrix[d,N] y;
matrix[N,N] K_sqExp;
vector[d] w;
y ~ multigp(K_sqExp,w);

Any comments or suggestions on this? It might be more broadly useful
and, thus, should have a different name perhaps? Should it be
generalized to handle parameter/kernel sharing types?


The other issue is that of repeated matrix factorizations in the
mdivide_left and log_determinant varis. I'd like to create a single
"factorization" vari which can be used to solve linear systems,
compute determinants, etc and only ever factorize the matrix once.
Before that though, I'd just like to avoid repeatedly factorizing a
matrix at construction and in chain(). In principle this shouldn't be
hard, just stick an Eigen:LDLT<> object into the vari and we're done,
except it will never be deleted because of the vari memory management.
I have two proposals for this.

1) Call some kind of vari cleanup function. E.G., every vari has a
virtual cleanup() function which is called when memory is freed.
This is simple to do, but will incur an overhead which will be rarely
needed since the cleanup() function will probably empty 99% of the
time.

2) Alternatively, when a vari is constructed that needs to be cleaned
up, it is added to a special list which is traversed at cleanup time.
This has the advantage of almost zero overhead, but increased
implementation complexity.

Any thoughts on these two approaches? Is there an obvious alternative
that I'm missing?

I'm hoping to have most of this landed in the next few days unless
someone objects.

Cheers,
Marcus

[Bob Carpenter]

For timing, if it goes in independently and passes tests,
no big deal one way or the other. Michael's in the middle
of refactoring the samplers, but those are orthogonal to
new functions and reimplemented vari.

1. I have no opinion about which GP components are usefully
reusable.

As far as coding it, we're working on more general vector views.
But however you code it, we can translate it into the new framework
when we get there.

2. As to saving the results of a solver, by all means. Ben's
been wanting to do something like this for ages.

I prefer a separate stack of things to free. Otherwise, it's
way too much overhead on top of auto-diff.

A virtual cleanup function is the elegant way to do this. It'll
have a bit of extra overhead from the virtual function call over
just saving a stack of pointers to free. Do we need the extra
generality? Or is it worth it just to have the interfaces be
cleaner?

I don't think we'll be able to finagle Eigen's solvers to use
our memory, for instance by using a placement new, because
they're bound to be doing lower-level allocation for internal
matrices.

- Bob

[Marcus Brubaker]

Ok, I've pushed some changes. I've introduced a class called
chainable_alloc with a virtual destructor. When a chainable_alloc is
created, a pointer to it is put on a new stack which is then deleted
in recover_memory(). One important note is that a chainable_alloc
object should *never* be created on the stack, only ever on the heap.

Before doing the full factorization function, I instead implemented a
mdivide_left_spd() function for solving symmetric positive definite
systems and then made use of it in multi_normal_log(). This function
should probably also be used in the wishart and multi_student_t
distributions.

Next up, I'll finish implementing the multivariate GP function and
then I'll probably do the full factorization function for use in both
that and in multi_normal_log.

Cheers,
Marcus

> --
> You received this message because you are subscribed to the Google Groups
> "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Bob Carpenter]

Thanks for filling us in.

And even more thanks for figuring out a general coding standard
that'll let us have more complex vari instances!

- Bob