Possible ways to speed up spatio-temporal GP modeling

Following the example, I quickly tested HSGP with:

 gp = pm.gp.HSGP(m=[35, 35, 35], c=4.0, mean_func=mu, cov_func=cov_func)
 f = gp.prior("f", X=Xs)

But, I had the following error message:

“NotImplementedError: The power spectral density of products of covariance functions is not implemented.”

I guess I am multiplying two covariances here, but other than that, I don’t have a clue. I would appreciate it if any advice is offered.

Thank you.