finishing a covariance function implementation (for GPs)

[Seth Flaxman]

Hi folks. Daniel and I have been hoping to find time to finish up his
implementation of the squared exponential (aka RBF aka exponentiated
quadratic) kernel. Anyone have some time to help us? Basically we just
need to finish the tests--and the reason we haven't done that is I
haven't (yet) got my C++ skills up to speed while Daniel's been busy
with other stuff. Here's the branch:

https://github.com/stan-dev/math/tree/feature/issue-191-gp-squared_exponential

The tests we need are:

test/unit/math/prim/mat/fun/cov_sq_exp_test.cpp
test/unit/math/prim/scal/fun/squared_distance_test.cpp
test/unit/math/fwd/mat/fun/cov_sq_exp_test.cpp
test/unit/math/fwd/scal/fun/squared_distance_test.cpp
test/unit/math/rev/mat/fun/cov_sq_exp_test.cpp
test/unit/math/rev/scal/fun/squared_distance_test.cpp

(Help could certainly consist of someone spending time on a hangout
with me getting me up to speed with the C++ I don't know...)

And attached are some worked out numeric evaluations that we can use
to implement the tests.

Daniel, please correct anything I've gotten wrong.

Seth

[Andrew Gelman]

It would be great if there's someone who could help with this. Bob: what about that new summer student we have?
A

> --
> 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/d/optout.
> <cov_sq_exp_tests.m><derivatives.pdf>

[Bob Carpenter]

Devin doesn't know C++ yet and this one's hairy or Daniel would be
done with already. All these higher-order tests are also what's
bogging down the vectorization.

- Bob

[Aki Vehtari]

On Monday, May 23, 2016 at 5:13:19 AM UTC+3, Bob Carpenter wrote:
> Devin doesn't know C++ yet and this one's hairy or Daniel would be
> done with already. All these higher-order tests are also what's
> bogging down the vectorization.

And my guess is that with enough functions vectorized, the covariance matrix computation could be fully vectorized when written in Stan language so that the expression tree would be smallish and it would be quite fast even without internal C++ implementation of the covariance function...

Aki

[Daniel Lee]

Hi all,

I'm back! The code is complete without any optimization of the gradients. Before it's ready, it needs the tests Seth mentioned.

Seth and everyone else, writing tests isn't about C++ code. It's all about being willing to be careful and do a little bit of work. If all you did was in a document put down inputs and outputs, then coding those are trivial. It's generating enough test cases and the combinations under which we would instantiate the code that matters.

The reason I want this to be well tested first is we will want to make this faster. Once we get into the territory of writing custom derivative code for operators, it is immensely helpful to have existing tests that we know work. It's really hard to get code right the first time.

Seth forwarded the analytic gradients (sorry -- I meant to forward that to the list, but I sat on it). What we should do is in any format you want (but a test is best), write out:

- input. Something like: (x data = [1, 2, 3], sigma data = 1, lambda data = 1), (x var = [1, 2, 3], sigma var = 1, lambda var = 1)

- output. Something like: cov_mat = [3x3 matrix], dcov_mat[1,1] / d_x[1] = ?, dcov_mat[1,2] / d_x[2] = ?, ..., dcov_mat[1,1] / d_sigma, ...

And we have to hit the edge cases. So... sigma = 0. lambda = 0. sigma = inf...

I have a few more high priority things going on. I've tested it a little bit, but we need more tests. If you want to help, just write down sets of input and output and we can code that into tests easily. Would it help if I put in a skeleton test to show you how it's done? I don't know how much more I can do to get those that are interested in just writing down test cases.

Daniel

[Seth Flaxman]

A skeleton would certainly help me. Just to clarify though--what I
sent had both analytic and numeric gradients, didn't it? Happy to add
more edge cases of course.

[Daniel Lee]

Ah, I see them now.

That's a start, but there's a lot more to do. We also need the cases where there is only one input variable x.

Having it in a condensed form would be useful. For the tests, unless we're implementing the analytic gradient as a function, we only need the values. Maybe we should implement that analytic gradient in the test? That would be the most useful.

But there's a lot more like what happens when x1 is length 2 and x2 is length 3? Or the other way around? (I know what "should" happen, but does the code actually do that?... I hope so, but these are things that a quick typo in the code could mess things up.) What happens when x1 is length 0? And perhaps a case where x1 and x2 are fairly large would be useful.

Daniel

[Bob Carpenter]

But the current vectorization of functions doesn't cut
down expression tree size --- it's just a loop on the
inside of the function instead of int he Stan language,
which doesn't matter for efficiency.

The vectorized distributions are another matter --- those
do create just one new node (with an edge to each input --- the
complexity is a multiply-add per edge, and also a virtual function
call for each node).

- Bob

[Aki Vehtari]

On Tuesday, May 24, 2016 at 3:43:20 AM UTC+3, Bob Carpenter wrote:
> > On May 23, 2016, at 4:04 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
> >
> > On Monday, May 23, 2016 at 5:13:19 AM UTC+3, Bob Carpenter wrote:
> >> Devin doesn't know C++ yet and this one's hairy or Daniel would be
> >> done with already. All these higher-order tests are also what's
> >> bogging down the vectorization.
> >
> > And my guess is that with enough functions vectorized, the covariance matrix computation could be fully vectorized when written in Stan language so that the expression tree would be smallish and it would be quite fast even without internal C++ implementation of the covariance function...
>
> But the current vectorization of functions doesn't cut
> down expression tree size --- it's just a loop on the
> inside of the function instead of int he Stan language,
> which doesn't matter for efficiency.

Oh :(
Thanks for pointing this out.

Aki

[Krzysztof Sakrejda]

On Monday, May 23, 2016 at 10:43:54 AM UTC-4, Daniel Lee wrote:
> Hi all,
>
> I'm back! The code is complete without any optimization of the gradients. Before it's ready, it needs the tests Seth mentioned.

Probably silly question but: is this going to side-step the current need to construct the full covariance matrix for estimating a gaussian process?

[Daniel Lee]

No. This will just help construct a covariance matrix.

[Rob Trangucci]

I can lend a hand here. We certainly need more tests for these functions. Would anyone object to a cov_sq_exp(matrix[], real, real) and cov_sq_exp(matrix[], matrix[], real, real) implementations?

Rob

[Seth Flaxman]

That'd be great, Rob. Daniel and I discussed these awhile back, and
they seem perfectly reasonable. Eventually we'll also want to support
a vector of lengthscales, one per dimension for the non-isotropic case
(also known as ARD but I guess that's not the best name). In any case,
I'd say tests for what we have currently are the highest priority.

[Michael Betancourt]

I think what we really want down the road is
a function that returns the Cholesky decomposition
of the covariance matrix directly, and at that point
we could wrap everything up in a single (matrix-valued)
node on the expression graph, no?

The nice feature of this design is that if we get
structured covariance matrices then we can use
structured solvers behind the scenes (once we have
support for structured solvers, etc…)

[Daniel Lee]

+1! Should have thought of that earlier.

I think Aki had some thoughts on naming. I think we have some time before this is ready since it's so light on tests.

Daniel

[Aki Vehtari]

On Thursday, May 26, 2016 at 3:28:54 PM UTC+3, Daniel Lee wrote:
> I think Aki had some thoughts on naming. I think we have some time before this is ready since it's so light on tests.

1) squared exponential vs exponentiated quadratic
- squared exponential is the older name, but Neil Lawrence proposed some years ago that base on the equation it should be called exponentiated quadratic, and this has started to gain popularity also
2) abbreviation could be _se or _eq depending which names is chosen


Aki

[Daniel Lee]

Hi Aki,

Can you be more concrete? Right now, it's named: "cov_sq_exp". What should it be named instead? (Please provide the exact name of what you're proposing the function be called.)

