From 102650dfffc14f7a89a3d2ef66a7fe6b000feee6 Mon Sep 17 00:00:00 2001 From: federico Date: Thu, 24 Sep 2026 17:00:42 -0400 Subject: [PATCH 1/3] Add Problem::is_step_acceptable, a check on steps the line search accepts is_step_valid() runs on every trial point, which makes it the wrong place for checks that are expensive, or that can only veto a step without the energy reflecting why. Add is_step_acceptable(): it is called only once a trial step has passed the descent criteria, and returning false makes the line search keep backtracking. A problem may change its own discretization before rejecting (e.g. refine a quadrature where it found the step invalid, so the energy there is no longer finite). The line search therefore re-evaluates the energy and gradient at x0 after a rejection, so its criteria keep comparing like with like. If the old direction no longer descends on the changed problem, the line search fails through its usual path so the solver picks a new direction; the "energy at x unchanged" assert in that path is relaxed accordingly. Only Backtracking implements the loop, so Armijo, RobustArmijo and ResidualBacktracking inherit it. NoLineSearch never backtracks and does not call the hook. The default returns true, so existing problems are unaffected. Co-Authored-By: Claude Opus 5.5 --- src/polysolve/nonlinear/Problem.hpp | 21 ++++ .../nonlinear/line_search/Backtracking.cpp | 34 +++++- .../nonlinear/line_search/LineSearch.cpp | 6 +- .../nonlinear/line_search/LineSearch.hpp | 5 + tests/test_nonlinear_solver.cpp | 115 ++++++++++++++++++ 5 files changed, 175 insertions(+), 6 deletions(-) diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 92711ef7..746cffaf 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -102,6 +102,27 @@ 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 the + /// line search keeps backtracking. 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). The line + /// search therefore re-evaluates the energy and gradient at x0 after a + /// rejection, so its descent criteria keep comparing like with like. + /// + /// 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/line_search/Backtracking.cpp b/src/polysolve/nonlinear/line_search/Backtracking.cpp index 70f480c8..426c4971 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.cpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.cpp @@ -4,6 +4,8 @@ #include +#include + namespace polysolve::nonlinear::line_search { @@ -23,7 +25,12 @@ namespace polysolve::nonlinear::line_search { double step_size = starting_step_size; - init_compute_descent_step_size(delta_x, old_grad); + // Local copies: a rejection by is_step_acceptable() may change the + // problem, after which these are re-evaluated at x. + double energy0 = old_energy; + TVector grad0 = old_grad; + + init_compute_descent_step_size(delta_x, grad0); for (; step_size > current_min_step_size() && cur_iter < current_max_step_size_iter(); step_size *= step_ratio, ++cur_iter) { @@ -52,12 +59,31 @@ namespace polysolve::nonlinear::line_search continue; } - m_logger.trace("ls it: {} {}: {}, {}: {}", cur_iter, log::delta("E"), new_energy - old_energy, log::delta("x"), objFunc.step_norm(delta_x, norm_type)); + m_logger.trace("ls it: {} {}: {}, {}: {}", cur_iter, log::delta("E"), new_energy - energy0, 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, energy0, grad0, 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, so + // re-evaluate x before comparing further trial points against it. + m_logger.debug("Step of size {} rejected by the problem; re-evaluating and backtracking", step_size); + problem_changed = true; + objFunc.solution_changed(x); + energy0 = objFunc(x); + objFunc.gradient(x, grad0); + if (!std::isfinite(energy0) || !grad0.array().isFinite().all() || delta_x.dot(grad0) >= 0) + { + // delta_x was computed for the previous problem and need not + // descend on the new one. Fail through the caller's usual path + // (which restores x and ends the line search) so the solver + // picks a new direction. + m_logger.debug("Search direction is no longer a descent direction after the problem changed"); + return 0; } + init_compute_descent_step_size(delta_x, grad0); } return step_size; diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index 515a44fb..84e62e05 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -84,6 +84,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)) @@ -167,8 +168,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; diff --git a/src/polysolve/nonlinear/line_search/LineSearch.hpp b/src/polysolve/nonlinear/line_search/LineSearch.hpp index ad9ebbd2..b27ba5be 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.hpp @@ -99,6 +99,11 @@ namespace polysolve::nonlinear::line_search double step_ratio; int cur_iter; + /// Set when a step rejected by Problem::is_step_acceptable() forced a + /// re-evaluation of the problem during this line search: 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..d6999cb1 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,111 @@ 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 makes the line search keep backtracking") + { + problem.acceptable = [](const Problem::TVector &, const Problem::TVector &x1) { return x1(0) <= 1.2; }; + const double step = ls->line_search(x, dx, problem); + CHECK(step > 0); + CHECK(step * dx(0) <= 1.2); + CHECK(problem.n_acceptable_calls >= 2); + } + + SECTION("a rejection may change the problem") + { + // As polyfem does when it refines a quadrature: the first + // rejected step makes the energy beyond it infinite. The next + // trial (x(0) = 1) lies past the wall, so it has to be turned + // down by the line search's own finiteness check -- 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; + }; + 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); + } + + SECTION("a change that leaves no descent direction fails the line search") + { + problem.acceptable = [&](const Problem::TVector &, const Problem::TVector &) { + problem.c = -problem.c; // dx now points uphill + return false; + }; + CHECK(std::isnan(ls->line_search(x, dx, problem))); + CHECK(problem.n_acceptable_calls == 1); + } + } + } +} From 377cda3b512b725204a934e3237bf75f7991f6a7 Mon Sep 17 00:00:00 2001 From: federico Date: Thu, 24 Sep 2026 20:06:41 -0400 Subject: [PATCH 2/3] Take steps with one rounding per coefficient and cap them without fenv A problem that certifies the whole segment x + t*dx (e.g. polyfem's continuous inversion check) can only trust the configuration the solver takes if that configuration is within a known distance of the segment. - step_along(x, s, dx) forms x + s*dx with std::fma, so each coefficient is rounded once. Every trial and accepted step in the line searches and the final solver update goes through it. - product_rounded_down(a, b) replaces the fesetround(FE_DOWNWARD) block that scaled the step by the max step size. That block needed FENV_ACCESS or -frounding-math to be honored; the FMA-based version is correct in any rounding mode and under any optimization. - Tests: step_along matches a single rounding where a separate multiply and add differ; product_rounded_down never exceeds the exact product. Co-Authored-By: Claude Opus 5.5 --- src/polysolve/Utils.hpp | 30 ++++++++++++++ src/polysolve/nonlinear/Solver.cpp | 2 +- .../nonlinear/line_search/Backtracking.cpp | 2 +- .../nonlinear/line_search/LineSearch.cpp | 24 +++++------ tests/test_nonlinear_solver.cpp | 41 +++++++++++++++++++ 5 files changed, 83 insertions(+), 16 deletions(-) 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/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 426c4971..7588bef7 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.cpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.cpp @@ -34,7 +34,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 { diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index 84e62e05..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(); @@ -122,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); } @@ -156,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; @@ -196,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()) @@ -204,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 { @@ -232,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); @@ -244,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/tests/test_nonlinear_solver.cpp b/tests/test_nonlinear_solver.cpp index d6999cb1..420b566e 100644 --- a/tests/test_nonlinear_solver.cpp +++ b/tests/test_nonlinear_solver.cpp @@ -867,3 +867,44 @@ TEST_CASE("line-search-step-acceptable", "[solver][line_search]") } } } + +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); +} From 7ab0a759d82d8870c3e9edbfe5a74414e6803440 Mon Sep 17 00:00:00 2001 From: federico Date: Thu, 1 Oct 2026 10:55:56 -0400 Subject: [PATCH 3/3] End the line search when is_step_acceptable rejects a step A rejection may change the problem's discretization (polyfem refines a quadrature where the step was found invalid), and the search direction was computed for the old one. Rather than re-evaluating the energy and gradient at x and backtracking further along that direction, fail the line search through the caller's usual path, so the solver computes a new direction on the problem as it is now. Each line search thus works on a single, fixed problem, and the re-evaluation and descent re-check go away. Co-Authored-By: Claude Opus 5.5 --- src/polysolve/nonlinear/Problem.hpp | 15 ++++---- .../nonlinear/line_search/Backtracking.cpp | 37 ++++++------------- .../nonlinear/line_search/LineSearch.hpp | 7 ++-- tests/test_nonlinear_solver.cpp | 33 ++++++----------- 4 files changed, 35 insertions(+), 57 deletions(-) diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 746cffaf..026640ae 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -105,16 +105,17 @@ namespace polysolve::nonlinear /// @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 the - /// line search keeps backtracking. 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. + /// 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). The line - /// search therefore re-evaluates the energy and gradient at x0 after a - /// rejection, so its descent criteria keep comparing like with like. + /// 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. diff --git a/src/polysolve/nonlinear/line_search/Backtracking.cpp b/src/polysolve/nonlinear/line_search/Backtracking.cpp index 7588bef7..66582a61 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.cpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.cpp @@ -4,8 +4,6 @@ #include -#include - namespace polysolve::nonlinear::line_search { @@ -25,12 +23,7 @@ namespace polysolve::nonlinear::line_search { double step_size = starting_step_size; - // Local copies: a rejection by is_step_acceptable() may change the - // problem, after which these are re-evaluated at x. - double energy0 = old_energy; - TVector grad0 = old_grad; - - init_compute_descent_step_size(delta_x, grad0); + init_compute_descent_step_size(delta_x, old_grad); for (; step_size > current_min_step_size() && cur_iter < current_max_step_size_iter(); step_size *= step_ratio, ++cur_iter) { @@ -59,31 +52,23 @@ namespace polysolve::nonlinear::line_search continue; } - m_logger.trace("ls it: {} {}: {}, {}: {}", cur_iter, log::delta("E"), new_energy - energy0, log::delta("x"), objFunc.step_norm(delta_x, norm_type)); + 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, energy0, grad0, 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, so - // re-evaluate x before comparing further trial points against it. - m_logger.debug("Step of size {} rejected by the problem; re-evaluating and backtracking", 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; - objFunc.solution_changed(x); - energy0 = objFunc(x); - objFunc.gradient(x, grad0); - if (!std::isfinite(energy0) || !grad0.array().isFinite().all() || delta_x.dot(grad0) >= 0) - { - // delta_x was computed for the previous problem and need not - // descend on the new one. Fail through the caller's usual path - // (which restores x and ends the line search) so the solver - // picks a new direction. - m_logger.debug("Search direction is no longer a descent direction after the problem changed"); - return 0; - } - init_compute_descent_step_size(delta_x, grad0); + return 0; } return step_size; diff --git a/src/polysolve/nonlinear/line_search/LineSearch.hpp b/src/polysolve/nonlinear/line_search/LineSearch.hpp index b27ba5be..2c9de04a 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.hpp @@ -99,9 +99,10 @@ namespace polysolve::nonlinear::line_search double step_ratio; int cur_iter; - /// Set when a step rejected by Problem::is_step_acceptable() forced a - /// re-evaluation of the problem during this line search: its energy at - /// x may then differ from the one the line search started from. + /// 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: diff --git a/tests/test_nonlinear_solver.cpp b/tests/test_nonlinear_solver.cpp index 420b566e..1d64e9b8 100644 --- a/tests/test_nonlinear_solver.cpp +++ b/tests/test_nonlinear_solver.cpp @@ -826,44 +826,35 @@ TEST_CASE("line-search-step-acceptable", "[solver][line_search]") CHECK(problem.n_acceptable_calls == 1); } - SECTION("a rejection makes the line search keep backtracking") + SECTION("a rejection ends the line search") { problem.acceptable = [](const Problem::TVector &, const Problem::TVector &x1) { return x1(0) <= 1.2; }; - const double step = ls->line_search(x, dx, problem); - CHECK(step > 0); - CHECK(step * dx(0) <= 1.2); - CHECK(problem.n_acceptable_calls >= 2); + CHECK(std::isnan(ls->line_search(x, dx, problem))); + CHECK(problem.n_acceptable_calls == 1); } - SECTION("a rejection may change the problem") + SECTION("a rejection may change the problem for the next line search") { - // As polyfem does when it refines a quadrature: the first - // rejected step makes the energy beyond it infinite. The next - // trial (x(0) = 1) lies past the wall, so it has to be turned - // down by the line search's own finiteness check -- only the - // step it finally accepts reaches the hook again. + // 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); } - - SECTION("a change that leaves no descent direction fails the line search") - { - problem.acceptable = [&](const Problem::TVector &, const Problem::TVector &) { - problem.c = -problem.c; // dx now points uphill - return false; - }; - CHECK(std::isnan(ls->line_search(x, dx, problem))); - CHECK(problem.n_acceptable_calls == 1); - } } } }