My package is just using Pymc3 as the only backend (for now, hope to expand it to GPFlow), so it can’t really do anything that Pymc3 can’t. I just implemented the bwengels example there (linear model of coregionalization, I believe?), but you’re right that it’s not immediately obvious how/whether you could use structurally distinct but correlated kernels for the two outputs.