Optimizer: compute the CubicInterp roots without cancellation

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;

src/test/unit/optimization/bfgs_linesearch_test.cpp

@@ -76,6 +76,21 @@ 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);
+
+  // 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 {
  public:
   double operator()(const Eigen::Matrix<double, Eigen::Dynamic, 1> &x) {