Implement RobustArmijo and refactor Armijo Line Search

CMakeLists.txt

@@ -173,7 +173,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}")
 
 

non-linear-solver-spec.json → nonlinear-solver-spec.json

@@ -1,5 +1,4 @@
-[
-    {
+[{
         "pointer": "/",
         "default": null,
         "type": "object",
@@ -277,17 +276,19 @@
             "max_step_size_iter_final",
             "default_init_step_size",
             "step_ratio",
-            "Armijo"
+            "Armijo",
+            "RobustArmijo"
         ],
         "doc": "Settings for line-search in the nonlinear solver"
     },
     {
         "pointer": "/line_search/method",
-        "default": "Backtracking",
+        "default": "RobustArmijo",
         "type": "string",
         "options": [
             "Armijo",
             "ArmijoAlt",
+            "RobustArmijo",
             "Backtracking",
             "MoreThuente",
             "None"
@@ -347,10 +348,27 @@
     },
     {
         "pointer": "/line_search/Armijo/c",
-        "default": 0.5,
+        "default": 1e-4,
         "type": "float",
+        "min_value": 0,
         "doc": "Armijo c parameter."
     },
+    {
+        "pointer": "/line_search/RobustArmijo",
+        "default": null,
+        "type": "object",
+        "optional": [
+            "delta_relative_tolerance"
+        ],
+        "doc": "Options for RobustArmijo."
+    },
+    {
+        "pointer": "/line_search/RobustArmijo/delta_relative_tolerance",
+        "default": 0.1,
+        "type": "float",
+        "min_value": 0,
+        "doc": "Relative tolerance on E to switch to approximate."
+    },
     {
         "pointer": "/box_constraints",
         "type": "object",

src/polysolve/nonlinear/Solver.cpp

@@ -218,7 +218,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(

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;

src/polysolve/nonlinear/line_search/Armijo.cpp

@@ -2,10 +2,6 @@
 
 #include "Armijo.hpp"
 
-#include <polysolve/Utils.hpp>
-
-#include <spdlog/spdlog.h>
-
 namespace polysolve::nonlinear::line_search
 {
     Armijo::Armijo(const json &params, spdlog::logger &logger)
@@ -14,58 +10,25 @@ namespace polysolve::nonlinear::line_search
         c = params["line_search"]["Armijo"]["c"];
     }
 
-    double Armijo::compute_descent_step_size(
-        const TVector &x,
-        const TVector &searchDir,
+    void Armijo::init_compute_descent_step_size(
+        const TVector &delta_x,
+        const TVector &old_grad)
+    {
+        armijo_criteria = c * delta_x.dot(old_grad);
+        assert(armijo_criteria <= 0);
+    }
+
+    bool Armijo::criteria(
+        const TVector &delta_x,
         Problem &objFunc,
         const bool use_grad_norm,
-        const double old_energy_in,
-        const double starting_step_size)
+        const double old_energy,
+        const TVector &old_grad,
+        const TVector &new_x,
+        const double new_energy,
+        const double step_size) const
     {
-        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++;
-        }
-
-        return alpha;
+        return new_energy <= old_energy + step_size * armijo_criteria;
     }
 
 }; // namespace polysolve::nonlinear::line_search

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,21 @@ 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,
+            const TVector &old_grad) override;
+
+        virtual bool criteria(
             const TVector &delta_x,
             Problem &objFunc,
             const bool use_grad_norm,
-            const double old_energy_in,
-            const double starting_step_size) override;
+            const double old_energy,
+            const TVector &old_grad,
+            const TVector &new_x,
+            const double new_energy,
+            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

src/polysolve/nonlinear/line_search/Backtracking.cpp

@@ -4,8 +4,6 @@
 
 #include <spdlog/spdlog.h>
 
-#include <cfenv>
-
 namespace polysolve::nonlinear::line_search
 {
 
@@ -20,58 +18,68 @@
         Problem &objFunc,
         const bool use_grad_norm,
         const double old_energy,
+        const TVector &old_grad,
         const double starting_step_size)
     {
         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);
+                POLYSOLVE_SCOPED_STOPWATCH("solution changed - constraint set update in LS", constraint_set_update_time, m_logger);
                 objFunc.solution_changed(new_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);
-
-            if (!std::isfinite(cur_energy) || cur_energy >= old_energy || !is_step_valid)
+            if (!std::isfinite(new_energy))
             {
-                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, objFunc, use_grad_norm, old_energy, old_grad, new_x, new_energy, step_size))
             {
-                break;
+                break; // found a good step size
             }
-            this->cur_iter++;
         }
 
         return step_size;
     }
+
+    bool Backtracking::criteria(
+        const TVector &delta_x,
+        Problem &objFunc,
+        const bool use_grad_norm,
+        const double old_energy,
+        const TVector &old_grad,
+        const TVector &new_x,
+        const double new_energy,
+        const double step_size) const
+    {
+        if (use_grad_norm)
+        {
+            TVector new_grad;
+            objFunc.gradient(new_x, new_grad);
+            return new_grad.norm() < old_grad.norm(); // TODO: cache old_grad.norm()
+        }
+        return new_energy < old_energy;
+    }
+
 } // namespace polysolve::nonlinear::line_search

src/polysolve/nonlinear/line_search/Backtracking.hpp

@@ -15,13 +15,28 @@ 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,
+            Problem &objFunc,
+            const bool use_grad_norm,
+            const double old_energy,
+            const TVector &old_grad,
+            const TVector &new_x,
+            const double new_energy,
+            const double step_size) const;
     };
 } // namespace polysolve::nonlinear::line_search

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})

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<Problem, 1>::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

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: