That looks really great! Definitely seems to do a good job. I for one would like to see this added to the GP examples notebook. Way too slow though! I see the same speed issue when I ran your notebook.
You may have noticed the work being done by @jordan-melendez over in another thread. One way to get a speed boost here is to use the MatrixNormal model he’s working on. Instead of modeling the full covariance matrix (which is 64^2 x 64^2 !), you could have a separate covariance over the x and y dimensions – one 64 x 64 covariance matrix per dimension. This would be much much faster. What I’m thinking of is in this paper, in particular sec. 4.7.