One bottleneck with the nested Laplace approximation is doing a high-dimensional convex optimization. We don’t quite have the tools for this in Stan yet. I’m using a 500-dimensional system, which arises in the latent gaussian poisson process (see again this conversation) as a testing case. I’ve gone back to only focusing on C++ code and all I’m doing is solving an algebraic equation and computing sensitivities. So far:
- Powell dogleg solver requires ~ 2000 s
- general Newton solver ~ 350 s
- general Newton using precision matrices ~127 s
- specialized Newton solver with analytical gradients ~43 s
Here’s the code for the performance tests.
The precision matrix trick means doing everything in terms of the precision matrix instead of using the covariance matrix. If you look at the code, it makes a lot of sense. I tried rewriting the specialized Newton solver with precision matrices, but here, something slows down the whole process, and I’m not sure why. Also, computing sensitivities with precision matrices turns out to be quite expensive, but in theory, it shouldn’t be.
My guess is these unexpected costs are related to memory management.