Compatibility with Stheno DerivedGPs (or similar)?

[brendanwallison]

Hello, I love this project and your related Juliacon presentation. I was playing around and notice I am able to construct GPODEFunctions from PosteriorGPs, but the moment that I wrap that posteriorGP in a Stheno atomic() constructor I am no longer able to make a GPODEFunction. This is in spite of the fact that at this point it is identical to the original GP, just with an additional bookeeping GPC() object. Would there be much required to build in Stheno compatibility?

For context, my original intention was to use your pullEuler solver on a Stheno DerivedGP to get a feel for how uncertainty propagates. Alternately, I wonder if perhaps the functionality already exists in some manner not requiring Stheno and I'm missing it. Specifically, is there a way to use this package in the context of universal differential equations where the GP supplies some of the terms but other terms are mechanistic? This could include non-additive interactions, for instance dx = (a + NN(x))*x. I was using Stheno for these, but maybe there is a way to implement this using just AbstractGPs in a way that plays well with your work.

[Crown421]

Hey, very glad that you like this project. It is still very much work in progress, and for now I just did some naive implementation, as a few interfaces with AbstractGPs already need wrappers to make them work.

I will have a look at Stheno, and what might be missing.

Your context is one that I am also very interested in, but have not had time to explore. For additional terms, I think that this should be done via the GP mean function, so instead of a zero mean, use a custom mean function of the existing terms. That needs to be researched and considered in more detail though.

The example you give is much more interesting, because in an expression like $a+GP(x))*x$ you now get the product of $GP(x)$ and $x$, both of which will be random variables (in the PULL solver setting). The issue here is that the product of random variables becomes immediately very complicated, even in the case where both factors are normally distributed. I would be happy to have a chat about this if you would like.

[brendanwallison]

In the multiplicative case, I didn't think about the implications of x becoming a random variable itself as you step through the ODE. So this is less a compatibility issue than an open research question.

I was referencing this type of DerivedGP in Stheno: https://juliagaussianprocesses.github.io/Stheno.jl/stable/internals/#Multiplication

But as they note it only works if x is not itself a random variable, which becomes untrue after a single timestep.

For the simpler additive case implemented via a non-zero mean function, you are thinking this is still tricky? In other words, though I haven't tried this yet, let's say I provide an abstractGP with non-zero mean to the GPODEFunc and the code happens to work--I should not blindly trust the result.

[Crown421]

I think the additive case with the custom mean function is the correct way, but I have no evidence for this, just some intuition. I think that should give some decent results, and if you try it I would be very curious.

I will also be adding a sampling -based method as a (much) more computationally expensive, but probably more reliable, way way to get the empirical distribution of the trajectories, that can also be used for comparison.

[brendanwallison]

Definitely, I will try to make some time to explore and share whatever I find with you.

A more rigorous sampling method/baseline also sounds great, particularly for more complex cases like the multiplicative scenario. Looking forward to it!

Thanks again for all the interesting work you're putting into the package.

[brendanwallison]

Hey, just a quick follow-up/confirmation: is the pullEuler solver just for single-output GPs at the moment? I noticed that your spiralODE example, which is multi-output, calls outdated code that just propagates the mean. On the other hand when I follow the gpode example instead I get an error about the initial condition: "MethodError: Cannot convert an object of type Vector{Float64} to an object of type Measurement{Float64}". This seems pretty clear-cut, but I wanted to double-check in case I'm mistaken.

[Crown421]

Hey, this current version only works for 1D (as multi-dim is not quite so easy as I initially thought, essentially order of multiplication "suddenly" matters).

I am right at this moment working out what will hopefully be the last issues, then add tests and update the example.

[Crown421]

@brendanwallison My apologies, it took a bit longer than expected. The PULLEuler now also works for higher dimensions, and the sampling methods are there, but not properly documented yet.