From 5b110dfbcfc114fc1856bb1bd017ccd7803f7023 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 6 Dec 2023 22:30:20 -0500 Subject: [PATCH] Implement RobustArmijo and refactor Armijo LS --- CMakeLists.txt | 2 +- ...er-spec.json => nonlinear-solver-spec.json | 7 +- src/polysolve/nonlinear/Solver.cpp | 2 +- .../nonlinear/descent_strategies/Newton.hpp | 5 +- .../nonlinear/line_search/Armijo.cpp | 70 +++++-------------- .../nonlinear/line_search/Armijo.hpp | 24 ++++--- .../nonlinear/line_search/Backtracking.cpp | 58 ++++++++------- .../nonlinear/line_search/Backtracking.hpp | 17 ++++- .../nonlinear/line_search/CMakeLists.txt | 2 + .../nonlinear/line_search/CppOptArmijo.cpp | 3 +- .../nonlinear/line_search/CppOptArmijo.hpp | 1 + .../nonlinear/line_search/LineSearch.cpp | 55 +++++++++------ .../nonlinear/line_search/LineSearch.hpp | 3 +- .../nonlinear/line_search/MoreThuente.cpp | 3 +- .../nonlinear/line_search/MoreThuente.hpp | 1 + .../nonlinear/line_search/NoLineSearch.cpp | 1 + .../nonlinear/line_search/NoLineSearch.hpp | 1 + .../nonlinear/line_search/RobustArmijo.cpp | 41 +++++++++++ .../nonlinear/line_search/RobustArmijo.hpp | 28 ++++++++ tests/test_nonlinear_solver.cpp | 30 ++++---- 20 files changed, 216 insertions(+), 138 deletions(-) rename non-linear-solver-spec.json => nonlinear-solver-spec.json (99%) create mode 100644 src/polysolve/nonlinear/line_search/RobustArmijo.cpp create mode 100644 src/polysolve/nonlinear/line_search/RobustArmijo.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 93b363b..9fe0065 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -175,7 +175,7 @@ if(POLYSOLVE_LARGE_INDEX) endif() target_compile_definitions(polysolve_linear PRIVATE POLYSOLVE_LINEAR_SPEC="${PROJECT_SOURCE_DIR}/linear-solver-spec.json") -target_compile_definitions(polysolve PRIVATE POLYSOLVE_NON_LINEAR_SPEC="${PROJECT_SOURCE_DIR}/non-linear-solver-spec.json") +target_compile_definitions(polysolve PRIVATE POLYSOLVE_NON_LINEAR_SPEC="${PROJECT_SOURCE_DIR}/nonlinear-solver-spec.json") target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_JSON_SPEC_DIR="${PROJECT_SOURCE_DIR}") diff --git a/non-linear-solver-spec.json b/nonlinear-solver-spec.json similarity index 99% rename from non-linear-solver-spec.json rename to nonlinear-solver-spec.json index 74665d2..67cd68e 100644 --- a/non-linear-solver-spec.json +++ b/nonlinear-solver-spec.json @@ -1,5 +1,4 @@ -[ - { +[{ "pointer": "/", "default": null, "type": "object", @@ -188,6 +187,7 @@ "options": [ "Armijo", "ArmijoAlt", + "RobustArmijo", "Backtracking", "MoreThuente", "None" @@ -247,8 +247,9 @@ }, { "pointer": "/line_search/Armijo/c", - "default": 0.5, + "default": 1e-4, "type": "float", + "min_value": 0, "doc": "Armijo c parameter." }, { diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index daf8811..c34aa39 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -184,7 +184,7 @@ namespace polysolve::nonlinear const auto g_norm_tol = this->m_stop.gradNorm; this->m_stop.gradNorm = first_grad_norm_tol; - StopWatch stop_watch("non-linear solver", this->total_time, m_logger); + StopWatch stop_watch("nonlinear solver", this->total_time, m_logger); stop_watch.start(); m_logger.debug( diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.hpp b/src/polysolve/nonlinear/descent_strategies/Newton.hpp index 74e809c..590fec9 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.hpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.hpp @@ -104,7 +104,10 @@ namespace polysolve::nonlinear const double characteristic_length, spdlog::logger &logger); - std::string name() const override { return internal_name() + "RegularizedNewton"; } + std::string name() const override + { + return fmt::format("{}RegularizedNewton (reg_weight={:g})", internal_name(), reg_weight); + } void reset(const int ndof) override; bool handle_error() override; diff --git a/src/polysolve/nonlinear/line_search/Armijo.cpp b/src/polysolve/nonlinear/line_search/Armijo.cpp index de45452..4c5b551 100644 --- a/src/polysolve/nonlinear/line_search/Armijo.cpp +++ b/src/polysolve/nonlinear/line_search/Armijo.cpp @@ -2,10 +2,6 @@ #include "Armijo.hpp" -#include - -#include - namespace polysolve::nonlinear::line_search { Armijo::Armijo(const json ¶ms, spdlog::logger &logger) @@ -14,58 +10,24 @@ namespace polysolve::nonlinear::line_search c = params["line_search"]["Armijo"]["c"]; } - double Armijo::compute_descent_step_size( - const TVector &x, - const TVector &searchDir, - Problem &objFunc, - const bool use_grad_norm, - const double old_energy_in, - const double starting_step_size) + void Armijo::init_compute_descent_step_size( + const TVector &delta_x, + const TVector &old_grad) { - TVector grad(x.rows()); - double f_in = old_energy_in; - double alpha = starting_step_size; - - bool valid; - - TVector x1 = x + alpha * searchDir; - { - POLYSOLVE_SCOPED_STOPWATCH("constraint set update in LS", constraint_set_update_time, m_logger); - objFunc.solution_changed(x1); - } - - objFunc.gradient(x, grad); - - double f = use_grad_norm ? grad.squaredNorm() : objFunc.value(x1); - const double cache = c * grad.dot(searchDir); - valid = objFunc.is_step_valid(x, x1); - - while ((std::isinf(f) || std::isnan(f) || f > f_in + alpha * cache || !valid) && alpha > current_min_step_size() && cur_iter <= current_max_step_size_iter()) - { - alpha *= step_ratio; - x1 = x + alpha * searchDir; - - { - POLYSOLVE_SCOPED_STOPWATCH("constraint set update in LS", constraint_set_update_time, m_logger); - objFunc.solution_changed(x1); - } - - if (use_grad_norm) - { - objFunc.gradient(x1, grad); - f = grad.squaredNorm(); - } - else - f = objFunc.value(x1); - - valid = objFunc.is_step_valid(x, x1); - - m_logger.trace("ls it: {} f: {} (f_in + alpha * Cache): {} invalid: {} ", cur_iter, f, f_in + alpha * cache, !valid); - - cur_iter++; - } + armijo_criteria = c * delta_x.dot(old_grad); + assert(armijo_criteria <= 0); + } - return alpha; + bool Armijo::criteria( + const TVector &delta_x, + const double old_energy, + const TVector &old_grad, + const double new_energy, + const TVector &new_grad, + const double step_size) const + { + // TODO: Use use_grad_norm + return new_energy <= old_energy + step_size * armijo_criteria; } }; // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/Armijo.hpp b/src/polysolve/nonlinear/line_search/Armijo.hpp index bc3f122..5be85a4 100644 --- a/src/polysolve/nonlinear/line_search/Armijo.hpp +++ b/src/polysolve/nonlinear/line_search/Armijo.hpp @@ -1,13 +1,13 @@ #pragma once -#include "LineSearch.hpp" +#include "Backtracking.hpp" namespace polysolve::nonlinear::line_search { - class Armijo : public LineSearch + class Armijo : public Backtracking { public: - using Superclass = LineSearch; + using Superclass = Backtracking; using typename Superclass::Scalar; using typename Superclass::TVector; @@ -16,15 +16,19 @@ namespace polysolve::nonlinear::line_search virtual std::string name() override { return "Armijo"; } protected: - double compute_descent_step_size( - const TVector &x, + virtual void init_compute_descent_step_size( const TVector &delta_x, - Problem &objFunc, - const bool use_grad_norm, - const double old_energy_in, - const double starting_step_size) override; + const TVector &old_grad) override; + + virtual bool criteria( + const TVector &delta_x, + const double old_energy, + const TVector &old_grad, + const double new_energy, + const TVector &new_grad, + const double step_size) const override; - private: double c; + double armijo_criteria; ///< cached value: c * delta_x.dot(old_grad) }; } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/Backtracking.cpp b/src/polysolve/nonlinear/line_search/Backtracking.cpp index ae7f433..96a9bcd 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 { @@ -20,58 +18,64 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool use_grad_norm, const double old_energy, + const TVector &old_grad, const double starting_step_size) { + assert(!use_grad_norm); double step_size = starting_step_size; - TVector grad(x.rows()); + init_compute_descent_step_size(delta_x, old_grad); - // Find step that reduces the energy - double cur_energy = std::nan(""); - bool is_step_valid = false; - while (step_size > current_min_step_size() && cur_iter < current_max_step_size_iter()) + for (; step_size > current_min_step_size() && cur_iter < current_max_step_size_iter(); step_size *= step_ratio, ++cur_iter) { - TVector new_x = x + step_size * delta_x; + const TVector new_x = x + step_size * delta_x; try { - POLYSOLVE_SCOPED_STOPWATCH("solution changed - constraint set update in LS", this->constraint_set_update_time, m_logger); - objFunc.solution_changed(new_x); + POLYSOLVE_SCOPED_STOPWATCH("solution changed - constraint set update in LS", constraint_set_update_time, m_logger); + objFunc.solution_changed(x); } catch (const std::runtime_error &e) { m_logger.warn("Failed to take step due to \"{}\", reduce step size...", e.what()); - - step_size *= step_ratio; - this->cur_iter++; continue; } - if (use_grad_norm) + if (!objFunc.is_step_valid(x, new_x)) { - objFunc.gradient(new_x, grad); - cur_energy = grad.squaredNorm(); + continue; } - else - cur_energy = objFunc.value(new_x); - is_step_valid = objFunc.is_step_valid(x, new_x); + const double new_energy = objFunc.value(new_x); - m_logger.trace("ls it: {} delta: {} invalid: {} ", this->cur_iter, (cur_energy - old_energy), !is_step_valid); + TVector new_grad; // TODO: only compute gradient if needed + objFunc.gradient(new_x, new_grad); - if (!std::isfinite(cur_energy) || cur_energy >= old_energy || !is_step_valid) + if (!std::isfinite(new_energy) || !std::isfinite(new_grad.norm())) { - step_size *= step_ratio; - // max_step_size should return a collision free step - // assert(objFunc.is_step_collision_free(x, new_x)); + continue; } - else + + m_logger.trace("ls it: {} ΔE: {}", cur_iter, new_energy - old_energy); + + if (criteria(delta_x, old_energy, old_grad, new_energy, new_grad, step_size)) { - break; + break; // found a good step size } - this->cur_iter++; } return step_size; } + + bool Backtracking::criteria( + const TVector &delta_x, + const double old_energy, + const TVector &old_grad, + const double new_energy, + const TVector &new_grad, + const double step_size) const + { + return new_energy < old_energy; + } + } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/Backtracking.hpp b/src/polysolve/nonlinear/line_search/Backtracking.hpp index fec3bf8..a6cd14d 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.hpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.hpp @@ -15,13 +15,26 @@ namespace polysolve::nonlinear::line_search virtual std::string name() override { return "Backtracking"; } - public: double compute_descent_step_size( const TVector &x, const TVector &delta_x, Problem &objFunc, const bool use_grad_norm, - const double old_energy_in, + const double old_energy, + const TVector &old_grad, const double starting_step_size) override; + + protected: + virtual void init_compute_descent_step_size( + const TVector &delta_x, + const TVector &old_grad) {} + + virtual bool criteria( + const TVector &delta_x, + const double old_energy, + const TVector &old_grad, + const double new_energy, + const TVector &new_grad, + const double step_size) const; }; } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/CMakeLists.txt b/src/polysolve/nonlinear/line_search/CMakeLists.txt index 30ef6ee..4533d7a 100644 --- a/src/polysolve/nonlinear/line_search/CMakeLists.txt +++ b/src/polysolve/nonlinear/line_search/CMakeLists.txt @@ -11,6 +11,8 @@ set(SOURCES MoreThuente.hpp NoLineSearch.cpp NoLineSearch.hpp + RobustArmijo.cpp + RobustArmijo.hpp ) source_group(TREE "${CMAKE_CURRENT_SOURCE_DIR}" PREFIX "Source Files" FILES ${SOURCES}) diff --git a/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp b/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp index 25dadc3..f4144d1 100644 --- a/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp +++ b/src/polysolve/nonlinear/line_search/CppOptArmijo.cpp @@ -15,11 +15,12 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool use_grad_norm, const double old_energy, + const TVector &old_grad, const double starting_step_size) { double step_size = cppoptlib::Armijo::linesearch(x, delta_x, objFunc, starting_step_size); // this ensures no collisions and decrease in energy - return after_check.compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, step_size); + return after_check.compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, old_grad, step_size); } } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp b/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp index f9403e9..fc27713 100644 --- a/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp +++ b/src/polysolve/nonlinear/line_search/CppOptArmijo.hpp @@ -23,6 +23,7 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool use_grad_norm, const double old_energy, + const TVector &old_grad, const double starting_step_size) override; private: diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index 416c78d..5cf2eb9 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -2,6 +2,7 @@ #include "Armijo.hpp" #include "Backtracking.hpp" +#include "RobustArmijo.hpp" #include "CppOptArmijo.hpp" #include "MoreThuente.hpp" #include "NoLineSearch.hpp" @@ -16,6 +17,8 @@ namespace polysolve::nonlinear::line_search { + static constexpr double NaN = std::numeric_limits::quiet_NaN(); + std::shared_ptr LineSearch::create(const json ¶ms, spdlog::logger &logger) { const std::string name = params["line_search"]["method"]; @@ -27,6 +30,10 @@ namespace polysolve::nonlinear::line_search { return std::make_shared(params, logger); } + else if (name == "robust_armijo" || name == "RobustArmijo") + { + return std::make_shared(params, logger); + } else if (name == "bisection" || name == "Bisection") { logger.warn("{} linesearch was renamed to \"backtracking\"; using backtracking line-search", name); @@ -55,6 +62,7 @@ namespace polysolve::nonlinear::line_search { return {{"Armijo", "ArmijoAlt", + "RobustArmijo", "Backtracking", "MoreThuente", "None"}}; @@ -81,17 +89,25 @@ namespace polysolve::nonlinear::line_search // ---------------- // Begin linesearch // ---------------- - double old_energy, step_size; + double initial_energy, step_size; + TVector initial_grad; { POLYSOLVE_SCOPED_STOPWATCH("LS begin", m_logger); cur_iter = 0; - old_energy = objFunc.value(x); - if (std::isnan(old_energy)) + initial_energy = objFunc.value(x); + if (std::isnan(initial_energy)) { m_logger.error("Original energy in line search is nan!"); - return std::nan(""); + return NaN; + } + + objFunc.gradient(x, initial_grad); + if (!initial_grad.array().isFinite().all()) + { + m_logger.error("Original gradient in line search is nan!"); + return NaN; } step_size = default_init_step_size; @@ -99,7 +115,6 @@ namespace polysolve::nonlinear::line_search // TODO: removed feature // objFunc.heuristic_max_step(delta_x); } - const double initial_energy = old_energy; // ---------------------------- // Find finite energy step size @@ -108,7 +123,7 @@ namespace polysolve::nonlinear::line_search POLYSOLVE_SCOPED_STOPWATCH("LS compute finite energy step size", checking_for_nan_inf_time, m_logger); step_size = compute_nan_free_step_size(x, delta_x, objFunc, step_size, step_ratio); if (std::isnan(step_size)) - return std::nan(""); + return NaN; } const double nan_free_step_size = step_size; @@ -126,20 +141,16 @@ namespace polysolve::nonlinear::line_search m_logger.trace("Performing narrow-phase CCD"); step_size = compute_collision_free_step_size(x, delta_x, objFunc, step_size); if (std::isnan(step_size)) - return std::nan(""); + return NaN; } const double collision_free_step_size = step_size; - TVector grad(x.rows()); - objFunc.gradient(x, grad); - - if (grad.norm() < 1e-30) + if (initial_grad.norm() < 1e-30) return step_size; - const bool use_grad_norm = grad.norm() < use_grad_norm_tol; - if (use_grad_norm) - old_energy = grad.squaredNorm(); + // TODO: Fix this + const bool use_grad_norm = false; // initial_grad.norm() < use_grad_norm_tol; const double starting_step_size = step_size; // ---------------------- @@ -147,11 +158,11 @@ namespace polysolve::nonlinear::line_search // ---------------------- { POLYSOLVE_SCOPED_STOPWATCH("energy min in LS", classical_line_search_time, m_logger); - step_size = compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, step_size); + step_size = compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, initial_energy, initial_grad, step_size); if (std::isnan(step_size)) { // Superclass::save_sampled_values("failed-line-search-values.csv", x, delta_x, objFunc); - return std::nan(""); + return NaN; } } @@ -163,15 +174,15 @@ namespace polysolve::nonlinear::line_search { m_logger.log(is_final_strategy ? spdlog::level::warn : spdlog::level::debug, "Line search failed to find descent step (f(x)={:g} f(x+αΔx)={:g} α_CCD={:g} α={:g}, ||Δx||={:g} use_grad_norm={} iter={:d})", - old_energy, cur_energy, starting_step_size, + initial_energy, cur_energy, starting_step_size, step_size, delta_x.norm(), use_grad_norm, cur_iter); objFunc.solution_changed(x); -#ifndef NDEBUG + // tolerance for rounding error due to multithreading assert(abs(initial_energy - objFunc.value(x)) < 1e-15); -#endif + objFunc.line_search_end(); - return std::nan(""); + return NaN; } { @@ -220,7 +231,7 @@ namespace polysolve::nonlinear::line_search m_logger.log(is_final_strategy ? spdlog::level::err : spdlog::level::debug, "Line search failed to find a valid finite energy step (cur_iter={:d} step_size={:g})!", cur_iter, step_size); - return std::nan(""); + return NaN; } return step_size; @@ -242,7 +253,7 @@ namespace polysolve::nonlinear::line_search m_logger.log(is_final_strategy ? spdlog::level::err : spdlog::level::debug, "Line search failed because CCD produced a stepsize of zero!"); objFunc.line_search_end(); - return std::nan(""); + return NaN; } { diff --git a/src/polysolve/nonlinear/line_search/LineSearch.hpp b/src/polysolve/nonlinear/line_search/LineSearch.hpp index b1b0b3e..e195955 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.hpp @@ -86,7 +86,8 @@ namespace polysolve::nonlinear::line_search const TVector &delta_x, Problem &objFunc, const bool use_grad_norm, - const double old_energy_in, + const double old_energy, + const TVector &old_grad, const double starting_step_size) = 0; private: diff --git a/src/polysolve/nonlinear/line_search/MoreThuente.cpp b/src/polysolve/nonlinear/line_search/MoreThuente.cpp index bedd766..0c78bf9 100644 --- a/src/polysolve/nonlinear/line_search/MoreThuente.cpp +++ b/src/polysolve/nonlinear/line_search/MoreThuente.cpp @@ -15,10 +15,11 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool use_grad_norm, const double old_energy, + const TVector &old_grad, const double starting_step_size) { const double tmp = cppoptlib::MoreThuente::linesearch(x, delta_x, objFunc, starting_step_size); // this ensures no collisions and decrease in energy - return after_check.compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, std::min(tmp, starting_step_size)); + return after_check.compute_descent_step_size(x, delta_x, objFunc, use_grad_norm, old_energy, old_grad, std::min(tmp, starting_step_size)); } } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/MoreThuente.hpp b/src/polysolve/nonlinear/line_search/MoreThuente.hpp index 27e3d0b..96f7765 100644 --- a/src/polysolve/nonlinear/line_search/MoreThuente.hpp +++ b/src/polysolve/nonlinear/line_search/MoreThuente.hpp @@ -23,6 +23,7 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool use_grad_norm, const double old_energy, + const TVector &old_grad, const double starting_step_size) override; private: diff --git a/src/polysolve/nonlinear/line_search/NoLineSearch.cpp b/src/polysolve/nonlinear/line_search/NoLineSearch.cpp index 8311913..f0ff219 100644 --- a/src/polysolve/nonlinear/line_search/NoLineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/NoLineSearch.cpp @@ -13,6 +13,7 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool, const double, + const TVector &, const double starting_step_size) { objFunc.solution_changed(x + delta_x * starting_step_size); diff --git a/src/polysolve/nonlinear/line_search/NoLineSearch.hpp b/src/polysolve/nonlinear/line_search/NoLineSearch.hpp index 5733394..f7eb101 100644 --- a/src/polysolve/nonlinear/line_search/NoLineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/NoLineSearch.hpp @@ -21,6 +21,7 @@ namespace polysolve::nonlinear::line_search Problem &objFunc, const bool, const double, + const TVector &, const double starting_step_size) override; }; } // namespace polysolve::nonlinear::line_search \ No newline at end of file diff --git a/src/polysolve/nonlinear/line_search/RobustArmijo.cpp b/src/polysolve/nonlinear/line_search/RobustArmijo.cpp new file mode 100644 index 0000000..1734bc4 --- /dev/null +++ b/src/polysolve/nonlinear/line_search/RobustArmijo.cpp @@ -0,0 +1,41 @@ +#include "RobustArmijo.hpp" + +#include + +#include + +namespace polysolve::nonlinear::line_search +{ + RobustArmijo::RobustArmijo(const json ¶ms, spdlog::logger &logger) + : Superclass(params, logger) + { + } + + bool RobustArmijo::criteria( + const TVector &delta_x, + const double old_energy, + const TVector &old_grad, + const double new_energy, + const TVector &new_grad, + const double step_size) const + { + if (new_energy <= old_energy + step_size * this->armijo_criteria) // Try Armijo first + return true; + + if (std::abs(new_energy - old_energy) <= 0.1 * std::abs(old_energy)) + { + // TODO: Compute new_grad here as needed + // objFunc.gradient(new_x, new_grad); + assert(old_grad.size() == new_grad.size()); + + const double deltaE_approx = step_size / 2 * delta_x.dot(new_grad + old_grad); + const double abs_eps_est = step_size / 2 * std::abs(delta_x.dot(new_grad - old_grad)); + + if (deltaE_approx + abs_eps_est <= step_size * this->armijo_criteria) + return true; + } + + return false; + } + +}; // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/RobustArmijo.hpp b/src/polysolve/nonlinear/line_search/RobustArmijo.hpp new file mode 100644 index 0000000..987c0d6 --- /dev/null +++ b/src/polysolve/nonlinear/line_search/RobustArmijo.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include "Armijo.hpp" + +namespace polysolve::nonlinear::line_search +{ + class RobustArmijo : public Armijo + { + public: + using Superclass = Armijo; + using typename Superclass::Scalar; + using typename Superclass::TVector; + + RobustArmijo(const json ¶ms, spdlog::logger &logger); + + virtual std::string name() override { return "RobustArmijo"; } + + protected: + bool criteria( + const TVector &delta_x, + const double old_energy, + const TVector &old_grad, + const double new_energy, + const TVector &new_grad, + const double step_size) const override; + }; + +} // namespace polysolve::nonlinear::line_search diff --git a/tests/test_nonlinear_solver.cpp b/tests/test_nonlinear_solver.cpp index bd28962..f8e1b80 100644 --- a/tests/test_nonlinear_solver.cpp +++ b/tests/test_nonlinear_solver.cpp @@ -240,14 +240,18 @@ class Beale : public AnalyticTestProblem } }; - class InequalityConstraint : public Problem { public: - InequalityConstraint(const double upper_bound): upper_bound_(upper_bound) {} + InequalityConstraint(const double upper_bound) : upper_bound_(upper_bound) {} double value(const TVector &x) override { return x(0) - upper_bound_; } - void gradient(const TVector &x, TVector &gradv) override { gradv.setZero(x.size()); gradv(0) = 1; } + void gradient(const TVector &x, TVector &gradv) override + { + gradv.setZero(x.size()); + gradv(0) = 1; + } void hessian(const TVector &x, THessian &hessian) override {} + private: const double upper_bound_; }; @@ -351,18 +355,18 @@ void test_solvers(const std::vector &solvers, const int iters, cons } } -TEST_CASE("non-linear", "[solver]") +TEST_CASE("nonlinear", "[solver]") { test_solvers(Solver::available_solvers(), 1000, false); // test_solvers({"L-BFGS"}, 1000, false); } -TEST_CASE("non-linear-easier", "[solver]") +TEST_CASE("nonlinear-easier", "[solver]") { test_solvers(Solver::available_solvers(), 5000, true); } -TEST_CASE("non-linear-box-constraint", "[solver]") +TEST_CASE("nonlinear-box-constraint", "[solver]") { std::vector> problems; problems.push_back(std::make_unique()); @@ -395,7 +399,7 @@ TEST_CASE("non-linear-box-constraint", "[solver]") if (ls == "None" && solver_name != "MMA") continue; if (solver_name == "MMA" && ls != "None") - continue; + continue; solver_params["line_search"]["method"] = ls; auto solver = BoxConstraintSolver::create(solver_params, @@ -434,7 +438,7 @@ TEST_CASE("non-linear-box-constraint", "[solver]") } } -TEST_CASE("non-linear-box-constraint-input", "[solver]") +TEST_CASE("nonlinear-box-constraint-input", "[solver]") { std::vector> problems; problems.push_back(std::make_unique()); @@ -461,7 +465,7 @@ TEST_CASE("non-linear-box-constraint-input", "[solver]") if (ls == "None" && solver_name != "MMA") continue; if (solver_name == "MMA" && ls != "None") - continue; + continue; solver_params["line_search"]["method"] = ls; QuadraticProblem::TVector x(prob->size()); @@ -528,13 +532,11 @@ TEST_CASE("MMA", "[solver]") solver_params["solver"] = "MMA"; solver_params["line_search"]["method"] = "None"; - auto solver = BoxConstraintSolver::create(solver_params, - linear_solver_params, - characteristic_length, - *logger); + auto solver = BoxConstraintSolver::create( + solver_params, linear_solver_params, characteristic_length, *logger); auto c = std::make_shared(solver_params["box_constraints"]["bounds"][1]); - dynamic_cast(*solver).add_constraint(c); + dynamic_cast(*solver).add_constraint(c); QuadraticProblem::TVector x(prob->size()); x.setConstant(3);