This would be a great fit for what you need: Kronecker Structured Covariances — PyMC3 3.11.5 documentation
But, the number of data points in t might be uncomfortably large. There might be something you could do with inducing points in that dimension, but I haven’t coded it up before. Feel free to post back here if it’s still too slow!