diff --git a/src/polysolve/Utils.hpp b/src/polysolve/Utils.hpp index e6141a1a..1f2e6eb9 100644 --- a/src/polysolve/Utils.hpp +++ b/src/polysolve/Utils.hpp @@ -3,6 +3,9 @@ #include "Types.hpp" #include +#include +#include +#include #define POLYSOLVE_SCOPED_STOPWATCH(...) polysolve::StopWatch __polysolve_stopwatch(__VA_ARGS__) @@ -62,6 +65,33 @@ namespace polysolve double extract_param(const std::string &key, const std::string &name, const json &json); + /// @brief x + s * dx, each coefficient rounded once (fused multiply-add). + /// + /// Every place the nonlinear solver forms a trial or accepted step goes + /// through this, so the configuration it finally takes is bit-for-bit the + /// one it checked, and differs from the exact point on the ray by at most + /// half an ulp per coefficient. That is what lets a problem certify a + /// whole segment and trust any point the line search picks on it. + /// std::fma rounds once on every platform (in hardware where available). + inline Eigen::VectorXd step_along(const Eigen::VectorXd &x, const double s, const Eigen::VectorXd &dx) + { + assert(x.size() == dx.size()); + Eigen::VectorXd out(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) + out[i] = std::fma(s, dx[i], x[i]); + return out; + } + + /// @brief The largest double not above a * b (a, b >= 0), in any rounding + /// mode: the fused multiply-add gives the exact rounding error of the + /// product, and when the rounded product overshoots it steps down one + /// representable value. + inline double product_rounded_down(const double a, const double b) + { + const double p = a * b; + return std::fma(a, b, -p) < 0 ? std::nextafter(p, -std::numeric_limits::infinity()) : p; + } + namespace log { #ifdef POLYSOLVE_WITH_UNICODE diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 92711ef7..026640ae 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -102,6 +102,28 @@ namespace polysolve::nonlinear virtual bool after_line_search_custom_operation(const TVector &x0, const TVector &x1) { return false; } + /// @brief Last word on a step the line search is about to accept. + /// + /// Called once a trial step has passed the line search's descent + /// criteria, before it is accepted. Returning false rejects it and + /// ends the line search as failed, so the solver computes a new + /// direction. Unlike is_step_valid(), which runs on every trial point, + /// this runs only on steps that would otherwise be taken, so it suits + /// checks too expensive to run per trial. + /// + /// Before returning false the problem may change its own + /// discretization (e.g. refine a quadrature where the step was found + /// invalid, so that the energy there is no longer finite). Ending the + /// line search keeps each one on a single, fixed problem: the change + /// only affects the directions and line searches that follow. + /// + /// Line searches that never backtrack (NoLineSearch) cannot honor a + /// rejection and do not call this. + /// @param x0 Current solution. + /// @param x1 Step the line search is about to accept. + /// @return True to accept the step, false to reject it. + virtual bool is_step_acceptable(const TVector &x0, const TVector &x1) { return true; } + /// @brief Callback function used to determine if the solver should stop. /// @param state Current state of the solver. /// @param x Current solution. diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index b7c3e8b3..1277e819 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -493,7 +493,7 @@ namespace polysolve::nonlinear } { - TVector x1 = x + rate * delta_x; + TVector x1 = step_along(x, rate, delta_x); if (objFunc.after_line_search_custom_operation(x, x1)) objFunc.solution_changed(x1); x = x1; diff --git a/src/polysolve/nonlinear/line_search/Backtracking.cpp b/src/polysolve/nonlinear/line_search/Backtracking.cpp index 70f480c8..66582a61 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.cpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.cpp @@ -27,7 +27,7 @@ namespace polysolve::nonlinear::line_search for (; step_size > current_min_step_size() && cur_iter < current_max_step_size_iter(); step_size *= step_ratio, ++cur_iter) { - const TVector new_x = x + step_size * delta_x; + const TVector new_x = step_along(x, step_size, delta_x); try { @@ -54,10 +54,21 @@ namespace polysolve::nonlinear::line_search m_logger.trace("ls it: {} {}: {}, {}: {}", cur_iter, log::delta("E"), new_energy - old_energy, log::delta("x"), objFunc.step_norm(delta_x, norm_type)); - if (criteria(delta_x, objFunc, use_grad_norm, old_energy, old_grad, new_x, new_energy, step_size)) - { + if (!criteria(delta_x, objFunc, use_grad_norm, old_energy, old_grad, new_x, new_energy, step_size)) + continue; + + if (objFunc.is_step_acceptable(x, new_x)) break; // found a good step size - } + + // Rejected. The problem may have changed its discretization, and + // delta_x was computed for the old one, so do not keep searching + // along it: fail through the caller's usual path (which restores + // x and ends the line search), and the solver computes a new + // direction on the problem as it is now. Each line search thus + // works on a single, fixed problem. + m_logger.debug("Step of size {} rejected by the problem; ending the line search", step_size); + problem_changed = true; + return 0; } return step_size; diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index 515a44fb..6efb56a6 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -16,8 +16,6 @@ #include #endif -#include - namespace polysolve::nonlinear::line_search { static constexpr double NaN = std::numeric_limits::quiet_NaN(); @@ -84,6 +82,7 @@ namespace polysolve::nonlinear::line_search POLYSOLVE_SCOPED_STOPWATCH("LS begin", m_logger); cur_iter = 0; + problem_changed = false; initial_energy = objFunc(x); if (std::isnan(initial_energy)) @@ -121,7 +120,7 @@ namespace polysolve::nonlinear::line_search // ----------------------------- { POLYSOLVE_SCOPED_STOPWATCH("Line Search Begin - CCD broad-phase", broad_phase_ccd_time, m_logger); - TVector new_x = x + step_size * delta_x; + TVector new_x = step_along(x, step_size, delta_x); objFunc.line_search_begin(x, new_x); } @@ -155,7 +154,7 @@ namespace polysolve::nonlinear::line_search } } - const double cur_energy = objFunc(x + step_size * delta_x); + const double cur_energy = objFunc(step_along(x, step_size, delta_x)); const double descent_step_size = step_size; @@ -167,8 +166,9 @@ namespace polysolve::nonlinear::line_search step_size, objFunc.step_norm(delta_x, norm_type), use_grad_norm, cur_iter); objFunc.solution_changed(x); - // tolerance for rounding error due to multithreading - assert(abs(initial_energy - objFunc(x)) < 1e-15); + // tolerance for rounding error due to multithreading; the energy at + // x only has to match if the problem did not change underneath us + assert(problem_changed || abs(initial_energy - objFunc(x)) < 1e-15); objFunc.line_search_end(); return NaN; @@ -194,7 +194,7 @@ namespace polysolve::nonlinear::line_search const double rate) { double step_size = starting_step_size; - TVector new_x = x + step_size * delta_x; + TVector new_x = step_along(x, step_size, delta_x); // Find step that does not result in nan or infinite energy while (step_size > current_min_step_size() && cur_iter < current_max_step_size_iter()) @@ -202,7 +202,7 @@ namespace polysolve::nonlinear::line_search if (!objFunc.is_step_valid(x, new_x) || !std::isfinite(objFunc(new_x))) { step_size *= rate; - new_x = x + step_size * delta_x; + new_x = step_along(x, step_size, delta_x); } else { @@ -230,7 +230,7 @@ namespace polysolve::nonlinear::line_search const double starting_step_size) { double step_size = starting_step_size; - TVector new_x = x + step_size * delta_x; + TVector new_x = step_along(x, step_size, delta_x); // Find step that is collision free double max_step_size = objFunc.max_step_size(x, new_x); @@ -242,13 +242,11 @@ namespace polysolve::nonlinear::line_search return NaN; } - { - // #pragma STDC FENV_ACCESS ON - const int current_round = std::fegetround(); - std::fesetround(FE_DOWNWARD); - step_size *= max_step_size; // TODO: check me if correct - std::fesetround(current_round); - } + // Round the capped step down, so it never exceeds the certified + // fraction of the segment. Unlike switching the FPU to FE_DOWNWARD + // around the multiply, this does not depend on the compiler honoring + // the rounding mode (no FENV_ACCESS or -frounding-math needed). + step_size = product_rounded_down(step_size, max_step_size); return step_size; } diff --git a/src/polysolve/nonlinear/line_search/LineSearch.hpp b/src/polysolve/nonlinear/line_search/LineSearch.hpp index ad9ebbd2..2c9de04a 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.hpp @@ -99,6 +99,12 @@ namespace polysolve::nonlinear::line_search double step_ratio; int cur_iter; + /// Set when Problem::is_step_acceptable() rejected a step, which ends + /// the line search. The problem may have changed in the process (e.g. + /// refined its discretization), so its energy at x may then differ + /// from the one the line search started from. + bool problem_changed = false; + private: /// @brief Compute step size that avoids nan/infinite energy /// @param x Current input (n x 1) diff --git a/tests/test_nonlinear_solver.cpp b/tests/test_nonlinear_solver.cpp index 805ba7d3..1d64e9b8 100644 --- a/tests/test_nonlinear_solver.cpp +++ b/tests/test_nonlinear_solver.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -11,6 +12,12 @@ #include #include +#include +#include "nonlinear-solver-spec.hpp" + +#include +#include + ////////////////////////////////////////////////////////////////////////// using namespace polysolve; @@ -752,3 +759,143 @@ TEST_CASE("iteration-callback", "[solver][callback]") problem.gradient(x, g); CHECK(g.norm() > 1e-7); } + +namespace +{ + // f(x) = |x - c|^2, minimised at c. From x = 0 the Newton step is c + // itself, so a line search along it accepts the full step unless the + // problem objects. The hook and an optional "wall" (infinite energy past + // x(0) > wall) let each case below object in a different way. + class AcceptanceProblem : public Problem + { + public: + TVector c = (TVector(2) << 2., 1.).finished(); + double wall = std::numeric_limits::infinity(); + std::function acceptable = [](const TVector &, const TVector &) { return true; }; + int n_acceptable_calls = 0; + + double value(const TVector &x) override + { + if (x(0) > wall) + return std::numeric_limits::infinity(); + return (x - c).squaredNorm(); + } + void gradient(const TVector &x, TVector &grad) override { grad = 2 * (x - c); } + void hessian(const TVector &x, THessian &hessian) override + { + hessian = sparse_identity(x.size(), x.size()); + hessian *= 2; + } + bool is_step_acceptable(const TVector &x0, const TVector &x1) override + { + ++n_acceptable_calls; + return acceptable(x0, x1); + } + }; + + json line_search_params(const std::string &method) + { + json params; + params["line_search"]["method"] = method; + jse::JSE jse; + const json rules = jse::embed::non_linear_spec::solver_spec::spec(); + REQUIRE(jse.verify_json(params, rules)); + return jse.inject_defaults(params, rules); + } +} // namespace + +TEST_CASE("line-search-step-acceptable", "[solver][line_search]") +{ + static std::shared_ptr logger = spdlog::stdout_color_mt("test-logger-acceptable"); + logger->set_level(spdlog::level::warn); + + // Every line search that backtracks; NoLineSearch never backtracks and so + // cannot honor a rejection. + for (const std::string method : {"Backtracking", "Armijo", "RobustArmijo"}) + { + DYNAMIC_SECTION(method) + { + const auto ls = line_search::LineSearch::create(line_search_params(method), *logger); + AcceptanceProblem problem; + const Problem::TVector x = Problem::TVector::Zero(2); + const Problem::TVector dx = problem.c - x; + + SECTION("default accepts the full step and is asked once") + { + CHECK(ls->line_search(x, dx, problem) == 1.0); + CHECK(problem.n_acceptable_calls == 1); + } + + SECTION("a rejection ends the line search") + { + problem.acceptable = [](const Problem::TVector &, const Problem::TVector &x1) { return x1(0) <= 1.2; }; + CHECK(std::isnan(ls->line_search(x, dx, problem))); + CHECK(problem.n_acceptable_calls == 1); + } + + SECTION("a rejection may change the problem for the next line search") + { + // As polyfem does when it refines a quadrature: the rejected + // step makes the energy beyond it infinite. That line search + // ends; the next one, on the changed problem, turns down the + // trials past the wall through its own finiteness check, so + // only the step it finally accepts reaches the hook again. + problem.acceptable = [&](const Problem::TVector &, const Problem::TVector &x1) { + if (x1(0) <= 0.8) + return true; + problem.wall = 0.8; + return false; + }; + CHECK(std::isnan(ls->line_search(x, dx, problem))); + CHECK(problem.n_acceptable_calls == 1); + + const double step = ls->line_search(x, dx, problem); + CHECK(step > 0); + CHECK(step * dx(0) <= 0.8); + CHECK(std::isfinite(problem.value(x + step * dx))); + CHECK(problem.n_acceptable_calls == 2); + } + } + } +} + +TEST_CASE("step-along-single-rounding", "[solver][line_search]") +{ + // Matches the fused multiply-add coefficient by coefficient... + Eigen::VectorXd x = Eigen::VectorXd::Random(50), dx = Eigen::VectorXd::Random(50); + const double s = 0.3712; + const Eigen::VectorXd y = step_along(x, s, dx); + for (Eigen::Index i = 0; i < x.size(); ++i) + CHECK(y[i] == std::fma(s, dx[i], x[i])); + + // ...which is not what multiply-then-add gives: here the product rounds, + // and cancelling against x exposes the lost bits. + const double e = std::ldexp(1., -52); + Eigen::VectorXd x1(1), dx1(1); + x1 << -(1 + 2 * e); + dx1 << 1 + e; + // volatile keeps the compiler from contracting this into an FMA itself + // (clang does so by default within an expression). + volatile const double product = (1 + e) * dx1[0]; + const double sep = x1[0] + product; + CHECK(sep == 0.); + CHECK(step_along(x1, 1 + e, dx1)[0] == std::ldexp(1., -104)); +} + +TEST_CASE("product-rounded-down", "[solver][line_search]") +{ + srand(0); + for (int i = 0; i < 10000; ++i) + { + const double a = std::abs(Eigen::VectorXd::Random(1)[0]); + const double b = std::abs(Eigen::VectorXd::Random(1)[0]); + const double p = product_rounded_down(a, b); + // Never above the exact product... + CHECK(std::fma(a, b, -p) >= 0); + // ...and the next double up is. + CHECK(std::fma(a, b, -std::nextafter(p, 2.)) < 0); + } + // Exact products are left alone, e.g. the power-of-two steps + // backtracking produces with the default ratio. + CHECK(product_rounded_down(0.25, 0.8125) == 0.203125); +}