Rob and I looked through the branch this morning. It looks like I was a good citizen and wrote tests months ago! We need forward mode and edge cases, but that's a lot less than I thought.

Daniel

[Bob Carpenter]

Neil Lawrence is right that it makes more sense, so if
there's agreement among those who teach it (he and Aki being
leaders in that endeavour), then by all means lets go with that.
So that'd be one of

cov_exp_quad OR exp_quad_cov

I do not think we should call it cov_eq or eq_cov, because "eq"
is too confusable with "equals" and we've used "exp" and
"quad" and "cov" elsewhere.

- Bob

[Krzysztof Sakrejda]

On Thursday, May 26, 2016 at 2:50:27 PM UTC-4, Bob Carpenter wrote:
> Neil Lawrence is right that it makes more sense, so if
> there's agreement among those who teach it (he and Aki being
> leaders in that endeavour), then by all means lets go with that.
> So that'd be one of
>
> cov_exp_quad OR exp_quad_cov
>
> I do not think we should call it cov_eq or eq_cov, because "eq"
> is too confusable with "equals" and we've used "exp" and
> "quad" and "cov" elsewhere.

+1 The short acronyms are too hard to remember. Krzysztof

[Ben Goodrich]

On Thursday, May 26, 2016 at 2:50:27 PM UTC-4, Bob Carpenter wrote:

So that'd be one of

  cov_exp_quad  OR  exp_quad_cov

 
cov_exp_quad. If we eventually have a few of these covariance functions, we want them to appear contiguously in the index.

Ben

[Michael Betancourt]

+1

[Robert Trangucci]

Sounds like a consensus on cov_exp_quad. I'll make the change on the branch. 

Rob