Skip to content

Do not pre-size Laplace matrices that are reassigned #3464

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jachymb wants to merge 11 commits into
base: develop
Choose a base branch
from
Open

Do not pre-size Laplace matrices that are reassigned #3464

Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
57c11f4
laplace_marginal_density: do not pre-size matrices that are reassigned
Oct 4, 2026
04da9a0
laplace_marginal_density: test that solvers 2 and 3 allocate no unuse…
Oct 4, 2026
a738be1
Merge branch 'stan-dev:develop' into bugfix/laplace-arena-reallocation
jachymb Oct 14, 2026
f47f7c6
Merge branch 'develop' into bugfix/laplace-arena-reallocation
Oct 14, 2026
15a2b49
laplace arena test: measure with approx_bytes_used, run solver 3 for …
Oct 14, 2026
900c752
laplace_marginal_density: size R with a ternary, declare LU_solve_cov…
Oct 14, 2026
2e4faf3
Merge remote-tracking branch 'origin/bugfix/laplace-arena-reallocatio…
Oct 14, 2026
d45cfa5
laplace_marginal_density: keep LU_solve_covariance in the enclosing s…
Oct 15, 2026
843d77d
laplace_marginal_density: pre-size R for every solver and write it in…
Oct 15, 2026
9f7da5f
laplace_marginal_density: name the latent dimension once
Oct 15, 2026
9a2fbd0
laplace_marginal_density: remove the arena allocation test
Oct 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 15 additions & 18 deletions stan/math/mix/functor/laplace_marginal_density.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,17 +197,13 @@ inline auto laplace_marginal_density(LLFun&& ll_fun, LLTupleArgs&& ll_args,
ll_fun, ll_args_vals, covariance_val, options, msgs);
auto ll_args_filter = internal::filter_var_scalar_types(ll_args_copy);
// tuple of references to var types
// Solver 1, 2
const bool solver_1_or_2
= md_est.solver_used == 1 || md_est.solver_used == 2;
arena_t<Eigen::MatrixXd> R(md_est.theta.size() * solver_1_or_2,
md_est.theta.size() * solver_1_or_2);
// Solver 3
arena_t<Eigen::MatrixXd> LU_solve_covariance(
covariance.rows() * (md_est.solver_used == 3),
covariance.cols() * (md_est.solver_used == 3));
const Eigen::Index theta_size = md_est.theta.size();
// Every solver writes R in place; no right-hand side below aliases it
arena_t<Eigen::MatrixXd> R(theta_size, theta_size);
// Solver 3 assigns it
arena_t<Eigen::MatrixXd> LU_solve_covariance;
// Solver 1, 2, 3
arena_t<Eigen::VectorXd> s2(md_est.theta.size());
arena_t<Eigen::VectorXd> s2(theta_size);
using stan::math::internal::ZeroOut;
if (md_est.solver_used == 1) {
if (options.hessian_block_size == 1) {
Expand Down Expand Up @@ -247,13 +243,14 @@ inline auto laplace_marginal_density(LLFun&& ll_fun, LLTupleArgs&& ll_args,
internal::copy_compute_s2<ZeroOut>(partial_parm, ll_args_filter);
}
} else if (md_est.solver_used == 2) {
R = md_est.W_r
- md_est.W_r * md_est.K_root
* md_est.L.transpose()
.template triangularView<Eigen::Upper>()
.solve(
md_est.L.template triangularView<Eigen::Lower>()
.solve(md_est.K_root.transpose() * md_est.W_r));
R.deep_copy(md_est.W_r);
R.noalias()
-= md_est.W_r * md_est.K_root
* md_est.L.transpose()
.template triangularView<Eigen::Upper>()
.solve(
md_est.L.template triangularView<Eigen::Lower>().solve(
md_est.K_root.transpose() * md_est.W_r));

arena_t<Eigen::MatrixXd> C
= md_est.L.template triangularView<Eigen::Lower>().solve(
Expand All @@ -268,7 +265,7 @@ inline auto laplace_marginal_density(LLFun&& ll_fun, LLTupleArgs&& ll_args,
auto I_minus_BinvKW
= Eigen::MatrixXd::Identity(md_est.W_r.rows(), md_est.W_r.cols())
- LU_solve_covariance * md_est.W_r;
R = md_est.W_r * I_minus_BinvKW; // == W - W B^{-1} K W
R.noalias() = md_est.W_r * I_minus_BinvKW; // == W - W B^{-1} K W
arena_t<Eigen::MatrixXd> A
= covariance_val - covariance_val * md_est.W_r * LU_solve_covariance;
auto s2_tmp = laplace_likelihood::compute_s2(ll_fun, md_est.theta, A,
Expand Down
Loading