Multiple Gaussian processes without for loop?

Hmm, I don’t think that’s quite what I’m looking for. That shows multi-dimensional inputs, but my inputs are in one dimension, with multiple outputs.

I also know that PyMC implements Kronecker-structure covariance, but I don’t think I can have multiple length scales in that case (like the MWE).