Optimizer: spend fewer evaluations in failing Wolfe line searches

src/stan/optimization/bfgs_linesearch.hpp

@@ -108,24 +108,39 @@ 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_failed
+ * failed function evaluations.
  **/
 template <typename FunctorType, typename Scalar, typename XType>
 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_failed = std::numeric_limits<Scalar>::infinity()) {
   Scalar d1, d2, newDFp;
   int itNum(0);
+  int nFailed(0);
+  Scalar width_1(std::numeric_limits<Scalar>::infinity());
+  Scalar width_2(width_1);
 
   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
+        || 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
@@ -135,18 +150,29 @@ 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)) {
+      ++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)
+          || alpha == alpha_prev || nFailed >= max_failed)
         return 1;
       newX = x + alpha * p;
     }
@@ -220,7 +246,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 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.
@@ -234,6 +261,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 +304,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 +314,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;
     }
 

src/test/unit/optimization/bfgs_linesearch_test.cpp

@@ -277,3 +277,102 @@ TEST(OptimizationBfgsLinesearch, WolfLSZoom_failing_evaluations) {
   EXPECT_EQ(1, ret);
   EXPECT_LT(func.n_evals_, 1000);
 }
+
+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<double, -1, 1> 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_);
+}
+
+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<double, -1, 1> &x, double &f,
+                  Eigen::Matrix<double, -1, 1> &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<double, -1, 1> 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<double, -1, 1> &x, double &f,
+                  Eigen::Matrix<double, -1, 1> &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<double, -1, 1> 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<double, -1, 1> &x, double &f,
+                  Eigen::Matrix<double, -1, 1> &g) {
+    ++n_evals;
+    f = -x[0];
+    g.resize(1);
+    g[0] = -1;
+    return 0;
+  };
+  Eigen::Matrix<double, -1, 1> 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);
+}