Add Problem::is_step_acceptable and take steps with one rounding per coefficient

src/polysolve/Utils.hpp

@@ -3,6 +3,9 @@
 #include "Types.hpp"
 
 #include <spdlog/spdlog.h>
+#include <limits>
+#include <cmath>
+#include <cassert>
 
 #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<double>::infinity()) : p;
+    }
+
     namespace log
     {
 #ifdef POLYSOLVE_WITH_UNICODE

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.

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;

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;

src/polysolve/nonlinear/line_search/LineSearch.cpp

@@ -16,8 +16,6 @@
 #include <spdlog/fmt/bundled/color.h>
 #endif
 
-#include <cfenv>
-
 namespace polysolve::nonlinear::line_search
 {
     static constexpr double NaN = std::numeric_limits<double>::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,15 +194,15 @@ 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())
         {
             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;
     }

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)

tests/test_nonlinear_solver.cpp

@@ -3,6 +3,7 @@
 #include <polysolve/nonlinear/Solver.hpp>
 #include <polysolve/nonlinear/BoxConstraintSolver.hpp>
 #include <polysolve/nonlinear/Problem.hpp>
+#include <polysolve/nonlinear/line_search/LineSearch.hpp>
 #include <polysolve/Utils.hpp>
 #include <polysolve/Types.hpp>
 #include <polysolve/linear/Solver.hpp>
@@ -11,6 +12,12 @@
 #include <polysolve/JSONUtils.hpp>
 #include <catch2/catch.hpp>
 
+#include <jse/jse.h>
+#include "nonlinear-solver-spec.hpp"
+
+#include <functional>
+#include <limits>
+
 //////////////////////////////////////////////////////////////////////////
 
 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<double>::infinity();
+        std::function<bool(const TVector &, const TVector &)> 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<double>::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<spdlog::logger> 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);
+}