From 744cfcfb859f5fafe1c272010d80ec1ea79e33d6 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sat, 26 Sep 2026 09:52:01 +0200 Subject: [PATCH 1/6] Use a symmetric eigensolver for the Laplace block Hessian square root block_matrix_sqrt used RealSchur and rejected any negative entry on the Schur diagonal. The block is a negative Hessian of a log-concave likelihood: symmetric positive semi-definite, and rank deficient whenever there are more latent variables than observations (common with a dense coefficient block). Its zero eigenvalues come out of floating point with either sign, so solver 1 failed on every such call and fell back to solver 2 after computing one Hessian for nothing. Symmetrise the block, use SelfAdjointEigenSolver, clamp eigenvalues that are negative only at rounding level and reject a block with a genuinely negative eigenvalue. The non-finite check also only fired when no entry was finite; it now fires when any entry is not. Co-Authored-By: Claude Fable 5.1 --- .../laplace_marginal_density_estimator.hpp | 69 ++++++++------ .../math/laplace/block_matrix_sqrt_test.cpp | 93 +++++++++++++++++++ 2 files changed, 135 insertions(+), 27 deletions(-) create mode 100644 test/unit/math/laplace/block_matrix_sqrt_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 61b805e4980..300077a9a17 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -10,8 +10,9 @@ #include #include #include -#include +#include #include +#include #include #include @@ -219,57 +220,71 @@ struct laplace_density_estimates { }; /** - * Returns the principal square root of a block diagonal matrix. + * Returns the principal square root of a symmetric positive semi-definite + * block diagonal matrix. + * + * Each block is symmetrised and decomposed with a symmetric eigensolver. + * Eigenvalues that are negative only at rounding level, as the zero + * eigenvalues of a rank-deficient block are (for example the negative + * Hessian of a likelihood with more latent variables than observations), + * are clamped to zero. An eigenvalue below + * `-block_size * epsilon * max(|eigenvalues|, 1)` means the block is not + * positive semi-definite. + * * @tparam WRootMat A type inheriting from `Eigen::EigenBase`. * @param W_root The output matrix to store the square root. * @param W The input block diagonal matrix. * @param block_size The size of each block in the block diagonal matrix. + * @throw std::domain_error if a block has non-finite entries or is not + * positive semi-definite. */ template inline void block_matrix_sqrt(WRootMat& W_root, const Eigen::SparseMatrix& W, const Eigen::Index block_size) { - int n_block = W.cols() / block_size; + const Eigen::Index n_block = W.cols() / block_size; Eigen::MatrixXd local_block(block_size, block_size); Eigen::MatrixXd local_block_sqrt(block_size, block_size); - Eigen::MatrixXd sqrt_t_mat = Eigen::MatrixXd::Zero(block_size, block_size); + Eigen::SelfAdjointEigenSolver eigensolver; // No block operation available for sparse matrices, so we have to loop // See https://eigen.tuxfamily.org/dox/group__TutorialSparse.html#title7 - for (int i = 0; i < n_block; i++) { - sqrt_t_mat.setZero(); + for (Eigen::Index i = 0; i < n_block; i++) { local_block = W.block(i * block_size, i * block_size, block_size, block_size); - if (!local_block.array().isFinite().any()) { + if (!local_block.array().isFinite().all()) { throw std::domain_error( std::string("Error in block_matrix_sqrt: " - "NaNs detected in block diagonal starting at (") + "non-finite values detected in block diagonal " + "starting at (") + std::to_string(i) + ", " + std::to_string(i) + ")"); } - // Issue here, sqrt is done over T of the complex schur - Eigen::RealSchur schurOfA(local_block); - // Compute Schur decomposition of arg - const auto& t_mat = schurOfA.matrixT(); - const auto& u_mat = schurOfA.matrixU(); - // Check if diagonal of schur is not positive - if ((t_mat.diagonal().array() < 0).any()) { + local_block_sqrt = 0.5 * (local_block + local_block.transpose()); + eigensolver.compute(local_block_sqrt); + if (eigensolver.info() != Eigen::Success) { throw std::domain_error( std::string("Error in block_matrix_sqrt: " - "values less than 0 detected in block diagonal's schur " - "decomposition starting at (") + "eigendecomposition failed for block diagonal " + "starting at (") + std::to_string(i) + ", " + std::to_string(i) + ")"); } - try { - // Compute square root of T - Eigen::matrix_sqrt_quasi_triangular(t_mat, sqrt_t_mat); - // Compute square root of arg - local_block_sqrt = u_mat * sqrt_t_mat * u_mat.adjoint(); - } catch (const std::exception& e) { + const Eigen::VectorXd eigenvalues = eigensolver.eigenvalues(); + const double tolerance + = block_size * std::numeric_limits::epsilon() + * std::max(eigenvalues.cwiseAbs().maxCoeff(), 1.0); + if (eigenvalues.minCoeff() < -tolerance) { throw std::domain_error( - "Error in block_matrix_sqrt: " - "The matrix is not positive definite"); + std::string("Error in block_matrix_sqrt: block diagonal starting " + "at (") + + std::to_string(i) + ", " + std::to_string(i) + + ") is not positive semi-definite (smallest eigenvalue " + + std::to_string(eigenvalues.minCoeff()) + ")"); } - for (int k = 0; k < block_size; k++) { - for (int j = 0; j < block_size; j++) { + local_block_sqrt.noalias() + = eigensolver.eigenvectors() + * eigenvalues.cwiseMax(0.0).cwiseSqrt().asDiagonal() + * eigensolver.eigenvectors().transpose(); + for (Eigen::Index k = 0; k < block_size; k++) { + for (Eigen::Index j = 0; j < block_size; j++) { W_root.coeffRef(i * block_size + j, i * block_size + k) = local_block_sqrt(j, k); } diff --git a/test/unit/math/laplace/block_matrix_sqrt_test.cpp b/test/unit/math/laplace/block_matrix_sqrt_test.cpp new file mode 100644 index 00000000000..9fd447fcaf4 --- /dev/null +++ b/test/unit/math/laplace/block_matrix_sqrt_test.cpp @@ -0,0 +1,93 @@ +#include +#include +#include +#include + +namespace { + +// sparse block-diagonal matrix with the block pattern the Laplace solver +// reserves for W_r (see CholeskyWSolverBlock) +Eigen::SparseMatrix block_pattern(int n_blocks, int block_size) { + const int n = n_blocks * block_size; + Eigen::SparseMatrix m(n, n); + m.reserve(Eigen::VectorXi::Constant(n, block_size)); + for (int b = 0; b < n_blocks; ++b) { + for (int k = 0; k < block_size; ++k) { + for (int j = 0; j < block_size; ++j) { + m.insert(b * block_size + j, b * block_size + k) = 1.0; + } + } + } + m.makeCompressed(); + return m; +} + +Eigen::SparseMatrix block_diag( + const std::vector& blocks) { + const int block_size = blocks[0].rows(); + Eigen::SparseMatrix w = block_pattern(blocks.size(), block_size); + for (std::size_t b = 0; b < blocks.size(); ++b) { + for (int k = 0; k < block_size; ++k) { + for (int j = 0; j < block_size; ++j) { + w.coeffRef(b * block_size + j, b * block_size + k) = blocks[b](j, k); + } + } + } + return w; +} + +void expect_principal_sqrt(const Eigen::SparseMatrix& w, + int block_size) { + Eigen::SparseMatrix w_root = block_pattern(w.rows() / block_size, + block_size); + EXPECT_NO_THROW( + stan::math::internal::block_matrix_sqrt(w_root, w, block_size)); + const Eigen::MatrixXd root = w_root; + EXPECT_TRUE(root.isApprox(root.transpose(), 1e-12)); + Eigen::SelfAdjointEigenSolver eig(root); + EXPECT_GE(eig.eigenvalues().minCoeff(), -1e-12); + EXPECT_TRUE((root * root).isApprox(Eigen::MatrixXd(w), 1e-10)); +} + +} // namespace + +TEST(LaplaceBlockMatrixSqrt, PositiveDefiniteBlocks) { + Eigen::MatrixXd m(3, 3); + m << 1.0, 0.3, -0.2, 0.5, 2.0, 0.1, -0.4, 0.2, 1.5; + Eigen::MatrixXd a = m * m.transpose() + Eigen::MatrixXd::Identity(3, 3); + Eigen::MatrixXd b = 2.0 * Eigen::MatrixXd::Identity(3, 3); + expect_principal_sqrt(block_diag({a, b}), 3); +} + +// A negative Hessian with more latent variables than observations is only +// positive semi-definite: its zero eigenvalues come out of floating point +// as tiny values of either sign and must not be rejected. +TEST(LaplaceBlockMatrixSqrt, RankDeficientBlockIsAccepted) { + Eigen::MatrixXd z(3, 6); + z << 1.0, 0.5, -0.3, 0.8, 0.1, -0.6, 0.2, -1.1, 0.4, 0.3, 0.9, 0.7, -0.5, 0.6, + 1.2, -0.2, 0.4, 0.3; + Eigen::MatrixXd w = z.transpose() * z; // rank 3 of 6 + expect_principal_sqrt(block_diag({w}), 6); + + Eigen::MatrixXd ones = Eigen::MatrixXd::Ones(2, 2); // rank 1 of 2 + Eigen::MatrixXd spd = Eigen::MatrixXd::Identity(2, 2); + expect_principal_sqrt(block_diag({spd, ones}), 2); +} + +TEST(LaplaceBlockMatrixSqrt, IndefiniteBlockThrows) { + Eigen::MatrixXd indefinite(2, 2); + indefinite << 1.0, 0.0, 0.0, -1.0; + Eigen::SparseMatrix w = block_diag({indefinite}); + Eigen::SparseMatrix w_root = block_pattern(1, 2); + EXPECT_THROW(stan::math::internal::block_matrix_sqrt(w_root, w, 2), + std::domain_error); +} + +TEST(LaplaceBlockMatrixSqrt, NonFiniteBlockThrows) { + Eigen::MatrixXd nan_block = Eigen::MatrixXd::Identity(2, 2); + nan_block(0, 1) = std::numeric_limits::quiet_NaN(); + Eigen::SparseMatrix w = block_diag({nan_block}); + Eigen::SparseMatrix w_root = block_pattern(1, 2); + EXPECT_THROW(stan::math::internal::block_matrix_sqrt(w_root, w, 2), + std::domain_error); +} From 2b22b75330c619b7aab30bc7c891b5a71f47f897 Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Sat, 26 Sep 2026 09:52:04 +0200 Subject: [PATCH 2/6] Fix Laplace solver fallback reporting log_solver_fallback threw whenever the message stream was null, even with allow_fallthrough set, so callers without a stream lost the fallback entirely. And since the call sits inside std::call_once, a failure with allow_fallthrough unset after an earlier logged fallback was not reported at all and ended in a misleading "solver not valid" error. Decide first: throw with the failing solver's reason when fallthrough is disabled, otherwise log the fallback once if a stream is available. Co-Authored-By: Claude Fable 5.1 --- .../laplace_marginal_density_estimator.hpp | 49 ++++++++---- .../laplace/laplace_solver_fallback_test.cpp | 74 +++++++++++++++++++ 2 files changed, 109 insertions(+), 14 deletions(-) create mode 100644 test/unit/math/laplace/laplace_solver_fallback_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 300077a9a17..2fc7efa36dc 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -1036,8 +1036,25 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, } /** - * @brief Log a solver fallback event to the provided stream. - * @param[in] allow_fallthrough If false, throw instead of logging + * @brief Throw for a solver failure when falling through to the next solver + * is not allowed. + * @param[in] context Context string for the message + * @param[in] iter Current iteration number + * @param[in] failed_solver Name of the solver that failed + * @param[in] e Exception that caused the failure + */ +[[noreturn]] inline void throw_solver_failure(std::string_view context, + Eigen::Index iter, + std::string_view failed_solver, + const std::exception& e) { + std::ostringstream os; + os << context << ": " << failed_solver << " failed at iteration " << iter + << " and allow_fallthrough is false. Reason: " << e.what(); + throw std::domain_error(os.str()); +} + +/** + * @brief Log a solver fallback event to the provided stream, if any. * @param[in,out] msgs Output stream (may be nullptr) * @param[in] context Context string for the log * @param[in] iter Current iteration number @@ -1045,16 +1062,17 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, * @param[in] next_solver Name of the solver being attempted next * @param[in] e Exception that caused the fallback */ -inline void log_solver_fallback(const bool allow_fallthrough, - std::ostream* msgs, std::string_view context, +inline void log_solver_fallback(std::ostream* msgs, std::string_view context, Eigen::Index iter, std::string_view failed_solver, std::string_view next_solver, const std::exception& e) { + if (!msgs) { + return; + } // Build once so we don't interleave with other logs. std::ostringstream os; - std::string msg_type = allow_fallthrough ? "WARNING" : "ERROR"; - os << "[" << context << "] " << msg_type << ": solver fallback\n" + os << "[" << context << "] WARNING: solver fallback\n" << " " << std::left << std::setw(12) << "iteration:" << iter << "\n" << " " << std::left << std::setw(12) << "failed:" << failed_solver << "\n" << " " << std::left << std::setw(12) << "reason:" << e.what() << "\n" @@ -1062,11 +1080,7 @@ inline void log_solver_fallback(const bool allow_fallthrough, << "trying " << next_solver << "\n" << "note: this warning message will only be displayed once." << "\n"; - if (allow_fallthrough && msgs) { - (*msgs) << os.str(); - } else { - throw std::domain_error(std::string("[") + std::string(context) + "]"); - } + (*msgs) << os.str(); } template @@ -1228,13 +1242,16 @@ inline auto laplace_marginal_density_est( const std::string solver_type = (options.hessian_block_size == 1) ? "Diagonal" : "Block"; std::string failed = "solver 1 (" + solver_type + " Hessian-root Cholesky)"; + if (!options.allow_fallthrough) { + throw_solver_failure("laplace_marginal_density", step_iter, failed, e); + } std::call_once( fallback_warning, [](auto&&... args) { log_solver_fallback(std::forward(args)...); }, - options.allow_fallthrough, msgs, "laplace_marginal_density", step_iter, - std::move(failed), "solver 2 (Covariance-root Cholesky)", e); + msgs, "laplace_marginal_density", step_iter, std::move(failed), + "solver 2 (Covariance-root Cholesky)", e); } try { if (options.solver == 2 || options.allow_fallthrough) { @@ -1243,12 +1260,16 @@ inline auto laplace_marginal_density_est( covariance, update_fun, msgs); } } catch (const std::exception& e) { + if (!options.allow_fallthrough) { + throw_solver_failure("laplace_marginal_density", step_iter, + "solver 2 (Covariance-root Cholesky)", e); + } std::call_once( fallback_warning, [](auto&&... args) { log_solver_fallback(std::forward(args)...); }, - options.allow_fallthrough, msgs, "laplace_marginal_density", step_iter, + msgs, "laplace_marginal_density", step_iter, "solver 2 (Covariance-root Cholesky)", "solver 3 (General LU solver)", e); } diff --git a/test/unit/math/laplace/laplace_solver_fallback_test.cpp b/test/unit/math/laplace/laplace_solver_fallback_test.cpp new file mode 100644 index 00000000000..53530cdef01 --- /dev/null +++ b/test/unit/math/laplace/laplace_solver_fallback_test.cpp @@ -0,0 +1,74 @@ +#include +#include +#include +#include +#include +#include + +namespace { + +struct IdentityCovariance2 { + Eigen::MatrixXd operator()(std::ostream* /*msgs*/) const { + return Eigen::MatrixXd::Identity(2, 2); + } +}; + +// Concave in the first coordinate, mildly convex in the second: the negative +// Hessian of the log likelihood has a negative diagonal entry, so solver 1 +// (diagonal Hessian-root Cholesky) rejects it, while the posterior curvature +// 1 - 0.2 stays positive and solver 2 converges. +struct PartlyConvexLikelihood { + template + auto operator()(const Theta& theta, std::ostream* /*msgs*/) const { + return -0.5 * stan::math::square(theta(0) - 1.0) + + 0.1 * stan::math::square(theta(1)); + } +}; + +double run(bool allow_fallthrough, std::ostream* msgs) { + Eigen::VectorXd theta0 = Eigen::VectorXd::Zero(2); + return stan::math::laplace_marginal_tol( + PartlyConvexLikelihood{}, std::tuple<>{}, 1, IdentityCovariance2{}, + std::tuple<>{}, + std::make_tuple(theta0, 1e-10, 100, 1, 100, allow_fallthrough ? 1 : 0), + msgs); +} + +} // namespace + +TEST(LaplaceSolverFallback, FallbackIsLoggedToStream) { + std::ostringstream msgs; + const double result = run(true, &msgs); + EXPECT_TRUE(std::isfinite(result)); + EXPECT_NE(msgs.str().find("solver fallback"), std::string::npos); + EXPECT_NE(msgs.str().find("not positive definite"), std::string::npos); +} + +// A null message stream must not turn an allowed fallback into an error. +TEST(LaplaceSolverFallback, FallbackWithNullStreamDoesNotThrow) { + double result = std::numeric_limits::quiet_NaN(); + EXPECT_NO_THROW(result = run(true, nullptr)); + EXPECT_TRUE(std::isfinite(result)); +} + +// With fallthrough disabled the solver failure is reported as such, also +// after an earlier fallback in the same process has already been logged. +TEST(LaplaceSolverFallback, DisabledFallthroughThrowsWithReason) { + std::ostringstream msgs; + try { + run(false, &msgs); + FAIL() << "expected std::domain_error"; + } catch (const std::domain_error& e) { + const std::string what = e.what(); + EXPECT_NE(what.find("allow_fallthrough"), std::string::npos) << what; + EXPECT_NE(what.find("not positive definite"), std::string::npos) << what; + } + try { + run(false, nullptr); + FAIL() << "expected std::domain_error"; + } catch (const std::domain_error& e) { + EXPECT_NE(std::string(e.what()).find("allow_fallthrough"), + std::string::npos) + << e.what(); + } +} From 8f8e8cdd380da8befb9289c8c9e8daa9f0e30e8a Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 26 Sep 2026 15:19:24 -0400 Subject: [PATCH 3/6] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/mix/functor/laplace_marginal_density_estimator.hpp | 5 ++--- test/unit/math/laplace/block_matrix_sqrt_test.cpp | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index 2fc7efa36dc..21c881d9e2a 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -268,9 +268,8 @@ inline void block_matrix_sqrt(WRootMat& W_root, + std::to_string(i) + ", " + std::to_string(i) + ")"); } const Eigen::VectorXd eigenvalues = eigensolver.eigenvalues(); - const double tolerance - = block_size * std::numeric_limits::epsilon() - * std::max(eigenvalues.cwiseAbs().maxCoeff(), 1.0); + const double tolerance = block_size * std::numeric_limits::epsilon() + * std::max(eigenvalues.cwiseAbs().maxCoeff(), 1.0); if (eigenvalues.minCoeff() < -tolerance) { throw std::domain_error( std::string("Error in block_matrix_sqrt: block diagonal starting " diff --git a/test/unit/math/laplace/block_matrix_sqrt_test.cpp b/test/unit/math/laplace/block_matrix_sqrt_test.cpp index 9fd447fcaf4..36c0a861fa6 100644 --- a/test/unit/math/laplace/block_matrix_sqrt_test.cpp +++ b/test/unit/math/laplace/block_matrix_sqrt_test.cpp @@ -38,8 +38,8 @@ Eigen::SparseMatrix block_diag( void expect_principal_sqrt(const Eigen::SparseMatrix& w, int block_size) { - Eigen::SparseMatrix w_root = block_pattern(w.rows() / block_size, - block_size); + Eigen::SparseMatrix w_root + = block_pattern(w.rows() / block_size, block_size); EXPECT_NO_THROW( stan::math::internal::block_matrix_sqrt(w_root, w, block_size)); const Eigen::MatrixXd root = w_root; From eeb7792e9bc7e3ab4123b11befa72ce936260b61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Mon, 28 Sep 2026 22:09:18 +0200 Subject: [PATCH 4/6] Update stan/math/mix/functor/laplace_marginal_density_estimator.hpp Co-authored-by: Steve Bronder --- stan/math/mix/functor/laplace_marginal_density_estimator.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index 21c881d9e2a..f3bf20d5e66 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -268,8 +268,8 @@ inline void block_matrix_sqrt(WRootMat& W_root, + std::to_string(i) + ", " + std::to_string(i) + ")"); } const Eigen::VectorXd eigenvalues = eigensolver.eigenvalues(); - const double tolerance = block_size * std::numeric_limits::epsilon() - * std::max(eigenvalues.cwiseAbs().maxCoeff(), 1.0); +const double tolerance = block_size * std::numeric_limits::epsilon() + * eigenvalues.cwiseAbs().maxCoeff(); if (eigenvalues.minCoeff() < -tolerance) { throw std::domain_error( std::string("Error in block_matrix_sqrt: block diagonal starting " From 5f22ff1a8a66650fd1ce478df5f8dfeac7c32bca Mon Sep 17 00:00:00 2001 From: "Jachym.Barvinek" Date: Mon, 28 Sep 2026 22:10:44 +0200 Subject: [PATCH 5/6] Match the block_matrix_sqrt documentation and indentation to the relative tolerance Co-Authored-By: Claude Fable 5.1 --- .../math/mix/functor/laplace_marginal_density_estimator.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index f3bf20d5e66..b9bc52ccd06 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -228,7 +228,7 @@ struct laplace_density_estimates { * eigenvalues of a rank-deficient block are (for example the negative * Hessian of a likelihood with more latent variables than observations), * are clamped to zero. An eigenvalue below - * `-block_size * epsilon * max(|eigenvalues|, 1)` means the block is not + * `-block_size * epsilon * max(|eigenvalues|)` means the block is not * positive semi-definite. * * @tparam WRootMat A type inheriting from `Eigen::EigenBase`. @@ -268,8 +268,8 @@ inline void block_matrix_sqrt(WRootMat& W_root, + std::to_string(i) + ", " + std::to_string(i) + ")"); } const Eigen::VectorXd eigenvalues = eigensolver.eigenvalues(); -const double tolerance = block_size * std::numeric_limits::epsilon() - * eigenvalues.cwiseAbs().maxCoeff(); + const double tolerance = block_size * std::numeric_limits::epsilon() + * eigenvalues.cwiseAbs().maxCoeff(); if (eigenvalues.minCoeff() < -tolerance) { throw std::domain_error( std::string("Error in block_matrix_sqrt: block diagonal starting " From 45ee5217e535bd0df13c196c264d734d1aea8806 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 29 Sep 2026 13:06:50 -0400 Subject: [PATCH 6/6] adds two std::once_flag values so that each solver bump can throw twice. Hides throwing conditionals behind cold compiler hint macro --- .../laplace_marginal_density_estimator.hpp | 55 +++++++++++-------- 1 file changed, 31 insertions(+), 24 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index b9bc52ccd06..7d7a4b3a893 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -251,32 +251,38 @@ inline void block_matrix_sqrt(WRootMat& W_root, for (Eigen::Index i = 0; i < n_block; i++) { local_block = W.block(i * block_size, i * block_size, block_size, block_size); - if (!local_block.array().isFinite().all()) { - throw std::domain_error( - std::string("Error in block_matrix_sqrt: " - "non-finite values detected in block diagonal " - "starting at (") - + std::to_string(i) + ", " + std::to_string(i) + ")"); + if (unlikely(!local_block.array().isFinite().all())) { + [](auto i) STAN_COLD_PATH { + throw std::domain_error( + std::string("Error in block_matrix_sqrt: " + "non-finite values detected in block diagonal " + "starting at (") + + std::to_string(i) + ", " + std::to_string(i) + ")"); + }(i); } local_block_sqrt = 0.5 * (local_block + local_block.transpose()); eigensolver.compute(local_block_sqrt); - if (eigensolver.info() != Eigen::Success) { - throw std::domain_error( - std::string("Error in block_matrix_sqrt: " - "eigendecomposition failed for block diagonal " - "starting at (") - + std::to_string(i) + ", " + std::to_string(i) + ")"); + if (unlikely(eigensolver.info() != Eigen::Success)) { + [](auto i) STAN_COLD_PATH { + throw std::domain_error( + std::string("Error in block_matrix_sqrt: " + "eigendecomposition failed for block diagonal " + "starting at (") + + std::to_string(i) + ", " + std::to_string(i) + ")"); + }(i); } const Eigen::VectorXd eigenvalues = eigensolver.eigenvalues(); const double tolerance = block_size * std::numeric_limits::epsilon() * eigenvalues.cwiseAbs().maxCoeff(); - if (eigenvalues.minCoeff() < -tolerance) { - throw std::domain_error( - std::string("Error in block_matrix_sqrt: block diagonal starting " - "at (") - + std::to_string(i) + ", " + std::to_string(i) - + ") is not positive semi-definite (smallest eigenvalue " - + std::to_string(eigenvalues.minCoeff()) + ")"); + if (unlikely(eigenvalues.minCoeff() < -tolerance)) { + [](auto&& i, auto&& eigenvalues) { + throw std::domain_error( + std::string("Error in block_matrix_sqrt: block diagonal starting " + "at (") + + std::to_string(i) + ", " + std::to_string(i) + + ") is not positive semi-definite (smallest eigenvalue " + + std::to_string(eigenvalues.minCoeff()) + ")"); + }(i, eigenvalues); } local_block_sqrt.noalias() = eigensolver.eigenvectors() @@ -1142,7 +1148,8 @@ inline auto create_update_fun(ObjFun&& obj_fun, ThetaGradFun&& theta_grad_f, }; } -static STAN_THREADS_DEF std::once_flag fallback_warning; +static STAN_THREADS_DEF std::once_flag fallback_warning_1_2; +static STAN_THREADS_DEF std::once_flag fallback_warning_2_3; /** * For a latent Gaussian model with hyperparameters phi and * latent variables theta, and observations y, this function computes @@ -1245,8 +1252,8 @@ inline auto laplace_marginal_density_est( throw_solver_failure("laplace_marginal_density", step_iter, failed, e); } std::call_once( - fallback_warning, - [](auto&&... args) { + fallback_warning_1_2, + [](auto&&... args) STAN_COLD_PATH { log_solver_fallback(std::forward(args)...); }, msgs, "laplace_marginal_density", step_iter, std::move(failed), @@ -1264,8 +1271,8 @@ inline auto laplace_marginal_density_est( "solver 2 (Covariance-root Cholesky)", e); } std::call_once( - fallback_warning, - [](auto&&... args) { + fallback_warning_2_3, + [](auto&&... args) STAN_COLD_PATH { log_solver_fallback(std::forward(args)...); }, msgs, "laplace_marginal_density", step_iter,