From 4d663bbc0fd48d2ce792f8e95ecaedc2f61745c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Sat, 17 Oct 2026 16:07:06 +0200 Subject: [PATCH 1/2] Compute the CubicInterp roots without cancellation CubicInterp found the stationary points of the interpolating cubic as -(c2 +- sqrt(c2^2 - 2 c1 c3)) / c3. For nearly quadratic data c3 is rounding noise and one of the two roots cancels: CubicInterp returned an end of the interval, or a point up to 28% off, instead of the minimiser. The roots are now -q / c3 and -2 c1 / q with q = c2 + sign(c2) sqrt(...). --- src/stan/optimization/bfgs_linesearch.hpp | 8 ++++++-- src/test/unit/optimization/bfgs_linesearch_test.cpp | 11 +++++++++++ 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/src/stan/optimization/bfgs_linesearch.hpp b/src/stan/optimization/bfgs_linesearch.hpp index a883845d050..bfce16d9992 100644 --- a/src/stan/optimization/bfgs_linesearch.hpp +++ b/src/stan/optimization/bfgs_linesearch.hpp @@ -38,9 +38,13 @@ Scalar CubicInterp(const Scalar &df0, const Scalar &x1, const Scalar &f1, const Scalar c2(-(4 * df0 + 2 * df1) / x1 + 6 * f1 / (x1 * x1)); const Scalar &c1(df0); + // Roots of g'(x) = c3 x^2 / 2 + c2 x + c1 in the cancellation-free form: + // for (near-)quadratic data c3 is rounding noise, and the form that + // subtracts, -(c2 - sign(c2) t_s) / c3, would lose the root -c1 / c2. const Scalar t_s = std::sqrt(c2 * c2 - 2.0 * c1 * c3); - const Scalar s1 = -(c2 + t_s) / c3; - const Scalar s2 = -(c2 - t_s) / c3; + const Scalar q = c2 + std::copysign(t_s, c2); + const Scalar s1 = -q / c3; + const Scalar s2 = -2.0 * c1 / q; Scalar tmpF; Scalar minF, minX; diff --git a/src/test/unit/optimization/bfgs_linesearch_test.cpp b/src/test/unit/optimization/bfgs_linesearch_test.cpp index 42cedef1f7f..4eccb286288 100644 --- a/src/test/unit/optimization/bfgs_linesearch_test.cpp +++ b/src/test/unit/optimization/bfgs_linesearch_test.cpp @@ -76,6 +76,17 @@ TEST(OptimizationBfgsLinesearch, CubicInterp_6arg) { } } +TEST(OptimizationBfgsLinesearch, CubicInterp_quadratic_data) { + // Data at 0 and x1 = 0.4 from phi(x) = -x + x^2 / (2 m), minimised at + // m = 0.1. The cubic coefficient is rounding noise (about -1e-14), and the + // root -(c2 - sqrt(c2^2 - 2 c1 c3)) / c3 cancelled: it came out as 0.128. + using stan::optimization::CubicInterp; + const double m = 0.1, x1 = 0.4; + const double f1 = -x1 + x1 * x1 / (2 * m); + const double df1 = -1 + x1 / m; + EXPECT_NEAR(m, CubicInterp(-1.0, x1, f1, df1, 1e-12, 1.0), 1e-13); +} + class linesearch_testfunc { public: double operator()(const Eigen::Matrix &x) { From d4506dfe99c25523d0dff5877643595f92ec5e22 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= Date: Mon, 19 Oct 2026 11:23:37 +0200 Subject: [PATCH 2/2] Test CubicInterp on exactly quadratic data The existing test relies on c3 being rounding noise, whose value FMA contraction or another compiler can change. Data from -x + x^2 at x1 = 1 make every coefficient exact and c3 = 0 on any platform; develop returned the lower bound there, since the root was 0 / 0. --- src/test/unit/optimization/bfgs_linesearch_test.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/test/unit/optimization/bfgs_linesearch_test.cpp b/src/test/unit/optimization/bfgs_linesearch_test.cpp index 4eccb286288..808dfebecb9 100644 --- a/src/test/unit/optimization/bfgs_linesearch_test.cpp +++ b/src/test/unit/optimization/bfgs_linesearch_test.cpp @@ -85,6 +85,10 @@ TEST(OptimizationBfgsLinesearch, CubicInterp_quadratic_data) { const double f1 = -x1 + x1 * x1 / (2 * m); const double df1 = -1 + x1 / m; EXPECT_NEAR(m, CubicInterp(-1.0, x1, f1, df1, 1e-12, 1.0), 1e-13); + + // m = 0.5, x1 = 1: f1 = 0 and df1 = 1, every coefficient is exact and + // c3 = 0 on any platform. The root was 0 / 0 and the lower bound came out. + EXPECT_DOUBLE_EQ(0.5, CubicInterp(-1.0, 1.0, 0.0, 1.0, 1e-12, 1.0)); } class linesearch_testfunc {