Until we get MPI and/or GPU-accelerated cholesky decomposition, GPs involving large (>100x100) covariance matrices will be a headache. In these cases I’ve been reverting to splines as an approximation. Demo here: Spline fitting demo inc comparison of sparse vs non-sparse
The last time I ran comparisons, cov_exp_quad was indeed faster than hand-coded computation of the covariance matrix. If you were dealing with GPs on a grid, you could get nearly the same performance as cov_exp_quad by using tricks to eliminate redundant computations, but still not quite as fast.