diff --git a/src/stan/optimization/bfgs_linesearch.hpp b/src/stan/optimization/bfgs_linesearch.hpp index a883845d050..41fd57d993c 100644 --- a/src/stan/optimization/bfgs_linesearch.hpp +++ b/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 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::infinity()) { Scalar d1, d2, newDFp; int itNum(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) + 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; } diff --git a/src/test/unit/optimization/bfgs_linesearch_test.cpp b/src/test/unit/optimization/bfgs_linesearch_test.cpp index 42cedef1f7f..e7f795fb5da 100644 --- a/src/test/unit/optimization/bfgs_linesearch_test.cpp +++ b/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 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 &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); +}