From 552737c1db3acbdde9e606cad514c25b77fcf3fc Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 19 Mar 2026 14:49:41 +0200 Subject: [PATCH 1/7] fix laplace_marginal optimization --- .../laplace_marginal_density_estimator.hpp | 91 +++++++++++++++++-- stan/math/mix/functor/wolfe_line_search.hpp | 28 +++++- 2 files changed, 109 insertions(+), 10 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index 541dbe96fcc..f57ce942594 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -380,6 +380,31 @@ struct NewtonState { wolfe_status.num_backtracks_ = -1; // Safe initial value for BB step } + /** + * @brief Constructs Newton state with a consistent (a_init, theta_init) pair. + * + * When the caller supplies a non-zero theta_init, a_init = Sigma^{-1} * + * theta_init must be provided to maintain the invariant theta = Sigma * a. + * + * @param theta_size Dimension of the latent space + * @param obj_fun Objective function: (a, theta) -> double + * @param theta_grad_f Gradient function: theta -> grad + * @param a_init Initial a value consistent with theta_init + * @param theta_init Initial theta value + */ + template + NewtonState(int theta_size, ObjFun&& obj_fun, ThetaGradFun&& theta_grad_f, + const Eigen::VectorXd& a_init, + ThetaInitializer&& theta_init) + : wolfe_info(std::forward(obj_fun), a_init, + std::forward(theta_init), + std::forward(theta_grad_f), 0), + b(theta_size), + B(theta_size, theta_size), + prev_g(theta_size) { + wolfe_status.num_backtracks_ = -1; // Safe initial value for BB step + } + /** * @brief Access the current step state (mutable). * @return Reference to current WolfeStep @@ -426,9 +451,13 @@ inline void llt_with_jitter(LLT& llt_B, B_t& B, double min_jitter = 1e-10, double max_jitter = 1e-5) { llt_B.compute(B); if (llt_B.info() != Eigen::Success) { + double prev_jitter = 0.0; double jitter_try = min_jitter; for (; jitter_try < max_jitter; jitter_try *= 10) { - B.diagonal().array() += jitter_try; + // Remove previously added jitter before adding the new (larger) amount, + // so that the total diagonal perturbation is exactly jitter_try. + B.diagonal().array() += (jitter_try - prev_jitter); + prev_jitter = jitter_try; llt_B.compute(B); if (llt_B.info() == Eigen::Success) { break; @@ -935,6 +964,9 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, scratch.alpha() = 1.0; update_fun(scratch, state.curr(), state.prev(), scratch.eval_, state.wolfe_info.p_); + // Save the full Newton step objective before the Wolfe line search + // overwrites scratch with intermediate trial points. + const double full_newton_obj = scratch.eval_.obj(); if (scratch.alpha() <= options.line_search.min_alpha) { state.wolfe_status.accept_ = false; finish_update = true; @@ -953,15 +985,42 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, state.wolfe_status = internal::wolfe_line_search( state.wolfe_info, update_fun, options.line_search, msgs); } + // When the Wolfe line search rejects, don't immediately terminate. + // Instead, let the Newton loop try at least one more iteration. + // The original code compared the stale curr.obj() (which equalled + // prev.obj() after the swap in update_next_step) and would always + // terminate on ANY Wolfe rejection — even on the very first Newton + // step. Now we only declare search_failed if the full Newton step + // itself didn't improve the objective. + bool search_failed; + if (!state.wolfe_status.accept_) { + if (full_newton_obj > state.prev().obj()) { + // The full Newton step (evaluated before Wolfe ran) improved + // the objective. Re-evaluate scratch at the full Newton step + // so we can accept it as the current iterate. + scratch.eval_.alpha() = 1.0; + update_fun(scratch, state.curr(), state.prev(), scratch.eval_, + state.wolfe_info.p_); + state.curr().update(scratch); + state.wolfe_status.accept_ = true; + search_failed = false; + } else { + search_failed = true; + } + } else { + search_failed = false; + } /** - * Stop when objective change is small, or when a rejected Wolfe step - * fails to improve; finish_update then exits the Newton loop. + * Stop when objective change is small (absolute AND relative), or when + * a rejected Wolfe step fails to improve; finish_update then exits the + * Newton loop. */ + double obj_change = std::abs(state.curr().obj() - state.prev().obj()); bool objective_converged - = std::abs(state.curr().obj() - state.prev().obj()) - < options.tolerance; - bool search_failed = (!state.wolfe_status.accept_ - && state.curr().obj() <= state.prev().obj()); + = obj_change < options.tolerance + && obj_change + < options.tolerance + * std::abs(state.prev().obj()); finish_update = objective_converged || search_failed; } if (finish_update) { @@ -1152,7 +1211,23 @@ inline auto laplace_marginal_density_est( return laplace_likelihood::theta_grad(ll_fun, theta_val, ll_args, msgs); }; decltype(auto) theta_init = theta_init_impl(theta_size, options); - internal::NewtonState state(theta_size, obj_fun, theta_grad_f, theta_init); + // When the user supplies a non-zero theta_init, we must initialise a + // consistently so that the invariant theta = Sigma * a holds. Otherwise + // the prior term -0.5 * a'*theta vanishes (a=0 while theta!=0), inflating + // the initial objective and causing the Wolfe line search to reject the + // first Newton step. + auto make_state = [&](auto&& theta_0) { + if constexpr (InitTheta) { + Eigen::VectorXd a_init = covariance.llt().solve( + Eigen::VectorXd(theta_0)); + return internal::NewtonState(theta_size, obj_fun, theta_grad_f, + a_init, theta_0); + } else { + return internal::NewtonState(theta_size, obj_fun, theta_grad_f, + theta_0); + } + }; + auto state = make_state(theta_init); // Start with safe step size auto update_fun = create_update_fun( std::move(obj_fun), std::move(theta_grad_f), covariance, options); diff --git a/stan/math/mix/functor/wolfe_line_search.hpp b/stan/math/mix/functor/wolfe_line_search.hpp index 15ac0abe9e0..0b9adad6150 100644 --- a/stan/math/mix/functor/wolfe_line_search.hpp +++ b/stan/math/mix/functor/wolfe_line_search.hpp @@ -512,6 +512,29 @@ struct WolfeInfo { "theta and likelihood arguments."); } } + /** + * Construct WolfeInfo with a consistent (a_init, theta_init) pair. + * + * When the caller supplies a non-zero theta_init, the corresponding + * a_init = Sigma^{-1} * theta_init must be provided so that the + * invariant theta = Sigma * a holds at initialization. This avoids + * an inflated initial objective (the prior term -0.5 * a'*theta would + * otherwise vanish when a is zero but theta is not). + */ + template + WolfeInfo(ObjFun&& obj_fun, const Eigen::VectorXd& a_init, Theta0&& theta0, + ThetaGradF&& theta_grad_f, int /*tag*/) + : curr_(std::forward(obj_fun), a_init, + std::forward(theta0), + std::forward(theta_grad_f)), + prev_(curr_), + scratch_(a_init.size()) { + if (!std::isfinite(curr_.obj())) { + throw std::domain_error( + "laplace_marginal_density: log likelihood is not finite at initial " + "theta and likelihood arguments."); + } + } WolfeInfo(WolfeData&& curr, WolfeData&& prev) : curr_(std::move(curr)), prev_(std::move(prev)), @@ -902,9 +925,10 @@ inline WolfeStatus wolfe_line_search(Info& wolfe_info, UpdateFun&& update_fun, } else { // [3] high = mid; } + } else { + // [4] + high = mid; } - // [4] - high = mid; } else { // [5] high = mid; From 8dcfbafc12a90ce3782fccbe19690d0813ec5605 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 19 Mar 2026 11:57:57 -0400 Subject: [PATCH 2/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../laplace_marginal_density_estimator.hpp | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index f57ce942594..fae7ba938ce 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -394,8 +394,7 @@ struct NewtonState { */ template NewtonState(int theta_size, ObjFun&& obj_fun, ThetaGradFun&& theta_grad_f, - const Eigen::VectorXd& a_init, - ThetaInitializer&& theta_init) + const Eigen::VectorXd& a_init, ThetaInitializer&& theta_init) : wolfe_info(std::forward(obj_fun), a_init, std::forward(theta_init), std::forward(theta_grad_f), 0), @@ -1018,9 +1017,7 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, double obj_change = std::abs(state.curr().obj() - state.prev().obj()); bool objective_converged = obj_change < options.tolerance - && obj_change - < options.tolerance - * std::abs(state.prev().obj()); + && obj_change < options.tolerance * std::abs(state.prev().obj()); finish_update = objective_converged || search_failed; } if (finish_update) { @@ -1218,13 +1215,11 @@ inline auto laplace_marginal_density_est( // first Newton step. auto make_state = [&](auto&& theta_0) { if constexpr (InitTheta) { - Eigen::VectorXd a_init = covariance.llt().solve( - Eigen::VectorXd(theta_0)); - return internal::NewtonState(theta_size, obj_fun, theta_grad_f, - a_init, theta_0); - } else { - return internal::NewtonState(theta_size, obj_fun, theta_grad_f, + Eigen::VectorXd a_init = covariance.llt().solve(Eigen::VectorXd(theta_0)); + return internal::NewtonState(theta_size, obj_fun, theta_grad_f, a_init, theta_0); + } else { + return internal::NewtonState(theta_size, obj_fun, theta_grad_f, theta_0); } }; auto state = make_state(theta_init); From c99fbcf41208f3463314c060f2420ff104345eab Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Mon, 23 Mar 2026 14:49:56 -0400 Subject: [PATCH 3/7] Fix NewtonState and WolfeInfo constructors --- .../laplace_marginal_density_estimator.hpp | 40 +++---------------- stan/math/mix/functor/wolfe_line_search.hpp | 22 +++------- 2 files changed, 10 insertions(+), 52 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index fae7ba938ce..30a1605a55c 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -359,27 +359,6 @@ struct NewtonState { */ bool final_loop = false; - /** - * @brief Constructs Newton state with given dimensions and functors. - * - * @tparam ThetaInitializer Type of the initial theta provider - * @param theta_size Dimension of the latent space - * @param obj_fun Objective function: (a, theta) -> double - * @param theta_grad_f Gradient function: theta -> grad - * @param theta_init Initial theta value or provider - */ - template - NewtonState(int theta_size, ObjFun&& obj_fun, ThetaGradFun&& theta_grad_f, - ThetaInitializer&& theta_init) - : wolfe_info(std::forward(obj_fun), theta_size, - std::forward(theta_init), - std::forward(theta_grad_f)), - b(theta_size), - B(theta_size, theta_size), - prev_g(theta_size) { - wolfe_status.num_backtracks_ = -1; // Safe initial value for BB step - } - /** * @brief Constructs Newton state with a consistent (a_init, theta_init) pair. * @@ -392,12 +371,12 @@ struct NewtonState { * @param a_init Initial a value consistent with theta_init * @param theta_init Initial theta value */ - template + template NewtonState(int theta_size, ObjFun&& obj_fun, ThetaGradFun&& theta_grad_f, - const Eigen::VectorXd& a_init, ThetaInitializer&& theta_init) - : wolfe_info(std::forward(obj_fun), a_init, + CovarianceT&& covariance, ThetaInitializer&& theta_init) + : wolfe_info(std::forward(obj_fun), covariance.llt().solve(theta_init), std::forward(theta_init), - std::forward(theta_grad_f), 0), + std::forward(theta_grad_f)), b(theta_size), B(theta_size, theta_size), prev_g(theta_size) { @@ -1213,16 +1192,7 @@ inline auto laplace_marginal_density_est( // the prior term -0.5 * a'*theta vanishes (a=0 while theta!=0), inflating // the initial objective and causing the Wolfe line search to reject the // first Newton step. - auto make_state = [&](auto&& theta_0) { - if constexpr (InitTheta) { - Eigen::VectorXd a_init = covariance.llt().solve(Eigen::VectorXd(theta_0)); - return internal::NewtonState(theta_size, obj_fun, theta_grad_f, a_init, - theta_0); - } else { - return internal::NewtonState(theta_size, obj_fun, theta_grad_f, theta_0); - } - }; - auto state = make_state(theta_init); + auto state = NewtonState(theta_size, obj_fun, theta_grad_f, covariance, theta_init); // Start with safe step size auto update_fun = create_update_fun( std::move(obj_fun), std::move(theta_grad_f), covariance, options); diff --git a/stan/math/mix/functor/wolfe_line_search.hpp b/stan/math/mix/functor/wolfe_line_search.hpp index 0b9adad6150..e23f7498f44 100644 --- a/stan/math/mix/functor/wolfe_line_search.hpp +++ b/stan/math/mix/functor/wolfe_line_search.hpp @@ -499,19 +499,7 @@ struct WolfeInfo { Eigen::VectorXd p_; // Initial directional derivative double init_dir_; - template - WolfeInfo(ObjFun&& obj_fun, Eigen::Index n, Theta0&& theta0, - ThetaGradF&& theta_grad_f) - : curr_(std::forward(obj_fun), n, std::forward(theta0), - std::forward(theta_grad_f)), - prev_(curr_), - scratch_(n) { - if (!std::isfinite(curr_.obj())) { - throw std::domain_error( - "laplace_marginal_density: log likelihood is not finite at initial " - "theta and likelihood arguments."); - } - } + /** * Construct WolfeInfo with a consistent (a_init, theta_init) pair. * @@ -521,10 +509,10 @@ struct WolfeInfo { * an inflated initial objective (the prior term -0.5 * a'*theta would * otherwise vanish when a is zero but theta is not). */ - template - WolfeInfo(ObjFun&& obj_fun, const Eigen::VectorXd& a_init, Theta0&& theta0, - ThetaGradF&& theta_grad_f, int /*tag*/) - : curr_(std::forward(obj_fun), a_init, + template + WolfeInfo(ObjFun&& obj_fun, AInit&& a_init, Theta0&& theta0, + ThetaGradF&& theta_grad_f) + : curr_(std::forward(obj_fun), std::forward(a_init), std::forward(theta0), std::forward(theta_grad_f)), prev_(curr_), From 184f792764e6948a10c53ad8521b2193fcc36b35 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 24 Mar 2026 10:15:23 -0400 Subject: [PATCH 4/7] reuse initial full proposal step in laplace --- .../laplace_marginal_density_estimator.hpp | 113 +++++++-------- stan/math/mix/functor/wolfe_line_search.hpp | 6 + ...aplace_marginal_density_estimator_test.cpp | 133 ++++++++++++++++++ .../math/laplace/wolfe_line_search_test.cpp | 58 ++++++++ 4 files changed, 251 insertions(+), 59 deletions(-) create mode 100644 test/unit/math/laplace/laplace_marginal_density_estimator_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 30a1605a55c..b304ac94840 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -342,6 +342,9 @@ struct NewtonState { /** @brief Status of the most recent Wolfe line search */ WolfeStatus wolfe_status; + /** @brief Cached proposal evaluated before the Wolfe line search. */ + WolfeData proposal; + /** @brief Workspace vector: b = W * theta + grad(log_lik) */ Eigen::VectorXd b; @@ -377,6 +380,7 @@ struct NewtonState { : wolfe_info(std::forward(obj_fun), covariance.llt().solve(theta_init), std::forward(theta_init), std::forward(theta_grad_f)), + proposal(theta_size), b(theta_size), B(theta_size, theta_size), prev_g(theta_size) { @@ -407,9 +411,12 @@ struct NewtonState { */ const auto& prev() const& { return wolfe_info.prev_; } auto&& prev() && { return std::move(wolfe_info).prev(); } + auto& proposal_step() & { return proposal; } + const auto& proposal_step() const& { return proposal; } + auto&& proposal_step() && { return std::move(proposal); } template inline void update_next_step(const Options& options) { - this->prev().update(this->curr()); + this->prev().swap(this->curr()); this->curr().alpha() = std::clamp(this->curr().alpha(), 0.0, options.line_search.max_alpha); } @@ -485,7 +492,8 @@ struct CholeskyWSolverDiag { * @tparam LLFun Type of the log-likelihood functor * @tparam LLTupleArgs Type of the likelihood arguments tuple * @tparam CovarMat Type of the covariance matrix - * @param[in,out] state Shared Newton state (modified: B, b, curr().a()) + * @param[in,out] state Shared Newton state (modified: B, b, + * proposal_step().a()) * @param[in] ll_fun Log-likelihood functor * @param[in,out] ll_args Additional arguments for the likelihood * @param[in] covariance Prior covariance matrix Sigma @@ -521,12 +529,12 @@ struct CholeskyWSolverDiag { // 3. Factorize B with jittering fallback llt_with_jitter(llt_B, state.B); - // 4. Solve for curr.a + // 4. Solve for the raw Newton proposal in a-space. state.b.noalias() = (W_diag.array() * state.prev().theta().array()).matrix() + state.prev().theta_grad(); auto L = llt_B.matrixL(); auto LT = llt_B.matrixU(); - state.curr().a().noalias() + state.proposal_step().a().noalias() = state.b - W_r_diag.asDiagonal() * LT.solve( @@ -615,7 +623,8 @@ struct CholeskyWSolverBlock { * @tparam LLFun Type of the log-likelihood functor * @tparam LLTupleArgs Type of the likelihood arguments tuple * @tparam CovarMat Type of the covariance matrix - * @param[in,out] state Shared Newton state (modified: B, b, curr().a()) + * @param[in,out] state Shared Newton state (modified: B, b, + * proposal_step().a()) * @param[in] ll_fun Log-likelihood functor * @param[in,out] ll_args Additional arguments for the likelihood * @param[in] covariance Prior covariance matrix Sigma @@ -653,12 +662,12 @@ struct CholeskyWSolverBlock { // 4. Factorize B with jittering fallback llt_with_jitter(llt_B, state.B); - // 5. Solve for curr.a + // 5. Solve for the raw Newton proposal in a-space. state.b.noalias() = W_block * state.prev().theta() + state.prev().theta_grad(); auto L = llt_B.matrixL(); auto LT = llt_B.matrixU(); - state.curr().a().noalias() + state.proposal_step().a().noalias() = state.b - W_r * LT.solve(L.solve(W_r * (covariance * state.b))); } @@ -736,7 +745,7 @@ struct CholeskyKSolver { * @tparam LLFun Type of the log-likelihood functor * @tparam LLTupleArgs Type of the likelihood arguments tuple * @tparam CovarMat Type of the covariance matrix - * @param[in] state Shared Newton state (modified: B, b, curr().a()) + * @param[in] state Shared Newton state (modified: B, b, proposal_step().a()) * @param[in] ll_fun Log-likelihood functor * @param[in] ll_args Additional arguments for the likelihood * @param[in] covariance Prior covariance matrix Sigma @@ -763,12 +772,12 @@ struct CholeskyKSolver { // 3. Factorize B with jittering fallback llt_with_jitter(llt_B, state.B); - // 4. Solve for curr.a + // 4. Solve for the raw Newton proposal in a-space. state.b.noalias() = W_full * state.prev().theta() + state.prev().theta_grad(); auto L = llt_B.matrixL(); auto LT = llt_B.matrixU(); - state.curr().a().noalias() + state.proposal_step().a().noalias() = K_root.transpose().template triangularView().solve( LT.solve(L.solve(K_root.transpose() * state.b))); } @@ -833,7 +842,7 @@ struct LUSolver { * @tparam LLFun Type of the log-likelihood functor * @tparam LLTupleArgs Type of the likelihood arguments tuple * @tparam CovarMat Type of the covariance matrix - * @param[in,out] state Shared Newton state (modified: b, curr().a()) + * @param[in,out] state Shared Newton state (modified: b, proposal_step().a()) * @param[in] ll_fun Log-likelihood functor * @param[in,out] ll_args Additional arguments for the likelihood * @param[in] covariance Prior covariance matrix Sigma @@ -855,10 +864,10 @@ struct LUSolver { lu.compute(Eigen::MatrixXd::Identity(theta_size, theta_size) + covariance * W_full); - // 3. Solve for curr.a + // 3. Solve for the raw Newton proposal in a-space. state.b.noalias() = W_full * state.prev().theta() + state.prev().theta_grad(); - state.curr().a().noalias() + state.proposal_step().a().noalias() = state.b - W_full * lu.solve(covariance * state.b); } @@ -932,29 +941,32 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, solver.solve_step(state, ll_fun, ll_args, covariance, options.hessian_block_size, msgs); if (!state.final_loop) { - state.wolfe_info.p_ = state.curr().a() - state.prev().a(); + auto&& proposal = state.proposal_step(); + state.wolfe_info.p_ = proposal.a() - state.prev().a(); state.prev_g.noalias() = -covariance * state.prev().a() + covariance * state.prev().theta_grad(); state.wolfe_info.init_dir_ = state.prev_g.dot(state.wolfe_info.p_); // Flip direction if not ascending state.wolfe_info.flip_direction(); auto&& scratch = state.wolfe_info.scratch_; - scratch.alpha() = 1.0; - update_fun(scratch, state.curr(), state.prev(), scratch.eval_, - state.wolfe_info.p_); - // Save the full Newton step objective before the Wolfe line search - // overwrites scratch with intermediate trial points. - const double full_newton_obj = scratch.eval_.obj(); - if (scratch.alpha() <= options.line_search.min_alpha) { - state.wolfe_status.accept_ = false; - finish_update = true; + proposal.eval_.alpha() = 1.0; + const bool proposal_valid = update_fun( + proposal, state.curr(), state.prev(), proposal.eval_, + state.wolfe_info.p_); + const bool cached_proposal_ok + = proposal_valid && std::isfinite(proposal.obj()) + && std::isfinite(proposal.dir()) + && proposal.alpha() > options.line_search.min_alpha; + if (!cached_proposal_ok) { + state.wolfe_status + = WolfeStatus{WolfeReturn::StepTooSmall, 1, 0, false}; } else if (options.line_search.max_iterations == 0) { - state.curr().update(scratch); - state.wolfe_status.accept_ = true; + state.curr().update(proposal); + state.wolfe_status = WolfeStatus{WolfeReturn::Continue, 1, 0, true}; } else { - Eigen::VectorXd s = scratch.a() - state.prev().a(); + Eigen::VectorXd s = proposal.a() - state.prev().a(); auto full_step_grad - = (-covariance * scratch.a() + covariance * scratch.theta_grad()) + = (-covariance * proposal.a() + covariance * proposal.theta_grad()) .eval(); state.curr().alpha() = barzilai_borwein_step_size( s, full_step_grad, state.prev_g, state.prev().alpha(), @@ -963,47 +975,30 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, state.wolfe_status = internal::wolfe_line_search( state.wolfe_info, update_fun, options.line_search, msgs); } - // When the Wolfe line search rejects, don't immediately terminate. - // Instead, let the Newton loop try at least one more iteration. - // The original code compared the stale curr.obj() (which equalled - // prev.obj() after the swap in update_next_step) and would always - // terminate on ANY Wolfe rejection — even on the very first Newton - // step. Now we only declare search_failed if the full Newton step - // itself didn't improve the objective. - bool search_failed; - if (!state.wolfe_status.accept_) { - if (full_newton_obj > state.prev().obj()) { - // The full Newton step (evaluated before Wolfe ran) improved - // the objective. Re-evaluate scratch at the full Newton step - // so we can accept it as the current iterate. - scratch.eval_.alpha() = 1.0; - update_fun(scratch, state.curr(), state.prev(), scratch.eval_, - state.wolfe_info.p_); - state.curr().update(scratch); - state.wolfe_status.accept_ = true; - search_failed = false; - } else { - search_failed = true; - } - } else { + bool search_failed = !state.wolfe_status.accept_; + const bool proposal_armijo_ok + = cached_proposal_ok + && internal::check_armijo( + proposal.obj(), state.prev().obj(), proposal.alpha(), + state.wolfe_info.init_dir_, options.line_search); + if (search_failed && proposal_armijo_ok) { + state.curr().update(proposal); + state.wolfe_status = WolfeStatus{WolfeReturn::Armijo, + state.wolfe_status.num_evals_, + state.wolfe_status.num_backtracks_, + true}; search_failed = false; } - /** - * Stop when objective change is small (absolute AND relative), or when - * a rejected Wolfe step fails to improve; finish_update then exits the - * Newton loop. - */ - double obj_change = std::abs(state.curr().obj() - state.prev().obj()); bool objective_converged - = obj_change < options.tolerance - && obj_change < options.tolerance * std::abs(state.prev().obj()); + = state.wolfe_status.accept_ + && std::abs(state.curr().obj() - state.prev().obj()) + < options.tolerance; finish_update = objective_converged || search_failed; } if (finish_update) { if (!state.final_loop && state.wolfe_status.accept_) { // Do one final loop with exact wolfe conditions state.final_loop = true; - // NOTE: Swapping here so we need to swap prev and curr later state.update_next_step(options); continue; } diff --git a/stan/math/mix/functor/wolfe_line_search.hpp b/stan/math/mix/functor/wolfe_line_search.hpp index e23f7498f44..3a27d1ad7ec 100644 --- a/stan/math/mix/functor/wolfe_line_search.hpp +++ b/stan/math/mix/functor/wolfe_line_search.hpp @@ -461,6 +461,12 @@ struct WolfeData { a_.swap(other.a_); eval_ = other.eval_; } + void swap(WolfeData& other) { + theta_.swap(other.theta_); + theta_grad_.swap(other.theta_grad_); + a_.swap(other.a_); + std::swap(eval_, other.eval_); + } void update(WolfeData& other, const Eval& eval) { theta_.swap(other.theta_); a_.swap(other.a_); diff --git a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp new file mode 100644 index 00000000000..e2ab435a49b --- /dev/null +++ b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp @@ -0,0 +1,133 @@ +#include +#include +#include + +#include +#include +#include + +namespace stan::math { +namespace { + +struct IdentityCovariance { + template + Eigen::MatrixXd operator()(Stream* /*msgs*/) const { + return Eigen::MatrixXd::Identity(1, 1); + } +}; + +struct QuarticLikelihood { + template + auto operator()(const Theta& theta, std::ostream* /*msgs*/) const { + const auto& x = theta(0); + const auto x_sq = stan::math::square(x); + return 2.0 * x - 0.5 * x_sq - 0.5 * stan::math::square(x_sq); + } +}; + +struct TinyQuarticLikelihood { + template + auto operator()(const Theta& theta, std::ostream* /*msgs*/) const { + return 1e-8 * QuarticLikelihood{}(theta, nullptr); + } +}; + +struct StubNewtonSolver { + double proposal_a; + + template + void solve_step(NewtonStateT& state, const LLFun& /*ll_fun*/, + const LLTupleArgs& /*ll_args*/, + const CovarMat& /*covariance*/, int /*hessian_block_size*/, + std::ostream* /*msgs*/) const { + state.proposal_step().a()(0) = proposal_a; + } + + double compute_log_determinant() const { return 0.0; } + + template + double build_result(NewtonStateT& state, double /*log_det*/) const { + return state.prev().a()(0); + } +}; + +template +double run_laplace(const Likelihood& ll_fun, double theta0_value, + double tolerance, int max_num_steps, + int max_steps_line_search, std::ostream* msgs) { + Eigen::VectorXd theta0(1); + theta0 << theta0_value; + return stan::math::laplace_marginal_tol( + ll_fun, std::tuple<>{}, 1, IdentityCovariance{}, std::tuple<>{}, + std::make_tuple(theta0, tolerance, max_num_steps, 1, + max_steps_line_search, true), + msgs); +} + +TEST(LaplaceMarginalDensityEstimator, PublicLineSearchMatchesDirectStep) { + std::ostringstream no_search_msgs; + std::ostringstream wolfe_msgs; + + const double no_search = run_laplace(QuarticLikelihood{}, 2.0, 1e-12, 50, 0, + &no_search_msgs); + const double with_wolfe = run_laplace(QuarticLikelihood{}, 2.0, 1e-12, 50, + 1000, &wolfe_msgs); + + EXPECT_TRUE(std::isfinite(no_search)); + EXPECT_TRUE(std::isfinite(with_wolfe)); + EXPECT_NEAR(no_search, with_wolfe, 1e-8); +} + +TEST(LaplaceMarginalDensityEstimator, AbsoluteObjectiveToleranceStopsNearZero) { + std::ostringstream msgs; + + const double result + = run_laplace(TinyQuarticLikelihood{}, 0.0, 1e-8, 6, 1000, &msgs); + + EXPECT_TRUE(std::isfinite(result)); + EXPECT_EQ(msgs.str().find("max number of iterations"), std::string::npos); +} + +TEST(LaplaceMarginalDensityEstimator, + InvalidCachedProposalDoesNotTriggerArmijoFallback) { + Eigen::MatrixXd covariance = Eigen::MatrixXd::Identity(1, 1); + Eigen::VectorXd theta0 = Eigen::VectorXd::Zero(1); + auto obj_fun = [](const auto& /*a*/, const auto& /*theta*/) { + return -1.0; + }; + auto theta_grad_f = [](const auto& theta) { + return Eigen::VectorXd::Zero(theta.size()); + }; + internal::NewtonState state(1, obj_fun, theta_grad_f, covariance, theta0); + laplace_options_base options; + options.hessian_block_size = 1; + options.max_num_steps = 1; + options.tolerance = 1e-12; + options.line_search.max_iterations = 5; + options.line_search.min_alpha = 1e-8; + + StubNewtonSolver solver{5.0}; + Eigen::Index step_iter = 1; + auto failing_update = [min_alpha = options.line_search.min_alpha]( + auto& /*proposal*/, auto&& /*curr*/, + auto&& /*prev*/, auto& eval_in, auto&& /*p*/) { + eval_in.alpha() = 0.5 * min_alpha; + return false; + }; + auto unused_ll = [](const auto& /*theta*/, std::ostream* /*msgs*/) { + return 0.0; + }; + + const double result + = internal::run_newton_loop(solver, state, options, step_iter, unused_ll, + std::tuple<>{}, covariance, failing_update, + nullptr); + + EXPECT_DOUBLE_EQ(result, 0.0); + EXPECT_FALSE(state.wolfe_status.accept_); + EXPECT_EQ(state.wolfe_status.stop_, internal::WolfeReturn::StepTooSmall); +} + +} // namespace +} // namespace stan::math diff --git a/test/unit/math/laplace/wolfe_line_search_test.cpp b/test/unit/math/laplace/wolfe_line_search_test.cpp index f3b64917652..a5a1691c4ab 100644 --- a/test/unit/math/laplace/wolfe_line_search_test.cpp +++ b/test/unit/math/laplace/wolfe_line_search_test.cpp @@ -781,6 +781,64 @@ TEST(WolfeLineSearch, HonorsMaxAlphaBound) { EXPECT_NE(status.stop_, WolfeReturn::Fail); } +// Checks that zoom case 2 updates the low endpoint, not the high endpoint. +TEST(WolfeLineSearch, ZoomPreservesCaseTwoLowUpdate) { + using internal::Eval; + + WolfeInfo info(1); + info.prev_.a()(0) = 0.0; + info.prev_.theta()(0) = 0.0; + info.prev_.theta_grad()(0) = 0.0; + info.prev_.obj() = 0.0; + info.prev_.alpha() = 0.0; + info.curr_.a()(0) = 1.0; + info.curr_.theta()(0) = 1.0; + info.curr_.theta_grad()(0) = 0.0; + info.curr_.obj() = 0.0; + info.curr_.alpha() = 0.25; // alpha_start = 0.5 + info.p_(0) = 1.0; + info.init_dir_ = 4.0; + + laplace_line_search_options opt{}; + opt.c1 = 0.5; + opt.c2 = 0.1; + opt.scale_up = 2.0; + opt.max_alpha = 1.0; + + auto scripted_update = [](auto& proposal, auto&&, auto&&, Eval& eval, + auto&&) { + proposal.a()(0) = eval.alpha(); + proposal.theta()(0) = eval.alpha(); + proposal.theta_grad()(0) = 0.0; + const double alpha = eval.alpha(); + if (std::abs(alpha - 0.5) < 1e-12) { + eval.obj() = 0.5; + eval.dir() = -2.0; + } else if (std::abs(alpha - 0.25) < 1e-12) { + eval.obj() = 0.75; + eval.dir() = 1.0; + } else if (std::abs(alpha - 0.375) < 1e-12) { + eval.obj() = 0.95; + eval.dir() = 0.0; + } else if (std::abs(alpha - 0.125) < 1e-12) { + eval.obj() = 0.6; + eval.dir() = 0.0; + } else { + ADD_FAILURE() << "Unexpected alpha " << alpha; + eval.obj() = -1.0; + eval.dir() = -2.0; + } + }; + + auto status = wolfe_line_search(info, scripted_update, opt, + static_cast(nullptr)); + + EXPECT_TRUE(status.accept_); + EXPECT_EQ(status.stop_, WolfeReturn::Wolfe); + EXPECT_DOUBLE_EQ(info.curr_.alpha(), 0.375); + EXPECT_DOUBLE_EQ(info.curr_.a()(0), 0.375); +} + // Checks that the cubic-or-bisect chooser returns an interior maximiser. TEST(CubicOrBisect, ReturnsInteriorMaximiser) { using stan::math::internal::cubic_spline; From 52aaa067bf967484109858aa60521db6418c8ef8 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 24 Mar 2026 10:16:20 -0400 Subject: [PATCH 5/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../laplace_marginal_density_estimator.hpp | 26 +++++----- stan/math/mix/functor/wolfe_line_search.hpp | 3 +- ...aplace_marginal_density_estimator_test.cpp | 29 +++++------ .../math/laplace/wolfe_line_search_test.cpp | 48 +++++++++---------- 4 files changed, 52 insertions(+), 54 deletions(-) diff --git a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp index b304ac94840..2dbbae73cbc 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -374,10 +374,12 @@ struct NewtonState { * @param a_init Initial a value consistent with theta_init * @param theta_init Initial theta value */ - template + template NewtonState(int theta_size, ObjFun&& obj_fun, ThetaGradFun&& theta_grad_f, CovarianceT&& covariance, ThetaInitializer&& theta_init) - : wolfe_info(std::forward(obj_fun), covariance.llt().solve(theta_init), + : wolfe_info(std::forward(obj_fun), + covariance.llt().solve(theta_init), std::forward(theta_init), std::forward(theta_grad_f)), proposal(theta_size), @@ -950,9 +952,9 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, state.wolfe_info.flip_direction(); auto&& scratch = state.wolfe_info.scratch_; proposal.eval_.alpha() = 1.0; - const bool proposal_valid = update_fun( - proposal, state.curr(), state.prev(), proposal.eval_, - state.wolfe_info.p_); + const bool proposal_valid + = update_fun(proposal, state.curr(), state.prev(), proposal.eval_, + state.wolfe_info.p_); const bool cached_proposal_ok = proposal_valid && std::isfinite(proposal.obj()) && std::isfinite(proposal.dir()) @@ -979,14 +981,13 @@ inline auto run_newton_loop(SolverPolicy& solver, NewtonStateT& state, const bool proposal_armijo_ok = cached_proposal_ok && internal::check_armijo( - proposal.obj(), state.prev().obj(), proposal.alpha(), - state.wolfe_info.init_dir_, options.line_search); + proposal.obj(), state.prev().obj(), proposal.alpha(), + state.wolfe_info.init_dir_, options.line_search); if (search_failed && proposal_armijo_ok) { state.curr().update(proposal); - state.wolfe_status = WolfeStatus{WolfeReturn::Armijo, - state.wolfe_status.num_evals_, - state.wolfe_status.num_backtracks_, - true}; + state.wolfe_status + = WolfeStatus{WolfeReturn::Armijo, state.wolfe_status.num_evals_, + state.wolfe_status.num_backtracks_, true}; search_failed = false; } bool objective_converged @@ -1187,7 +1188,8 @@ inline auto laplace_marginal_density_est( // the prior term -0.5 * a'*theta vanishes (a=0 while theta!=0), inflating // the initial objective and causing the Wolfe line search to reject the // first Newton step. - auto state = NewtonState(theta_size, obj_fun, theta_grad_f, covariance, theta_init); + auto state + = NewtonState(theta_size, obj_fun, theta_grad_f, covariance, theta_init); // Start with safe step size auto update_fun = create_update_fun( std::move(obj_fun), std::move(theta_grad_f), covariance, options); diff --git a/stan/math/mix/functor/wolfe_line_search.hpp b/stan/math/mix/functor/wolfe_line_search.hpp index 3a27d1ad7ec..c52960ad499 100644 --- a/stan/math/mix/functor/wolfe_line_search.hpp +++ b/stan/math/mix/functor/wolfe_line_search.hpp @@ -515,7 +515,8 @@ struct WolfeInfo { * an inflated initial objective (the prior term -0.5 * a'*theta would * otherwise vanish when a is zero but theta is not). */ - template + template WolfeInfo(ObjFun&& obj_fun, AInit&& a_init, Theta0&& theta0, ThetaGradF&& theta_grad_f) : curr_(std::forward(obj_fun), std::forward(a_init), diff --git a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp index e2ab435a49b..67995dd5bcc 100644 --- a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp +++ b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp @@ -69,10 +69,10 @@ TEST(LaplaceMarginalDensityEstimator, PublicLineSearchMatchesDirectStep) { std::ostringstream no_search_msgs; std::ostringstream wolfe_msgs; - const double no_search = run_laplace(QuarticLikelihood{}, 2.0, 1e-12, 50, 0, - &no_search_msgs); - const double with_wolfe = run_laplace(QuarticLikelihood{}, 2.0, 1e-12, 50, - 1000, &wolfe_msgs); + const double no_search + = run_laplace(QuarticLikelihood{}, 2.0, 1e-12, 50, 0, &no_search_msgs); + const double with_wolfe + = run_laplace(QuarticLikelihood{}, 2.0, 1e-12, 50, 1000, &wolfe_msgs); EXPECT_TRUE(std::isfinite(no_search)); EXPECT_TRUE(std::isfinite(with_wolfe)); @@ -93,12 +93,9 @@ TEST(LaplaceMarginalDensityEstimator, InvalidCachedProposalDoesNotTriggerArmijoFallback) { Eigen::MatrixXd covariance = Eigen::MatrixXd::Identity(1, 1); Eigen::VectorXd theta0 = Eigen::VectorXd::Zero(1); - auto obj_fun = [](const auto& /*a*/, const auto& /*theta*/) { - return -1.0; - }; - auto theta_grad_f = [](const auto& theta) { - return Eigen::VectorXd::Zero(theta.size()); - }; + auto obj_fun = [](const auto& /*a*/, const auto& /*theta*/) { return -1.0; }; + auto theta_grad_f + = [](const auto& theta) { return Eigen::VectorXd::Zero(theta.size()); }; internal::NewtonState state(1, obj_fun, theta_grad_f, covariance, theta0); laplace_options_base options; options.hessian_block_size = 1; @@ -115,14 +112,12 @@ TEST(LaplaceMarginalDensityEstimator, eval_in.alpha() = 0.5 * min_alpha; return false; }; - auto unused_ll = [](const auto& /*theta*/, std::ostream* /*msgs*/) { - return 0.0; - }; + auto unused_ll + = [](const auto& /*theta*/, std::ostream* /*msgs*/) { return 0.0; }; - const double result - = internal::run_newton_loop(solver, state, options, step_iter, unused_ll, - std::tuple<>{}, covariance, failing_update, - nullptr); + const double result = internal::run_newton_loop( + solver, state, options, step_iter, unused_ll, std::tuple<>{}, covariance, + failing_update, nullptr); EXPECT_DOUBLE_EQ(result, 0.0); EXPECT_FALSE(state.wolfe_status.accept_); diff --git a/test/unit/math/laplace/wolfe_line_search_test.cpp b/test/unit/math/laplace/wolfe_line_search_test.cpp index a5a1691c4ab..90dbee9b765 100644 --- a/test/unit/math/laplace/wolfe_line_search_test.cpp +++ b/test/unit/math/laplace/wolfe_line_search_test.cpp @@ -805,30 +805,30 @@ TEST(WolfeLineSearch, ZoomPreservesCaseTwoLowUpdate) { opt.scale_up = 2.0; opt.max_alpha = 1.0; - auto scripted_update = [](auto& proposal, auto&&, auto&&, Eval& eval, - auto&&) { - proposal.a()(0) = eval.alpha(); - proposal.theta()(0) = eval.alpha(); - proposal.theta_grad()(0) = 0.0; - const double alpha = eval.alpha(); - if (std::abs(alpha - 0.5) < 1e-12) { - eval.obj() = 0.5; - eval.dir() = -2.0; - } else if (std::abs(alpha - 0.25) < 1e-12) { - eval.obj() = 0.75; - eval.dir() = 1.0; - } else if (std::abs(alpha - 0.375) < 1e-12) { - eval.obj() = 0.95; - eval.dir() = 0.0; - } else if (std::abs(alpha - 0.125) < 1e-12) { - eval.obj() = 0.6; - eval.dir() = 0.0; - } else { - ADD_FAILURE() << "Unexpected alpha " << alpha; - eval.obj() = -1.0; - eval.dir() = -2.0; - } - }; + auto scripted_update + = [](auto& proposal, auto&&, auto&&, Eval& eval, auto&&) { + proposal.a()(0) = eval.alpha(); + proposal.theta()(0) = eval.alpha(); + proposal.theta_grad()(0) = 0.0; + const double alpha = eval.alpha(); + if (std::abs(alpha - 0.5) < 1e-12) { + eval.obj() = 0.5; + eval.dir() = -2.0; + } else if (std::abs(alpha - 0.25) < 1e-12) { + eval.obj() = 0.75; + eval.dir() = 1.0; + } else if (std::abs(alpha - 0.375) < 1e-12) { + eval.obj() = 0.95; + eval.dir() = 0.0; + } else if (std::abs(alpha - 0.125) < 1e-12) { + eval.obj() = 0.6; + eval.dir() = 0.0; + } else { + ADD_FAILURE() << "Unexpected alpha " << alpha; + eval.obj() = -1.0; + eval.dir() = -2.0; + } + }; auto status = wolfe_line_search(info, scripted_update, opt, static_cast(nullptr)); From 5425e6b47657167771d2e2cb1b80e0f7abd7523d Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 24 Mar 2026 12:05:41 -0400 Subject: [PATCH 6/7] fix doxygen and cpplint --- .../math/mix/functor/laplace_marginal_density_estimator.hpp | 6 +++++- .../laplace/laplace_marginal_density_estimator_test.cpp | 2 +- 2 files changed, 6 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 2dbbae73cbc..074838a3edb 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -367,10 +367,14 @@ struct NewtonState { * * When the caller supplies a non-zero theta_init, a_init = Sigma^{-1} * * theta_init must be provided to maintain the invariant theta = Sigma * a. - * + * @tparam ObjFun A callable type for the objective function + * @tparam ThetaGradFun A callable type for the theta gradient function + * @tparam CovarianceT A matrix type for the covariance (must support LLT solve) + * @tparam ThetaInitializer A type for the initial theta (e.g., Eigen vector) * @param theta_size Dimension of the latent space * @param obj_fun Objective function: (a, theta) -> double * @param theta_grad_f Gradient function: theta -> grad + * @param covariance Covariance matrix for the latent variables * @param a_init Initial a value consistent with theta_init * @param theta_init Initial theta value */ diff --git a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp index 67995dd5bcc..6f3241508cc 100644 --- a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp +++ b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp @@ -28,7 +28,7 @@ struct QuarticLikelihood { struct TinyQuarticLikelihood { template auto operator()(const Theta& theta, std::ostream* /*msgs*/) const { - return 1e-8 * QuarticLikelihood{}(theta, nullptr); + return 1e-8 * QuarticLikelihood{}(theta, nullptr); // NOLINT } }; From 520d494bafcc1f8dec6554f93da8162d6a0eac15 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 24 Mar 2026 12:06:50 -0400 Subject: [PATCH 7/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/mix/functor/laplace_marginal_density_estimator.hpp | 3 ++- .../math/laplace/laplace_marginal_density_estimator_test.cpp | 2 +- 2 files changed, 3 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 074838a3edb..25d419bf6c4 100644 --- a/stan/math/mix/functor/laplace_marginal_density_estimator.hpp +++ b/stan/math/mix/functor/laplace_marginal_density_estimator.hpp @@ -369,7 +369,8 @@ struct NewtonState { * theta_init must be provided to maintain the invariant theta = Sigma * a. * @tparam ObjFun A callable type for the objective function * @tparam ThetaGradFun A callable type for the theta gradient function - * @tparam CovarianceT A matrix type for the covariance (must support LLT solve) + * @tparam CovarianceT A matrix type for the covariance (must support LLT + * solve) * @tparam ThetaInitializer A type for the initial theta (e.g., Eigen vector) * @param theta_size Dimension of the latent space * @param obj_fun Objective function: (a, theta) -> double diff --git a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp index 6f3241508cc..7f82181b32d 100644 --- a/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp +++ b/test/unit/math/laplace/laplace_marginal_density_estimator_test.cpp @@ -28,7 +28,7 @@ struct QuarticLikelihood { struct TinyQuarticLikelihood { template auto operator()(const Theta& theta, std::ostream* /*msgs*/) const { - return 1e-8 * QuarticLikelihood{}(theta, nullptr); // NOLINT + return 1e-8 * QuarticLikelihood{}(theta, nullptr); // NOLINT } };