It seems to me that the way we sample momentum at https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp#L61 is not quite true. That should be
z.p = z.inv_e_metric_.llt().matrixU().solve(u); # matrixU instead of matrixL
The current way gives Cov( p ) = (L^-1).(L^-t) = (Lt.L)^-1 while it should be Cov( p ) = (LLt)^-1.
Do I understand correctly? Thanks!
Hmm. I think you are understanding it correctly, but I am not sure I understand what it is supposed to be doing. @betanalpha?
inv_e_metric is the precision of a zero-mean normal distribution. This code is meant to sample from the corresponding normal distribution, which is my understanding of the correct linear algebra.
To sample with precision matrix, we take lower cholesky, then do triangle solve with its transpose. But here the code missed the transpose part I guess. (one of the derivation can be found here: https://math.stackexchange.com/questions/1845132/sample-points-from-a-multivariate-normal-distribution-using-only-the-precision-m/2489142).
I just checked what Eigen does and @fehiepsi is right, it’s an error.
Eigen computes A=LL^T, so A.llt().matrixL().solve(u); computes L^{-1}u and if u\sim N(0,I), then \mathbb{E}(L^{-1}u(L^{-1}u)^T) = L^{-1}L^{-T} = (L^T L)^{-1}\neq A^{-1}.
@betanalpha I can do a bug fix this evening if you need.