There should be 2 versions: + collapsed which integrates out the optimal q(u) [this is already done, but need to extend to multiple output dimensions] + uncollapsed which supports stochastic optimisation