Will Tebbutt

Hmm that's unfortunate -- it would be a shame to move to a mutating implementation. Another option would be to create a function `my_add` that defines `+` in the way...

This is an example of one `rrule` (for some function utilised in `kldivergence`) returning a `Diagonal` matrix (a "natural" tangent), and another digging into the PDMat and returning a `NamedTuple`...

This is one way to go. You might also be interested in the more recent [A unifying framework for Gaussian process pseudo-point approximations using power expectation propagation](https://www.jmlr.org/papers/volume18/16-603/16-603.pdf)

> As far as I can tell, this package currently implements a general approximation distribution q, which can be optimized "as a whole", which means as many optimization variables as...

Option 2 would be my preference. So everything other than the `ApproxPosteriorGP` struct, and lines 169-201 of the `sparse_approximations.jl` file would move to ApproximateGPs.jl. We've definitely discussed this elsewhere, but...

Yeah, we could probably do that. As the comment in the GaussianProcesses.jl implementation points out, you would just need to [modify the observation variance](https://www.jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf). Specifically for FITC, I think computing...

Ah, sorry, I guess I'm referring to the contents of equations 21-24b in the linked paper

Stheno.jl is pretty ideal for this -- if you're after a particular integral transform to construct inducing features, you can just add it there, and then you'll get interdomain stuff...

Yeah. I spoke with Alex the other day. We could rope him in to this at some point if that would be helpful.

Hmm it looks like something going wrong in ChainRules / Zygote, but I definitely couldn't tell you what 😬