From b3eef65eaee95857b7580c08a77b1de98cc69f7c Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 20:51:27 +0300 Subject: [PATCH 1/4] return the final Newton iterate as the Laplace latent mean --- stan/math/mix/functor/laplace_base_rng.hpp | 2 +- .../laplace/laplace_latent_solve_test.cpp | 67 +++++++++++++++++++ 2 files changed, 68 insertions(+), 1 deletion(-) diff --git a/stan/math/mix/functor/laplace_base_rng.hpp b/stan/math/mix/functor/laplace_base_rng.hpp index 0a59b778d78..81ffb17eb4b 100644 --- a/stan/math/mix/functor/laplace_base_rng.hpp +++ b/stan/math/mix/functor/laplace_base_rng.hpp @@ -56,7 +56,7 @@ inline auto laplace_base_rng(LLFunc&& ll_fun, LLArgs&& ll_args, std::forward(covar_args)); auto md_est = internal::laplace_marginal_density_est( ll_fun, std::forward(ll_args), covariance_train, options, msgs); - Eigen::VectorXd mean_train = covariance_train * md_est.theta_grad; + Eigen::VectorXd mean_train = md_est.theta; Eigen::MatrixXd Sigma; if (md_est.solver_used == 1) { Eigen::MatrixXd V_dec diff --git a/test/unit/math/laplace/laplace_latent_solve_test.cpp b/test/unit/math/laplace/laplace_latent_solve_test.cpp index 40a7b8bfe6f..e7ca3d4c990 100644 --- a/test/unit/math/laplace/laplace_latent_solve_test.cpp +++ b/test/unit/math/laplace/laplace_latent_solve_test.cpp @@ -242,3 +242,70 @@ TEST(laplace_latent_tol_rng, solver_2_draw_uses_laplace_covariance) { Eigen::VectorXd draw_ref = stan::math::multi_normal_rng(mean, Sigma, rng_ref); EXPECT_MATRIX_NEAR(draw_ref, draw, 1e-8); } + +namespace { +// Longitudinal part of a nonlinear joint model for tumor size: +// y ~ normal(m, sigma * m) with +// m(t) = bsld * (phi * exp(-d * t) + (1 - phi) * exp(g * t)), +// and 4 latent variables on the log or logit scale of (bsld, d, g, phi). +// The log likelihood is not concave in theta. +struct tumor_size_likelihood { + template + auto operator()(const Theta& theta, const Eigen::VectorXd& y, + const Eigen::VectorXd& t, const Eigen::VectorXd& mu, + double sigma, std::ostream* pstream) const { + using stan::math::exp; + using T = stan::value_type_t; + const T bsld = mu(0) * exp(theta(0)); + const T d = mu(1) * exp(theta(1)); + const T g = mu(2) * exp(theta(2)); + const T phi = stan::math::inv_logit(stan::math::logit(mu(3)) + theta(3)); + Eigen::Matrix m(y.size()); + for (Eigen::Index i = 0; i < y.size(); ++i) { + m(i) = bsld * (phi * exp(-d * t(i)) + (1 - phi) * exp(g * t(i))); + } + return stan::math::normal_lpdf(y, m, stan::math::multiply(sigma, m)); + } +}; + +struct diagonal_covariance { + Eigen::MatrixXd operator()(const Eigen::VectorXd& sd, + std::ostream* pstream) const { + return sd.array().square().matrix().asDiagonal(); + } +}; + +auto tumor_size_solve(double tolerance) { + Eigen::VectorXd y(12); + y << 201.04, 78.30, 48.01, 44.39, 52.91, 50.49, 48.63, 67.29, 52.21, 72.43, + 102.43, 105.06; + Eigen::VectorXd t = Eigen::VectorXd::LinSpaced(12, 0.0, 693.0); + Eigen::VectorXd mu(4); + mu << 117.4, 9.121e-3, 1.913e-3, 0.3781; + Eigen::VectorXd sd(4); + sd << 0.4925, 0.9217, 1.120, 1.299; + const double sigma = 0.09826; + Eigen::VectorXd theta_0(4); + theta_0 << 0.3636, -0.3385, 1.9794, -1.9322; + return stan::math::laplace_latent_solve_tol( + tumor_size_likelihood{}, std::forward_as_tuple(y, t, mu, sigma), 4, + diagonal_covariance{}, std::make_tuple(sd), + std::make_tuple(theta_0, tolerance, 100, 3, 100, false), nullptr); +} +} // namespace + +// The returned mean is the final Newton iterate. It used to be +// K * grad log p(y | theta), which differs from the iterate by +// K * grad log p(theta | y) before exact convergence. Here the posterior is +// much narrower than the prior, and with tolerance 1e-6 the old value was +// 0.53 posterior sd from the mode, while the final iterate is 0.0012 sd from +// it. +TEST(laplace_latent_solve, mean_is_final_iterate_tumor_size) { + auto [mean_ref, chol_ref] = tumor_size_solve(1e-14); + const Eigen::VectorXd post_sd + = (chol_ref * chol_ref.transpose()).diagonal().array().sqrt(); + auto [mean_est, chol_est] = tumor_size_solve(1e-6); + const double max_dist + = ((mean_est - mean_ref).array().abs() / post_sd.array()).maxCoeff(); + EXPECT_LT(max_dist, 0.01); +} From 2d8cec3213323ca67c50d2b662192a9b8e771247 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 20:51:34 +0300 Subject: [PATCH 2/4] fix the solver-3 log determinant for negative LU pivots and report det(B) <= 0 --- .../laplace_marginal_density_estimator.hpp | 22 +++- .../laplace_lu_log_determinant_test.cpp | 113 ++++++++++++++++++ 2 files changed, 131 insertions(+), 4 deletions(-) create mode 100644 test/unit/math/laplace/laplace_lu_log_determinant_test.cpp diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index 7d7a4b3a893..d26139b5e84 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -900,13 +900,27 @@ struct LUSolver { * @brief Compute log determinant from LU factorization. * * @note This uses the diagonal of the combined LU matrix produced by Eigen - * (equivalently the diagonal of U). It does not account for the sign of the - * permutation; callers assume `det(B) > 0` in the Laplace correction. + * (equivalently the diagonal of U). B = I + Sigma * W is not symmetric, so + * with partial pivoting the diagonal of U can have negative entries also + * when `det(B) > 0`. The sign of `det(B)` is the sign of the permutation + * times the signs of the diagonal entries of U. * - * @return Sum of log of the LU diagonal entries. + * @return log(det(B)) = sum(log(abs(diag(U)))) + * @throws std::domain_error if `det(B)` is not positive. Then + * `Sigma^-1 + W` is not positive definite at the final iterate, for example + * at a saddle point, and the Laplace approximation is not defined. */ double compute_log_determinant() const { - return lu.matrixLU().diagonal().array().log().sum(); + const auto& u_diag = lu.matrixLU().diagonal(); + const bool odd_sign = (lu.permutationP().determinant() < 0) + != ((u_diag.array() < 0.0).count() % 2 == 1); + if (odd_sign || (u_diag.array() == 0.0).any()) { + throw std::domain_error( + "laplace_marginal_density: det(I + Sigma * W) is not positive at " + "the final iterate, so the negative Hessian of the log density is " + "not positive definite there (for example at a saddle point)."); + } + return u_diag.array().abs().log().sum(); } /** diff --git a/test/unit/math/laplace/laplace_lu_log_determinant_test.cpp b/test/unit/math/laplace/laplace_lu_log_determinant_test.cpp new file mode 100644 index 00000000000..a96bcbac994 --- /dev/null +++ b/test/unit/math/laplace/laplace_lu_log_determinant_test.cpp @@ -0,0 +1,113 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace { +struct normal_likelihood { + template + auto operator()(const Theta& theta, const Eigen::VectorXd& y, double s, + std::ostream* pstream) const { + return stan::math::normal_lpdf(y, theta, s); + } +}; + +// theta(0) has two Cauchy observations y, which make log p(theta | y) +// bimodal in theta(0); theta(1) has one normal observation z. +struct cauchy_normal_likelihood { + template + auto operator()(const Theta& theta, const Eigen::VectorXd& y, double s, + double z, double sz, std::ostream* pstream) const { + return stan::math::cauchy_lpdf(y, theta(0), s) + + stan::math::normal_lpdf(z, theta(1), sz); + } +}; + +struct fixed_covariance { + Eigen::MatrixXd operator()(const Eigen::MatrixXd& K, + std::ostream* pstream) const { + return K; + } +}; +} // namespace + +// Solver 3 factorizes B = I + K W with partial pivoting. Here the rows are +// swapped and U has a negative diagonal entry although det(B) = 226 > 0. The +// log determinant used to be sum(log(diag(U))), which is NaN. For a normal +// likelihood the Laplace approximation is exact, so all solvers must return +// log N(y | 0, K + s^2 I). +TEST(laplace_marginal_tol, all_solvers_exact_with_negative_lu_pivot) { + Eigen::MatrixXd K(2, 2); + K << 1.0, 1.5, 1.5, 4.0; + Eigen::VectorXd y(2); + y << 0.3, -0.2; + const double s = std::sqrt(0.1); + const Eigen::MatrixXd B + = Eigen::MatrixXd::Identity(2, 2) + K / (s * s); // W = I / s^2 + Eigen::PartialPivLU lu(B); + ASSERT_LT(lu.matrixLU().diagonal().minCoeff(), 0.0); + ASSERT_GT(B.determinant(), 0.0); + const double exact = stan::math::multi_normal_lpdf( + y, Eigen::VectorXd::Zero(2), + (K + s * s * Eigen::MatrixXd::Identity(2, 2)).eval()); + Eigen::VectorXd theta_0 = Eigen::VectorXd::Zero(2); + for (int solver : {1, 2, 3}) { + SCOPED_TRACE("solver = " + std::to_string(solver)); + const double lml = stan::math::laplace_marginal_tol( + normal_likelihood{}, std::forward_as_tuple(y, s), 2, fixed_covariance{}, + std::make_tuple(K), + std::make_tuple(theta_0, 1e-12, 1000, solver, 1000, false), nullptr); + EXPECT_NEAR(exact, lml, 1e-8); + } +} + +// Started near the saddle point (0, 0) of log p(theta | y), the Newton +// iteration of solver 3 converges to it. There Sigma^-1 + W is not positive +// definite and the Laplace approximation is not defined. This must be +// reported as such: laplace_marginal_tol() used to return NaN or a finite +// value, and laplace_latent_solve_tol() threw from cholesky_decompose(). +// Started away from the saddle point, solver 3 finds a maximum. +TEST(laplace_marginal_tol, solver_3_saddle_point_throws) { + Eigen::MatrixXd K(2, 2); + K << 4.0, 1.0, 1.0, 1.0; + Eigen::VectorXd y(2); + y << -3.0, 3.0; + const double s = 0.3; + const double z = 0.0; + const double sz = 0.1; + const std::string reason = "det(I + Sigma * W) is not positive"; + Eigen::VectorXd theta_near(2); + theta_near << 0.01, 3.0; + auto ops_near = std::make_tuple(theta_near, 1e-10, 100, 3, 100, false); + try { + stan::math::laplace_marginal_tol( + cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2, + fixed_covariance{}, std::make_tuple(K), ops_near, nullptr); + FAIL() << "laplace_marginal_tol: expected std::domain_error"; + } catch (const std::domain_error& e) { + EXPECT_NE(std::string(e.what()).find(reason), std::string::npos) + << e.what(); + } + try { + stan::math::laplace_latent_solve_tol( + cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2, + fixed_covariance{}, std::make_tuple(K), ops_near, nullptr); + FAIL() << "laplace_latent_solve_tol: expected std::domain_error"; + } catch (const std::domain_error& e) { + EXPECT_NE(std::string(e.what()).find(reason), std::string::npos) + << e.what(); + } + Eigen::VectorXd theta_away(2); + theta_away << 0.3, 3.0; + auto [mean, chol] = stan::math::laplace_latent_solve_tol( + cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2, + fixed_covariance{}, std::make_tuple(K), + std::make_tuple(theta_away, 1e-10, 100, 3, 100, false), nullptr); + EXPECT_NEAR(2.9385, mean(0), 1e-4); + EXPECT_NEAR(0.0097, mean(1), 1e-4); +} From 24f160e18d147a99f5ea36757f86a2813412bb58 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 20:51:27 +0300 Subject: [PATCH 3/4] return the final Newton iterate as the Laplace latent mean --- stan/math/mix/functor/laplace_base_rng.hpp | 2 +- .../laplace/laplace_latent_solve_test.cpp | 67 +++++++++++++++++++ 2 files changed, 68 insertions(+), 1 deletion(-) diff --git a/stan/math/mix/functor/laplace_base_rng.hpp b/stan/math/mix/functor/laplace_base_rng.hpp index 0a59b778d78..81ffb17eb4b 100644 --- a/stan/math/mix/functor/laplace_base_rng.hpp +++ b/stan/math/mix/functor/laplace_base_rng.hpp @@ -56,7 +56,7 @@ inline auto laplace_base_rng(LLFunc&& ll_fun, LLArgs&& ll_args, std::forward(covar_args)); auto md_est = internal::laplace_marginal_density_est( ll_fun, std::forward(ll_args), covariance_train, options, msgs); - Eigen::VectorXd mean_train = covariance_train * md_est.theta_grad; + Eigen::VectorXd mean_train = md_est.theta; Eigen::MatrixXd Sigma; if (md_est.solver_used == 1) { Eigen::MatrixXd V_dec diff --git a/test/unit/math/laplace/laplace_latent_solve_test.cpp b/test/unit/math/laplace/laplace_latent_solve_test.cpp index 40a7b8bfe6f..e7ca3d4c990 100644 --- a/test/unit/math/laplace/laplace_latent_solve_test.cpp +++ b/test/unit/math/laplace/laplace_latent_solve_test.cpp @@ -242,3 +242,70 @@ TEST(laplace_latent_tol_rng, solver_2_draw_uses_laplace_covariance) { Eigen::VectorXd draw_ref = stan::math::multi_normal_rng(mean, Sigma, rng_ref); EXPECT_MATRIX_NEAR(draw_ref, draw, 1e-8); } + +namespace { +// Longitudinal part of a nonlinear joint model for tumor size: +// y ~ normal(m, sigma * m) with +// m(t) = bsld * (phi * exp(-d * t) + (1 - phi) * exp(g * t)), +// and 4 latent variables on the log or logit scale of (bsld, d, g, phi). +// The log likelihood is not concave in theta. +struct tumor_size_likelihood { + template + auto operator()(const Theta& theta, const Eigen::VectorXd& y, + const Eigen::VectorXd& t, const Eigen::VectorXd& mu, + double sigma, std::ostream* pstream) const { + using stan::math::exp; + using T = stan::value_type_t; + const T bsld = mu(0) * exp(theta(0)); + const T d = mu(1) * exp(theta(1)); + const T g = mu(2) * exp(theta(2)); + const T phi = stan::math::inv_logit(stan::math::logit(mu(3)) + theta(3)); + Eigen::Matrix m(y.size()); + for (Eigen::Index i = 0; i < y.size(); ++i) { + m(i) = bsld * (phi * exp(-d * t(i)) + (1 - phi) * exp(g * t(i))); + } + return stan::math::normal_lpdf(y, m, stan::math::multiply(sigma, m)); + } +}; + +struct diagonal_covariance { + Eigen::MatrixXd operator()(const Eigen::VectorXd& sd, + std::ostream* pstream) const { + return sd.array().square().matrix().asDiagonal(); + } +}; + +auto tumor_size_solve(double tolerance) { + Eigen::VectorXd y(12); + y << 201.04, 78.30, 48.01, 44.39, 52.91, 50.49, 48.63, 67.29, 52.21, 72.43, + 102.43, 105.06; + Eigen::VectorXd t = Eigen::VectorXd::LinSpaced(12, 0.0, 693.0); + Eigen::VectorXd mu(4); + mu << 117.4, 9.121e-3, 1.913e-3, 0.3781; + Eigen::VectorXd sd(4); + sd << 0.4925, 0.9217, 1.120, 1.299; + const double sigma = 0.09826; + Eigen::VectorXd theta_0(4); + theta_0 << 0.3636, -0.3385, 1.9794, -1.9322; + return stan::math::laplace_latent_solve_tol( + tumor_size_likelihood{}, std::forward_as_tuple(y, t, mu, sigma), 4, + diagonal_covariance{}, std::make_tuple(sd), + std::make_tuple(theta_0, tolerance, 100, 3, 100, false), nullptr); +} +} // namespace + +// The returned mean is the final Newton iterate. It used to be +// K * grad log p(y | theta), which differs from the iterate by +// K * grad log p(theta | y) before exact convergence. Here the posterior is +// much narrower than the prior, and with tolerance 1e-6 the old value was +// 0.53 posterior sd from the mode, while the final iterate is 0.0012 sd from +// it. +TEST(laplace_latent_solve, mean_is_final_iterate_tumor_size) { + auto [mean_ref, chol_ref] = tumor_size_solve(1e-14); + const Eigen::VectorXd post_sd + = (chol_ref * chol_ref.transpose()).diagonal().array().sqrt(); + auto [mean_est, chol_est] = tumor_size_solve(1e-6); + const double max_dist + = ((mean_est - mean_ref).array().abs() / post_sd.array()).maxCoeff(); + EXPECT_LT(max_dist, 0.01); +} From 697dce2dbecc24a9f0c5ee196a8df721666f9d7f Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 9 Oct 2026 20:51:34 +0300 Subject: [PATCH 4/4] fix the solver-3 log determinant for negative LU pivots and report det(B) <= 0 --- .../laplace_marginal_density_estimator.hpp | 22 +++- .../laplace_lu_log_determinant_test.cpp | 113 ++++++++++++++++++ 2 files changed, 131 insertions(+), 4 deletions(-) create mode 100644 test/unit/math/laplace/laplace_lu_log_determinant_test.cpp diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index 7d7a4b3a893..d26139b5e84 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -900,13 +900,27 @@ struct LUSolver { * @brief Compute log determinant from LU factorization. * * @note This uses the diagonal of the combined LU matrix produced by Eigen - * (equivalently the diagonal of U). It does not account for the sign of the - * permutation; callers assume `det(B) > 0` in the Laplace correction. + * (equivalently the diagonal of U). B = I + Sigma * W is not symmetric, so + * with partial pivoting the diagonal of U can have negative entries also + * when `det(B) > 0`. The sign of `det(B)` is the sign of the permutation + * times the signs of the diagonal entries of U. * - * @return Sum of log of the LU diagonal entries. + * @return log(det(B)) = sum(log(abs(diag(U)))) + * @throws std::domain_error if `det(B)` is not positive. Then + * `Sigma^-1 + W` is not positive definite at the final iterate, for example + * at a saddle point, and the Laplace approximation is not defined. */ double compute_log_determinant() const { - return lu.matrixLU().diagonal().array().log().sum(); + const auto& u_diag = lu.matrixLU().diagonal(); + const bool odd_sign = (lu.permutationP().determinant() < 0) + != ((u_diag.array() < 0.0).count() % 2 == 1); + if (odd_sign || (u_diag.array() == 0.0).any()) { + throw std::domain_error( + "laplace_marginal_density: det(I + Sigma * W) is not positive at " + "the final iterate, so the negative Hessian of the log density is " + "not positive definite there (for example at a saddle point)."); + } + return u_diag.array().abs().log().sum(); } /** diff --git a/test/unit/math/laplace/laplace_lu_log_determinant_test.cpp b/test/unit/math/laplace/laplace_lu_log_determinant_test.cpp new file mode 100644 index 00000000000..a96bcbac994 --- /dev/null +++ b/test/unit/math/laplace/laplace_lu_log_determinant_test.cpp @@ -0,0 +1,113 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace { +struct normal_likelihood { + template + auto operator()(const Theta& theta, const Eigen::VectorXd& y, double s, + std::ostream* pstream) const { + return stan::math::normal_lpdf(y, theta, s); + } +}; + +// theta(0) has two Cauchy observations y, which make log p(theta | y) +// bimodal in theta(0); theta(1) has one normal observation z. +struct cauchy_normal_likelihood { + template + auto operator()(const Theta& theta, const Eigen::VectorXd& y, double s, + double z, double sz, std::ostream* pstream) const { + return stan::math::cauchy_lpdf(y, theta(0), s) + + stan::math::normal_lpdf(z, theta(1), sz); + } +}; + +struct fixed_covariance { + Eigen::MatrixXd operator()(const Eigen::MatrixXd& K, + std::ostream* pstream) const { + return K; + } +}; +} // namespace + +// Solver 3 factorizes B = I + K W with partial pivoting. Here the rows are +// swapped and U has a negative diagonal entry although det(B) = 226 > 0. The +// log determinant used to be sum(log(diag(U))), which is NaN. For a normal +// likelihood the Laplace approximation is exact, so all solvers must return +// log N(y | 0, K + s^2 I). +TEST(laplace_marginal_tol, all_solvers_exact_with_negative_lu_pivot) { + Eigen::MatrixXd K(2, 2); + K << 1.0, 1.5, 1.5, 4.0; + Eigen::VectorXd y(2); + y << 0.3, -0.2; + const double s = std::sqrt(0.1); + const Eigen::MatrixXd B + = Eigen::MatrixXd::Identity(2, 2) + K / (s * s); // W = I / s^2 + Eigen::PartialPivLU lu(B); + ASSERT_LT(lu.matrixLU().diagonal().minCoeff(), 0.0); + ASSERT_GT(B.determinant(), 0.0); + const double exact = stan::math::multi_normal_lpdf( + y, Eigen::VectorXd::Zero(2), + (K + s * s * Eigen::MatrixXd::Identity(2, 2)).eval()); + Eigen::VectorXd theta_0 = Eigen::VectorXd::Zero(2); + for (int solver : {1, 2, 3}) { + SCOPED_TRACE("solver = " + std::to_string(solver)); + const double lml = stan::math::laplace_marginal_tol( + normal_likelihood{}, std::forward_as_tuple(y, s), 2, fixed_covariance{}, + std::make_tuple(K), + std::make_tuple(theta_0, 1e-12, 1000, solver, 1000, false), nullptr); + EXPECT_NEAR(exact, lml, 1e-8); + } +} + +// Started near the saddle point (0, 0) of log p(theta | y), the Newton +// iteration of solver 3 converges to it. There Sigma^-1 + W is not positive +// definite and the Laplace approximation is not defined. This must be +// reported as such: laplace_marginal_tol() used to return NaN or a finite +// value, and laplace_latent_solve_tol() threw from cholesky_decompose(). +// Started away from the saddle point, solver 3 finds a maximum. +TEST(laplace_marginal_tol, solver_3_saddle_point_throws) { + Eigen::MatrixXd K(2, 2); + K << 4.0, 1.0, 1.0, 1.0; + Eigen::VectorXd y(2); + y << -3.0, 3.0; + const double s = 0.3; + const double z = 0.0; + const double sz = 0.1; + const std::string reason = "det(I + Sigma * W) is not positive"; + Eigen::VectorXd theta_near(2); + theta_near << 0.01, 3.0; + auto ops_near = std::make_tuple(theta_near, 1e-10, 100, 3, 100, false); + try { + stan::math::laplace_marginal_tol( + cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2, + fixed_covariance{}, std::make_tuple(K), ops_near, nullptr); + FAIL() << "laplace_marginal_tol: expected std::domain_error"; + } catch (const std::domain_error& e) { + EXPECT_NE(std::string(e.what()).find(reason), std::string::npos) + << e.what(); + } + try { + stan::math::laplace_latent_solve_tol( + cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2, + fixed_covariance{}, std::make_tuple(K), ops_near, nullptr); + FAIL() << "laplace_latent_solve_tol: expected std::domain_error"; + } catch (const std::domain_error& e) { + EXPECT_NE(std::string(e.what()).find(reason), std::string::npos) + << e.what(); + } + Eigen::VectorXd theta_away(2); + theta_away << 0.3, 3.0; + auto [mean, chol] = stan::math::laplace_latent_solve_tol( + cauchy_normal_likelihood{}, std::forward_as_tuple(y, s, z, sz), 2, + fixed_covariance{}, std::make_tuple(K), + std::make_tuple(theta_away, 1e-10, 100, 3, 100, false), nullptr); + EXPECT_NEAR(2.9385, mean(0), 1e-4); + EXPECT_NEAR(0.0097, mean(1), 1e-4); +}