@karlnapf 
Currently my plan is to implement the ElogLik method this week.
In fact, I am almost done and I will soon submit a clean code.

And next week I will focus on the L-BFGS method mentioned in the #1971 and use this bound to do inference following the idea of the KL-method of Nickisch and Rasmussen 2008 for GPs. 
BTW, the original method of Nickisch and Rasmussen uses Newton method without any variational bound and the time complexity is quite high.

It turns out that the implementation is similar to my L-BFGS method for Laplacian.
The reason why I want to implement a complete method is that I can do complete unit test and do some observation.