From fc1693e2e867156c783882daa485121e6b7da366 Mon Sep 17 00:00:00 2001 From: Florence Bockting Date: Mon, 29 Jun 2026 14:20:53 +0300 Subject: [PATCH 1/4] refactor: support additional return argument dependening on ReturnMeanAndCovCholesky value --- stan/math/mix/functor/laplace_base_rng.hpp | 39 +++++++++++++++++----- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/stan/math/mix/functor/laplace_base_rng.hpp b/stan/math/mix/functor/laplace_base_rng.hpp index db14a894d94..1e746d2273a 100644 --- a/stan/math/mix/functor/laplace_base_rng.hpp +++ b/stan/math/mix/functor/laplace_base_rng.hpp @@ -5,6 +5,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -15,11 +16,19 @@ namespace math { * theta ~ Normal(theta | 0, Sigma(phi, x)) * y ~ pi(y | theta, eta) * - * returns a multivariate normal random variate sampled + * By default, returns a multivariate normal random variate sampled * from the Laplace approximation of p(theta_pred | y, phi, x_pred). + * If `ReturnMeanAndCovCholesky` is true, instead of drawing a sample this + * returns the posterior mean and the Cholesky factor of the posterior + * covariance of that same Laplace approximation, as a + * `std::tuple`. * Note that while the data is observed at x (train_tuple), the new samples * are drawn for covariates x_pred (pred_tuple). * To sample the "original" theta's, set pred_tuple = train_tuple. + * @tparam ReturnMeanAndCovCholesky If false (default), draw and return a + * random variate from the approximate posterior. If true, return a tuple + * containing the posterior mean and the lower-triangular Cholesky factor of + * the posterior covariance instead of a sample. * @tparam LLFunc Type of likelihood function. * @tparam LLArgs Tuple of arguments types of likelihood function. * \laplace_common_template_args @@ -31,13 +40,15 @@ namespace math { * \rng_arg * \msg_arg */ -template >* = nullptr> -inline Eigen::VectorXd laplace_base_rng( - LLFunc&& ll_fun, LLArgs&& ll_args, CovarFun&& covariance_function, - CovarArgs&& covar_args, const laplace_options& options, RNG& rng, - std::ostream* msgs) { +inline auto laplace_base_rng(LLFunc&& ll_fun, LLArgs&& ll_args, + CovarFun&& covariance_function, + CovarArgs&& covar_args, + const laplace_options& options, + RNG& rng, std::ostream* msgs) { Eigen::MatrixXd covariance_train = stan::math::apply( [msgs, &covariance_function](auto&&... args) { return covariance_function(std::forward(args)..., msgs); @@ -51,7 +62,12 @@ inline Eigen::VectorXd laplace_base_rng( = md_est.L.template triangularView().solve( md_est.W_r * covariance_train); Eigen::MatrixXd Sigma = covariance_train - V_dec.transpose() * V_dec; - return multi_normal_rng(std::move(mean_train), std::move(Sigma), rng); + if constexpr (ReturnMeanAndCovCholesky) { + Eigen::MatrixXd Sigma_chol = cholesky_decompose(Sigma); + return std::make_tuple(std::move(mean_train), std::move(Sigma_chol)); + } else { + return multi_normal_rng(std::move(mean_train), std::move(Sigma), rng); + } } else { Eigen::MatrixXd Sigma = covariance_train @@ -60,7 +76,12 @@ inline Eigen::VectorXd laplace_base_rng( - md_est.W_r * md_est.LU.solve(covariance_train * md_est.W_r)) * covariance_train; - return multi_normal_rng(std::move(mean_train), std::move(Sigma), rng); + if constexpr (ReturnMeanAndCovCholesky) { + Eigen::MatrixXd Sigma_chol = cholesky_decompose(Sigma); + return std::make_tuple(std::move(mean_train), std::move(Sigma_chol)); + } else { + return multi_normal_rng(std::move(mean_train), std::move(Sigma), rng); + } } } From 87e29725a9c56d7c46735fcce2457a72b56396af Mon Sep 17 00:00:00 2001 From: Florence Bockting Date: Mon, 29 Jun 2026 14:22:27 +0300 Subject: [PATCH 2/4] feat: add new laplace_latent(_tol)_solve() functions --- stan/math/mix/prob/laplace_latent_solve.hpp | 82 +++++++++++++++++++++ 1 file changed, 82 insertions(+) create mode 100644 stan/math/mix/prob/laplace_latent_solve.hpp diff --git a/stan/math/mix/prob/laplace_latent_solve.hpp b/stan/math/mix/prob/laplace_latent_solve.hpp new file mode 100644 index 00000000000..83a9c33bdc8 --- /dev/null +++ b/stan/math/mix/prob/laplace_latent_solve.hpp @@ -0,0 +1,82 @@ +#ifndef STAN_MATH_MIX_PROB_LAPLACE_LATENT_SOLVE_HPP +#define STAN_MATH_MIX_PROB_LAPLACE_LATENT_SOLVE_HPP + +#include + +namespace stan { +namespace math { + +/** + * In a latent gaussian model, + * + * theta ~ Normal(0, Sigma(phi)) + * y ~ p(y|theta,phi) + * + * returns the posterior mean and Cholesky factor from the Laplace + * approximation to p(theta|y,phi), where the log likelihood is given by L_f. + * @tparam LLFunc Type of likelihood function. + * @tparam LLArgs Tuple of arguments types of likelihood function. + * \laplace_common_template_args + * @param ll_fun Likelihood function. + * @param ll_args Arguments for likelihood function. + * \laplace_common_args + * @param[in] hessian_block_size Block size for the Hessian approximation with + * respect to the latent gaussian variable theta. + * \laplace_options + * \rng_arg + * \msg_arg + */ +template +inline auto laplace_latent_tol_solve(LLFunc&& ll_fun, LLArgs&& ll_args, + int hessian_block_size, + CovarFun&& covariance_function, + CovarArgs&& covar_args, OpsTuple&& ops, + RNG& rng, std::ostream* msgs) { + auto options + = internal::tuple_to_laplace_options(std::forward(ops)); + options.hessian_block_size = hessian_block_size; + return laplace_base_rng( + std::forward(ll_fun), std::forward(ll_args), + std::forward(covariance_function), + std::forward(covar_args), std::move(options), rng, msgs); +} + +/** + * In a latent gaussian model, + * + * theta ~ Normal(0, Sigma(phi)) + * y ~ p(y|theta,phi) + * + * returns the posterior mean and Cholesky factor + * from the Laplace approximation of p(theta | y, phi). + * @tparam LLFunc Type of likelihood function. + * @tparam LLArgs Tuple of arguments types of likelihood function. + * \laplace_common_template_args + * @tparam RNG A valid boost rng type + * @param ll_fun Likelihood function. + * @param ll_args Arguments for likelihood function. + * \laplace_common_args + * @param[in] hessian_block_size Block size for the Hessian approximation with + * respect to the latent gaussian variable theta. + * \rng_arg + * \msg_arg + */ +template +inline auto laplace_latent_solve(LLFunc&& ll_fun, LLArgs&& ll_args, + int hessian_block_size, + CovarFun&& covariance_function, + CovarArgs&& covar_args, RNG& rng, + std::ostream* msgs) { + auto options = laplace_options_default{hessian_block_size}; + return laplace_base_rng( + std::forward(ll_fun), std::forward(ll_args), + std::forward(covariance_function), + std::forward(covar_args), std::move(options), rng, msgs); +} + +} // namespace math +} // namespace stan + +#endif From 5f9fe8f3752bab0c9005cdc45e84f38ccd8213ec Mon Sep 17 00:00:00 2001 From: Florence Bockting Date: Mon, 29 Jun 2026 14:22:56 +0300 Subject: [PATCH 3/4] tests: add new tests for laplace_latent_solve --- .../laplace/laplace_latent_solve_test.cpp | 76 +++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 test/unit/math/laplace/laplace_latent_solve_test.cpp diff --git a/test/unit/math/laplace/laplace_latent_solve_test.cpp b/test/unit/math/laplace/laplace_latent_solve_test.cpp new file mode 100644 index 00000000000..eb7337a98fc --- /dev/null +++ b/test/unit/math/laplace/laplace_latent_solve_test.cpp @@ -0,0 +1,76 @@ +#include +#include +#include + +#include + +#include +#include +#include + +namespace { +struct poisson_log_likelihood { + template + auto operator()(const Theta& theta, const std::vector& y, + std::ostream* pstream) const { + return stan::math::poisson_log_lpmf(y, theta); + } +}; +} // namespace + +TEST_F(laplace_count_two_dim_diag_test, latent_solve_mean_and_cov) { + using stan::math::laplace_latent_solve; + auto [mean_est, chol_est] = laplace_latent_solve( + poisson_log_likelihood{}, std::forward_as_tuple(y), 1, + stan::math::test::diagonal_kernel_functor{}, + std::forward_as_tuple(phi(0), phi(1)), rng, nullptr); + constexpr double tol = 1e-6; + EXPECT_EQ(2, mean_est.size()); + EXPECT_NEAR(theta_root(0), mean_est(0), tol); + EXPECT_NEAR(theta_root(1), mean_est(1), tol); + EXPECT_NEAR(0.0, chol_est(0, 1), 1e-12); // check lower triangular matrix + Eigen::MatrixXd Sigma_est = chol_est * chol_est.transpose(); + EXPECT_NEAR(K_laplace(0, 0), Sigma_est(0, 0), tol); + EXPECT_NEAR(K_laplace(1, 1), Sigma_est(1, 1), tol); + EXPECT_NEAR(K_laplace(0, 1), Sigma_est(0, 1), tol); + EXPECT_NEAR(K_laplace(1, 0), Sigma_est(1, 0), tol); +} + +TEST_F(laplace_count_two_dim_diag_test, latent_tol_solve_mean_and_cov) { + using stan::math::laplace_latent_tol_solve; + constexpr double tolerance = 1e-12; + constexpr int max_num_steps = 1000; + constexpr int hessian_block_size = 1; + constexpr int solver = 1; + constexpr int max_steps_line_search = 0; + auto [mean_est, chol_est] = laplace_latent_tol_solve( + poisson_log_likelihood{}, std::forward_as_tuple(y), hessian_block_size, + stan::math::test::diagonal_kernel_functor{}, + std::forward_as_tuple(phi(0), phi(1)), + std::make_tuple(theta_0, tolerance, max_num_steps, solver, + max_steps_line_search, true), + rng, nullptr); + constexpr double tol = 1e-6; + EXPECT_EQ(2, mean_est.size()); + EXPECT_NEAR(theta_root(0), mean_est(0), tol); + EXPECT_NEAR(theta_root(1), mean_est(1), tol); + EXPECT_NEAR(0.0, chol_est(0, 1), 1e-12); // check lower triangular matrix + Eigen::MatrixXd Sigma_est = chol_est * chol_est.transpose(); + EXPECT_NEAR(K_laplace(0, 0), Sigma_est(0, 0), tol); + EXPECT_NEAR(K_laplace(1, 1), Sigma_est(1, 1), tol); + EXPECT_NEAR(K_laplace(0, 1), Sigma_est(0, 1), tol); + EXPECT_NEAR(K_laplace(1, 0), Sigma_est(1, 0), tol); +} + +TEST_F(laplace_count_two_dim_diag_test, + latent_solve_singular_covariance_throws) { + using stan::math::laplace_latent_solve; + EXPECT_THROW(({ + laplace_latent_solve( + poisson_log_likelihood{}, std::forward_as_tuple(y), 1, + stan::math::test::diagonal_kernel_functor{}, + std::forward_as_tuple(0.0, phi(1)), // singular covariance + rng, nullptr); + }), + std::domain_error); +} From ebff9a78fa02644509cdde1db27f7e47f5c994bb Mon Sep 17 00:00:00 2001 From: Florence Bockting Date: Mon, 29 Jun 2026 14:23:26 +0300 Subject: [PATCH 4/4] chore: add laplace_latent_solve to prob.hpp --- stan/math/mix/prob.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/mix/prob.hpp b/stan/math/mix/prob.hpp index a2fa67c4079..50b3dfdc0b3 100644 --- a/stan/math/mix/prob.hpp +++ b/stan/math/mix/prob.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include