From adf582c49a94f7040ad208198cc3cc6f039f1ac8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Sat, 17 Oct 2026 16:07:13 +0200 Subject: [PATCH 1/2] Avoid allocations per evaluation and per iteration in the optimizers ModelAdaptor computes the gradient with a member std::vector instead of calling log_prob_grad, which allocated and copied the parameters and the gradient on every call. s and y of BFGSMinimizer and the two-loop workspace of LBFGSUpdate are members, and WolfeLineSearch no longer copies the previous gradient, a copy that reached the caller only as the gradient output of a failed search. --- src/stan/optimization/bfgs.hpp | 62 ++++++++++++----------- src/stan/optimization/bfgs_linesearch.hpp | 2 - src/stan/optimization/lbfgs_update.hpp | 5 +- 3 files changed, 37 insertions(+), 32 deletions(-) diff --git a/src/stan/optimization/bfgs.hpp b/src/stan/optimization/bfgs.hpp index 86a5a20ccae..6a5ff89a150 100644 --- a/src/stan/optimization/bfgs.hpp +++ b/src/stan/optimization/bfgs.hpp @@ -58,7 +58,7 @@ class BFGSMinimizer { protected: FunctorType _func; - VectorT _gk, _gk_1, _xk_1, _xk, _pk, _pk_1; + VectorT _gk, _gk_1, _xk_1, _xk, _pk, _pk_1, _sk, _yk; Scalar _fk, _fk_1, _alphak_1; Scalar _alpha, _alpha0; size_t _itNum; @@ -164,7 +164,6 @@ class BFGSMinimizer { int step() { Scalar gradNorm, stepNorm; - VectorT sk, yk; int retCode(0); int resetB(0); @@ -226,21 +225,21 @@ class BFGSMinimizer { _gk.swap(_gk_1); _pk.swap(_pk_1); - sk.noalias() = _xk - _xk_1; - yk.noalias() = _gk - _gk_1; + _sk.noalias() = _xk - _xk_1; + _yk.noalias() = _gk - _gk_1; gradNorm = _gk.norm(); - stepNorm = sk.norm(); + stepNorm = _sk.norm(); // Update QN approximation if (resetB) { // If the QN approximation was reset, automatically scale it // and update the step-size accordingly - Scalar B0fact = _qn.update(yk, sk, true); + Scalar B0fact = _qn.update(_yk, _sk, true); _pk_1 /= B0fact; _alphak_1 = _alpha * B0fact; } else { - _qn.update(yk, sk); + _qn.update(_yk, _sk); _alphak_1 = _alpha; } // Compute search direction for next step @@ -288,7 +287,8 @@ class ModelAdaptor { M &_model; std::vector _params_i; std::ostream *_msgs; - std::vector _x, _g; + std::vector _x; + std::vector _ad; size_t _fevals; public: @@ -327,36 +327,40 @@ class ModelAdaptor { } int operator()(const Eigen::Matrix &x, double &f, Eigen::Matrix &g) { - using Eigen::Dynamic; - using Eigen::Matrix; - using stan::math::index_type; - using stan::model::log_prob_grad; - typedef typename index_type >::type idx_t; - - _x.resize(x.size()); - for (idx_t i = 0; i < x.size(); i++) - _x[i] = x[i]; - _fevals++; + // stan::model::log_prob_grad on a member vector of vars: + // no allocation per call and no copies through std::vector + _ad.resize(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) { + _ad[i] = stan::math::var(x.coeff(i)); + } try { - f = -log_prob_grad(_model, _x, _params_i, _g, _msgs); + stan::math::var lp + = _model.template log_prob(_ad, _params_i, _msgs); + f = -lp.val(); + lp.grad(); } catch (const std::domain_error &e) { + stan::math::recover_memory(); if (_msgs) (*_msgs) << e.what() << std::endl; return 1; + } catch (...) { + stan::math::recover_memory(); + throw; } + g.resize(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) { + g.coeffRef(i) = -_ad[i].adj(); + } + stan::math::recover_memory(); - g.resize(_g.size()); - for (size_t i = 0; i < _g.size(); i++) { - if (!std::isfinite(_g[i])) { - if (_msgs) - *_msgs << "Error evaluating model log probability: " - "Non-finite gradient." - << std::endl; - return 3; - } - g[i] = -_g[i]; + if (!g.allFinite()) { + if (_msgs) + *_msgs << "Error evaluating model log probability: " + "Non-finite gradient." + << std::endl; + return 3; } if (std::isfinite(f)) { diff --git a/src/stan/optimization/bfgs_linesearch.hpp b/src/stan/optimization/bfgs_linesearch.hpp index a883845d050..d3a67adc289 100644 --- a/src/stan/optimization/bfgs_linesearch.hpp +++ b/src/stan/optimization/bfgs_linesearch.hpp @@ -241,7 +241,6 @@ int WolfeLineSearch(FunctorType &func, Scalar &alpha, XType &x1, Scalar alpha1(alpha); Scalar prevF(f0); - XType prevDF(gradx0); Scalar prevDFp(dfp); Scalar newDFp; @@ -287,7 +286,6 @@ int WolfeLineSearch(FunctorType &func, Scalar &alpha, XType &x1, alpha0 = alpha1; prevF = func_val; - std::swap(prevDF, gradx1); prevDFp = newDFp; alpha1 *= 10.0; diff --git a/src/stan/optimization/lbfgs_update.hpp b/src/stan/optimization/lbfgs_update.hpp index 303a852f164..896f1a7270f 100644 --- a/src/stan/optimization/lbfgs_update.hpp +++ b/src/stan/optimization/lbfgs_update.hpp @@ -71,7 +71,8 @@ class LBFGSUpdate { * @param[in] gk Gradient direction. **/ inline void search_direction(VectorT &pk, const VectorT &gk) const { - std::vector alphas(_buf.size()); + std::vector &alphas = _alphas; + alphas.resize(_buf.size()); typename stan::util::ring_buffer::const_reverse_iterator buf_rit; typename stan::util::ring_buffer::const_iterator buf_it; typename std::vector::const_iterator alpha_it; @@ -105,6 +106,8 @@ class LBFGSUpdate { protected: stan::util::ring_buffer _buf; Scalar _gammak; + // workspace of search_direction, kept to avoid an allocation per call + mutable std::vector _alphas; }; } // namespace optimization } // namespace stan From 23662fc6cd11d8a3d7ce3253d9c241f58fc3195d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Mon, 19 Oct 2026 12:24:24 +0200 Subject: [PATCH 2/2] Make LBFGSUpdate::search_direction non-const now that it writes a workspace The two-loop workspace became a mutable member, so a const call was no longer safe to make concurrently. No caller uses a const LBFGSUpdate. --- src/stan/optimization/lbfgs_update.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stan/optimization/lbfgs_update.hpp b/src/stan/optimization/lbfgs_update.hpp index 896f1a7270f..0d7458d733b 100644 --- a/src/stan/optimization/lbfgs_update.hpp +++ b/src/stan/optimization/lbfgs_update.hpp @@ -70,7 +70,7 @@ class LBFGSUpdate { * direction gk. * @param[in] gk Gradient direction. **/ - inline void search_direction(VectorT &pk, const VectorT &gk) const { + inline void search_direction(VectorT &pk, const VectorT &gk) { std::vector &alphas = _alphas; alphas.resize(_buf.size()); typename stan::util::ring_buffer::const_reverse_iterator buf_rit; @@ -107,7 +107,7 @@ class LBFGSUpdate { stan::util::ring_buffer _buf; Scalar _gammak; // workspace of search_direction, kept to avoid an allocation per call - mutable std::vector _alphas; + std::vector _alphas; }; } // namespace optimization } // namespace stan