From b7eac4591b2374ee17e68e890ff1de779692fbe9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Sat, 17 Oct 2026 16:07:09 +0200 Subject: [PATCH 1/2] Spend fewer evaluations in failing Wolfe line searches WolfLSZoom stops after maxLSIts evaluations, failed ones included; a search that cannot succeed used to go on until the bracket collapsed to adjacent doubles. A finite cubic trial step within 1% of a bracket end, or outside the bracket, is moved to 1% inside instead of replaced by bisection, which shrinks the bracket 100-fold rather than 2-fold when the minimiser is close to that end; bisection remains on every fifth zoom iteration, for a non-finite trial, and for a trial at an end of a bracket under about 50 ulps wide, where the margin rounds away. WolfeLineSearch refuses an ascent or NaN direction at once, so that the caller resets to the gradient without first spending evaluations. --- src/stan/optimization/bfgs_linesearch.hpp | 44 ++++++++++++++----- .../optimization/bfgs_linesearch_test.cpp | 15 +++++++ 2 files changed, 48 insertions(+), 11 deletions(-) diff --git a/src/stan/optimization/bfgs_linesearch.hpp b/src/stan/optimization/bfgs_linesearch.hpp index a883845d050..cfcf7ff5461 100644 --- a/src/stan/optimization/bfgs_linesearch.hpp +++ b/src/stan/optimization/bfgs_linesearch.hpp @@ -108,21 +108,26 @@ Scalar CubicInterp(const Scalar &x0, const Scalar &f0, const Scalar &df0, /** * An internal utility function for implementing WolfeLineSearch() + * + * Gives up (returns 1) when the bracket collapses or after max_evals + * function evaluations, failed ones included. **/ template int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, FunctorType &func, const XType &x, const Scalar &f, const Scalar &dfp, const Scalar &c1dfp, const Scalar &c2dfp, const XType &p, Scalar alo, Scalar aloF, Scalar aloDFp, - Scalar ahi, Scalar ahiF, Scalar ahiDFp, - const Scalar &min_range) { + Scalar ahi, Scalar ahiF, Scalar ahiDFp, const Scalar &min_range, + Scalar max_evals = std::numeric_limits::infinity()) { Scalar d1, d2, newDFp; int itNum(0); + int nEvals(0); while (1) { itNum++; - if (std::fabs(alo - ahi) < min_range || std::nextafter(alo, ahi) == ahi) + if (std::fabs(alo - ahi) < min_range || std::nextafter(alo, ahi) == ahi + || nEvals >= max_evals) return 1; if (itNum % 5 == 0) { @@ -135,18 +140,28 @@ int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, d2 = -d2; alpha = ahi - (ahi - alo) * (ahiDFp + d2 - d1) / (ahiDFp - aloDFp + 2 * d2); - if (!std::isfinite(alpha) - || alpha < std::min(alo, ahi) + 0.01 * std::fabs(alo - ahi) - || alpha > std::max(alo, ahi) - 0.01 * std::fabs(alo - ahi)) + if (std::isfinite(alpha)) { + // Keep the trial at least 1% of the bracket away from both ends. When + // the minimizer is very close to one end (a first step far too long + // for the curvature), this shrinks the bracket 100-fold per trial, + // where falling back to bisection would only halve it. + const Scalar margin = 0.01 * std::fabs(alo - ahi); + alpha = std::min(std::max(alpha, std::min(alo, ahi) + margin), + std::max(alo, ahi) - margin); + } + // Bisect instead when the trial is not finite, or when it is an end of + // a bracket so narrow (below about 50 ulps) that the margin rounded away + if (!(std::min(alo, ahi) < alpha && alpha < std::max(alo, ahi))) { alpha = 0.5 * (alo + ahi); + } } newX = x + alpha * p; - while (func(newX, newF, newDF)) { + while (++nEvals, func(newX, newF, newDF)) { const Scalar alpha_prev = alpha; alpha = 0.5 * (alpha + std::min(alo, ahi)); if (std::fabs(std::min(alo, ahi) - alpha) < min_range - || alpha == alpha_prev) + || alpha == alpha_prev || nEvals >= max_evals) return 1; newX = x + alpha * p; } @@ -220,7 +235,8 @@ int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, * * @param minAlpha Smallest allowable step-size. * - * @param maxLSIts Maximum number line search iterations. + * @param maxLSIts Maximum number of line search iterations, both for + * expanding the step and for function evaluations while zooming in. * * @param maxLSRestarts Maximum number of times line search will * restart with \f$ f() \f$ failing. @@ -234,6 +250,12 @@ int WolfeLineSearch(FunctorType &func, Scalar &alpha, XType &x1, const Scalar &c1, const Scalar &c2, const Scalar &minAlpha, const Scalar &maxLSIts, const Scalar &maxLSRestarts) { const Scalar dfp(gradx0.dot(p)); + // The bracketing below assumes that f decreases along p. Refuse an ascent + // (or NaN) direction at once, so that the caller resets to -g instead of + // spending evaluations on a search that cannot bracket a decrease. + if (!(dfp <= 0)) { + return 1; + } const Scalar c1dfp(c1 * dfp); const Scalar c2dfp(c2 * dfp); @@ -271,7 +293,7 @@ int WolfeLineSearch(FunctorType &func, Scalar &alpha, XType &x1, if ((func_val > f0 + alpha * c1dfp) || (func_val >= prevF && nits > 0)) { retCode = WolfLSZoom(alpha, x1, func_val, gradx1, func, x0, f0, dfp, c1dfp, c2dfp, p, alpha0, prevF, prevDFp, alpha1, - func_val, newDFp, 1e-16); + func_val, newDFp, 1e-16, maxLSIts); break; } if (std::fabs(newDFp) <= -c2dfp) { @@ -281,7 +303,7 @@ int WolfeLineSearch(FunctorType &func, Scalar &alpha, XType &x1, if (newDFp >= 0) { retCode = WolfLSZoom(alpha, x1, func_val, gradx1, func, x0, f0, dfp, c1dfp, c2dfp, p, alpha1, func_val, newDFp, alpha0, - prevF, prevDFp, 1e-16); + prevF, prevDFp, 1e-16, maxLSIts); break; } diff --git a/src/test/unit/optimization/bfgs_linesearch_test.cpp b/src/test/unit/optimization/bfgs_linesearch_test.cpp index 42cedef1f7f..d51a0c3c6fd 100644 --- a/src/test/unit/optimization/bfgs_linesearch_test.cpp +++ b/src/test/unit/optimization/bfgs_linesearch_test.cpp @@ -277,3 +277,18 @@ TEST(OptimizationBfgsLinesearch, WolfLSZoom_failing_evaluations) { EXPECT_EQ(1, ret); EXPECT_LT(func.n_evals_, 1000); } + +TEST(OptimizationBfgsLinesearch, WolfLSZoom_max_evals_counts_failures) { + // every evaluation fails: the zoom gives up after max_evals of them + using stan::optimization::WolfLSZoom; + linesearch_jumpfunc func(10.0, 0.0, 20.0); + Eigen::Matrix x0(1), x1, p(1), gradx1; + x0 << 0.0; + p << 1.0; + const double f0 = 0.0, dfp = -1.0; + double alpha = 1.0, f1 = 0.0; + int ret = WolfLSZoom(alpha, x1, f1, gradx1, func, x0, f0, dfp, 1e-4 * dfp, + 0.9 * dfp, p, 0.0, f0, dfp, 1.0, 1.0, 1.0, 1e-16, 5.0); + EXPECT_EQ(1, ret); + EXPECT_EQ(5, func.n_evals_); +} From 42dedb1f6d532085bdea97ec5fc9f8a5f39bf7d3 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 10:32:37 +0300 Subject: [PATCH 2/2] Limit only failed zoom evaluations and bisect when the bracket shrinks slowly --- src/stan/optimization/bfgs_linesearch.hpp | 31 ++++--- .../optimization/bfgs_linesearch_test.cpp | 88 ++++++++++++++++++- 2 files changed, 107 insertions(+), 12 deletions(-) diff --git a/src/stan/optimization/bfgs_linesearch.hpp b/src/stan/optimization/bfgs_linesearch.hpp index cfcf7ff5461..41fd57d993c 100644 --- a/src/stan/optimization/bfgs_linesearch.hpp +++ b/src/stan/optimization/bfgs_linesearch.hpp @@ -109,8 +109,8 @@ Scalar CubicInterp(const Scalar &x0, const Scalar &f0, const Scalar &df0, /** * An internal utility function for implementing WolfeLineSearch() * - * Gives up (returns 1) when the bracket collapses or after max_evals - * function evaluations, failed ones included. + * Gives up (returns 1) when the bracket collapses or after max_failed + * failed function evaluations. **/ template int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, @@ -118,19 +118,29 @@ int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, const Scalar &dfp, const Scalar &c1dfp, const Scalar &c2dfp, const XType &p, Scalar alo, Scalar aloF, Scalar aloDFp, Scalar ahi, Scalar ahiF, Scalar ahiDFp, const Scalar &min_range, - Scalar max_evals = std::numeric_limits::infinity()) { + Scalar max_failed = std::numeric_limits::infinity()) { Scalar d1, d2, newDFp; int itNum(0); - int nEvals(0); + int nFailed(0); + Scalar width_1(std::numeric_limits::infinity()); + Scalar width_2(width_1); while (1) { itNum++; if (std::fabs(alo - ahi) < min_range || std::nextafter(alo, ahi) == ahi - || nEvals >= max_evals) + || nFailed >= max_failed) return 1; - if (itNum % 5 == 0) { + // Bisect when the last two trials did not shrink the bracket to 0.66 of + // its width (More and Thuente, 1994): a trial kept 1% from one end can + // remove only 1% of the bracket when the minimizer is near the other end. + const Scalar width = std::fabs(ahi - alo); + const bool slow = width > 0.66 * width_2; + width_2 = width_1; + width_1 = width; + + if (itNum % 5 == 0 || slow) { alpha = 0.5 * (alo + ahi); } else { // Perform cubic interpolation to determine next point to try @@ -157,11 +167,12 @@ int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, } newX = x + alpha * p; - while (++nEvals, func(newX, newF, newDF)) { + while (func(newX, newF, newDF)) { + ++nFailed; const Scalar alpha_prev = alpha; alpha = 0.5 * (alpha + std::min(alo, ahi)); if (std::fabs(std::min(alo, ahi) - alpha) < min_range - || alpha == alpha_prev || nEvals >= max_evals) + || alpha == alpha_prev || nFailed >= max_failed) return 1; newX = x + alpha * p; } @@ -235,8 +246,8 @@ int WolfLSZoom(Scalar &alpha, XType &newX, Scalar &newF, XType &newDF, * * @param minAlpha Smallest allowable step-size. * - * @param maxLSIts Maximum number of line search iterations, both for - * expanding the step and for function evaluations while zooming in. + * @param maxLSIts Maximum number of line search iterations for expanding + * the step, and of failed function evaluations while zooming in. * * @param maxLSRestarts Maximum number of times line search will * restart with \f$ f() \f$ failing. diff --git a/src/test/unit/optimization/bfgs_linesearch_test.cpp b/src/test/unit/optimization/bfgs_linesearch_test.cpp index d51a0c3c6fd..e7f795fb5da 100644 --- a/src/test/unit/optimization/bfgs_linesearch_test.cpp +++ b/src/test/unit/optimization/bfgs_linesearch_test.cpp @@ -278,8 +278,8 @@ TEST(OptimizationBfgsLinesearch, WolfLSZoom_failing_evaluations) { EXPECT_LT(func.n_evals_, 1000); } -TEST(OptimizationBfgsLinesearch, WolfLSZoom_max_evals_counts_failures) { - // every evaluation fails: the zoom gives up after max_evals of them +TEST(OptimizationBfgsLinesearch, WolfLSZoom_max_failed_evaluations) { + // every evaluation fails: the zoom gives up after max_failed of them using stan::optimization::WolfLSZoom; linesearch_jumpfunc func(10.0, 0.0, 20.0); Eigen::Matrix x0(1), x1, p(1), gradx1; @@ -292,3 +292,87 @@ TEST(OptimizationBfgsLinesearch, WolfLSZoom_max_evals_counts_failures) { EXPECT_EQ(1, ret); EXPECT_EQ(5, func.n_evals_); } + +TEST(OptimizationBfgsLinesearch, wolfeLineSearch_narrow_kink) { + // The slope changes from -1 to 10 over a width of 1e-13 at 1e-5, far + // below the first trial step 1. Only steps inside the kink meet the + // curvature condition. The zoom needs more than maxLSIts = 20 successful + // evaluations to narrow the bracket to the kink, so a limit of maxLSIts on + // all its evaluations makes the search fail. + using stan::optimization::WolfeLineSearch; + const double c1 = 1e-4, c2 = 0.9, m = 1e-5, w = 1e-13, s = 10; + auto func = [&](const Eigen::Matrix &x, double &f, + Eigen::Matrix &g) { + const double z = (x[0] - m) / w; + f = -x[0] + (1 + s) * w * (z > 30 ? z : std::log1p(std::exp(z))); + g.resize(1); + g[0] = -1 + (1 + s) / (1 + std::exp(-z)); + return 0; + }; + Eigen::Matrix x0(1), x1, p(1), gradx0, gradx1; + double f0, f1; + x0 << 0.0; + func(x0, f0, gradx0); + p << 1.0; + double alpha = 1.0; + int ret = WolfeLineSearch(func, alpha, x1, f1, gradx1, p, x0, f0, gradx0, c1, + c2, 1e-12, 20.0, 10.0); + EXPECT_EQ(0, ret); + EXPECT_LE(f1, f0 + c1 * alpha * p.dot(gradx0)); + EXPECT_LE(std::fabs(p.dot(gradx1)), c2 * std::fabs(p.dot(gradx0))); +} + +TEST(OptimizationBfgsLinesearch, wolfeLineSearch_stiff_first_step) { + // phi(x) = -x + 1e11 x^2 / 2: the first trial step 1e-3 is 1e8 times the + // minimizer 1e-11. A cubic trial kept 1% from the bracket end shrinks the + // bracket 100-fold per evaluation; bisection needs 22 evaluations. + using stan::optimization::WolfeLineSearch; + const double c1 = 1e-4, c2 = 0.9; + int n_evals = 0; + auto func = [&](const Eigen::Matrix &x, double &f, + Eigen::Matrix &g) { + ++n_evals; + f = -x[0] + 0.5e11 * x[0] * x[0]; + g.resize(1); + g[0] = -1 + 1e11 * x[0]; + return 0; + }; + Eigen::Matrix x0(1), x1, p(1), gradx0, gradx1; + double f0, f1; + x0 << 0.0; + func(x0, f0, gradx0); + n_evals = 0; + p << 1.0; + double alpha = 1e-3; + int ret = WolfeLineSearch(func, alpha, x1, f1, gradx1, p, x0, f0, gradx0, c1, + c2, 1e-12, 20.0, 10.0); + EXPECT_EQ(0, ret); + EXPECT_LE(f1, f0 + c1 * alpha * p.dot(gradx0)); + EXPECT_LE(std::fabs(p.dot(gradx1)), c2 * std::fabs(p.dot(gradx0))); + EXPECT_LE(n_evals, 10); +} + +TEST(OptimizationBfgsLinesearch, wolfeLineSearch_ascent_direction) { + // p is not a descent direction: the search fails without evaluating f + using stan::optimization::WolfeLineSearch; + int n_evals = 0; + auto func = [&](const Eigen::Matrix &x, double &f, + Eigen::Matrix &g) { + ++n_evals; + f = -x[0]; + g.resize(1); + g[0] = -1; + return 0; + }; + Eigen::Matrix x0(1), x1, p(1), gradx0, gradx1; + double f0, f1; + x0 << 0.0; + func(x0, f0, gradx0); + n_evals = 0; + p << -1.0; + double alpha = 1.0; + int ret = WolfeLineSearch(func, alpha, x1, f1, gradx1, p, x0, f0, gradx0, + 1e-4, 0.9, 1e-16, 20.0, 10.0); + EXPECT_EQ(1, ret); + EXPECT_EQ(0, n_evals); +}