Doing the decomposition for every momentum draw is inefficient

Summary:

For N parameters, there's an NxN Cholesky computed every time a new momentum is drawn. In reality we only need to recompute that when the metric changes.

Description:

Check the code here: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp#L54

That llt() only needs computed when z.inv_e_metric_ changes, not every time sample_p is called.

Reproducible Steps:

Sampling a simple model with a large number of parameters should be sufficient.

parameters {
  real x[500];
}

model {
  x ~ normal(0, 1);
}

Should do the trick. Run that model with:

./test sample num_warmup=0 adapt engaged=0 algorithm=hmc metric=dense_e

And compare the time with:

./test sample num_warmup=0 adapt engaged=0 algorithm=hmc metric=diag_e

And the difference should be noticeable once that cholesky is precomputed.

Current Output:

Output is fine, just slow.

Expected Output:

Same output.

Current Version:

v2.21.0

Start at src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp around the llt() call in sample_p. Reproduce the slowdown with the 500-parameter model using the provided dense_e and diag_e commands, then verify that the Cholesky work is only repeated when inv_e_metric_ changes and that output remains the same.