Merge pull request #43 from polyfem/refactoring

separated solver from descent strategy

non-linear-solver-spec.json

@@ -13,6 +13,7 @@
             "line_search",
             "allow_out_of_iterations",
             "LBFGS",
+            "LBFGSB",
             "Newton",
             "box_constraints"
         ],
@@ -97,19 +98,41 @@
         "type": "int",
         "doc": "The number of corrections to approximate the inverse Hessian matrix."
     },
+    {
+        "pointer": "/LBFGSB",
+        "default": null,
+        "type": "object",
+        "optional": [
+            "history_size"
+        ],
+        "doc": "Options for the boxed LBFGS."
+    },
+    {
+        "pointer": "/LBFGSB/history_size",
+        "default": 6,
+        "type": "int",
+        "doc": "The number of corrections to approximate the inverse Hessian matrix."
+    },
     {
         "pointer": "/Newton",
         "default": null,
         "type": "object",
         "optional": [
+            "residual_tolerance",
             "reg_weight_min",
             "reg_weight_max",
             "reg_weight_inc",
-            "reg_weight_dec",
-            "force_psd_projection"
+            "force_psd_projection",
+            "use_psd_projection"
         ],
         "doc": "Options for Newton."
     },
+    {
+        "pointer": "/Newton/residual_tolerance",
+        "default": 1e-5,
+        "type": "float",
+        "doc": "Tolerance of the linear system residual. If residual is above, the direction is rejected."
+    },
     {
         "pointer": "/Newton/reg_weight_min",
         "default": 1e-8,
@@ -128,18 +151,18 @@
         "type": "float",
         "doc": "Regulariztion weight increment."
     },
-    {
-        "pointer": "/Newton/reg_weight_dec",
-        "default": 2,
-        "type": "float",
-        "doc": "Regulariztion weight decrement."
-    },
     {
         "pointer": "/Newton/force_psd_projection",
         "default": false,
         "type": "bool",
         "doc": "Force the Hessian to be PSD when using second order solvers (i.e., Newton's method)."
     },
+    {
+        "pointer": "/Newton/use_psd_projection",
+        "default": true,
+        "type": "bool",
+        "doc": "Use PSD as fallback using second order solvers (i.e., Newton's method)."
+    },
     {
         "pointer": "/line_search",
         "default": null,

src/polysolve/nonlinear/BoxConstraintSolver.cpp

@@ -1,7 +1,7 @@
 #include "BoxConstraintSolver.hpp"
 
-#include "LBFGSB.hpp"
-#include "MMA.hpp"
+#include "descent_strategies/box_constraints/LBFGSB.hpp"
+#include "descent_strategies/box_constraints/MMA.hpp"
 
 #include <jse/jse.h>
 
@@ -39,17 +39,28 @@ namespace polysolve::nonlinear
 
         solver_params = jse.inject_defaults(solver_params, rules);
 
-        const std::string solver = solver_params["solver"];
+        const std::string solver_name = solver_params["solver"];
 
-        if (solver == "LBFGSB" || solver == "L-BFGS-B")
+        auto solver = std::make_unique<BoxConstraintSolver>(solver_name, solver_params, characteristic_length, logger);
+
+        if (solver_name == "LBFGSB" || solver_name == "L-BFGS-B")
         {
-            return std::make_unique<LBFGSB>(solver_params, characteristic_length, logger);
+            solver->add_strategy(std::make_unique<LBFGSB>(
+                solver_params, characteristic_length, logger));
         }
-        else if (solver == "MMA")
+        else if (solver_name == "MMA")
         {
-            return std::make_unique<MMA>(solver_params, characteristic_length, logger);
+            if (solver->line_search()->name() != "None")
+                log_and_throw_error(logger, "Invalid linesearch for MMA; MMA requires 'None' linesearch, instead got {}", solver->line_search()->name());
+
+            solver->add_strategy(std::make_unique<MMA>(
+                solver_params, characteristic_length, logger));
         }
-        throw std::runtime_error("Unrecognized solver type: " + solver);
+        else
+            throw std::runtime_error("Unrecognized solver type: " + solver_name);
+
+        solver->set_strategies_iterations(solver_params);
+        return solver;
     }
 
     std::vector<std::string> BoxConstraintSolver::available_solvers()
@@ -58,10 +69,11 @@ namespace polysolve::nonlinear
                 "MMA"};
     }
 
-    BoxConstraintSolver::BoxConstraintSolver(const json &solver_params,
+    BoxConstraintSolver::BoxConstraintSolver(const std::string &name,
+                                             const json &solver_params,
                                              const double characteristic_length,
                                              spdlog::logger &logger)
-        : Superclass(solver_params, characteristic_length, logger)
+        : Superclass(name, solver_params, characteristic_length, logger)
     {
         json box_constraint_params = solver_params["box_constraints"];
         if (box_constraint_params["max_change"] > 0)
@@ -103,6 +115,21 @@ namespace polysolve::nonlinear
         }
     }
 
+    bool BoxConstraintSolver::compute_update_direction(
+        Problem &objFunc,
+        const TVector &x,
+        const TVector &grad,
+        TVector &direction)
+    {
+        const TVector lower_bound = get_lower_bound(x);
+        const TVector upper_bound = get_upper_bound(x);
+
+        return m_strategies[m_descent_strategy]->compute_boxed_update_direction(
+            objFunc, x, grad,
+            lower_bound, upper_bound,
+            direction);
+    }
+
     double BoxConstraintSolver::compute_grad_norm(const Eigen::VectorXd &x,
                                                   const Eigen::VectorXd &grad) const
     {

src/polysolve/nonlinear/BoxConstraintSolver.hpp

@@ -1,10 +1,12 @@
 #pragma once
 
 #include "Solver.hpp"
+#include "descent_strategies/box_constraints/BoxedDescentStrategy.hpp"
 #include <polysolve/Utils.hpp>
 
 namespace polysolve::nonlinear
 {
+
     class BoxConstraintSolver : public Solver
     {
     public:
@@ -24,10 +26,9 @@ namespace polysolve::nonlinear
         static std::vector<std::string> available_solvers();
 
         using Superclass = Solver;
-        using typename Superclass::Scalar;
-        using typename Superclass::TVector;
 
-        BoxConstraintSolver(const json &solver_params,
+        BoxConstraintSolver(const std::string &name,
+                            const json &solver_params,
                             const double characteristic_length,
                             spdlog::logger &logger);
 
@@ -36,8 +37,22 @@ namespace polysolve::nonlinear
         Eigen::VectorXd get_upper_bound(const Eigen::VectorXd &x, bool consider_max_change = true) const;
         Eigen::VectorXd get_max_change(const Eigen::VectorXd &x) const;
 
+        void add_strategy(const std::shared_ptr<BoxedDescentStrategy> &s)
+        {
+            Superclass::add_strategy(s);
+            m_strategies.push_back(s);
+        }
+
+    protected:
+        bool compute_update_direction(
+            Problem &objFunc,
+            const TVector &x,
+            const TVector &grad,
+            TVector &direction) override;
+
     private:
         Eigen::MatrixXd bounds_;
+        std::vector<std::shared_ptr<BoxedDescentStrategy>> m_strategies;
 
         double max_change_val_ = 0;
         Eigen::VectorXd max_change_;

src/polysolve/nonlinear/CMakeLists.txt

@@ -1,26 +1,8 @@
 set(SOURCES
-	LBFGS.hpp
-	LBFGS.cpp
-	LBFGSB.cpp
-	LBFGSB.hpp
-	BFGS.cpp
-	BFGS.hpp
-	GradientDescent.cpp
-	GradientDescent.hpp
-	Solver.hpp
-	Solver.cpp
-	Newton.hpp
-	Newton.cpp
 	BoxConstraintSolver.cpp
 	BoxConstraintSolver.hpp
-	SparseNewton.hpp
-	SparseNewton.cpp
-	DenseNewton.hpp
-	DenseNewton.cpp
-	MMA.hpp
-	MMA.cpp
-	MMAAux.hpp
-	MMAAux.cpp
+	Solver.hpp
+	Solver.cpp
 	Problem.hpp
 	Problem.cpp
 )
@@ -31,4 +13,5 @@ target_sources(polysolve PRIVATE ${SOURCES})
 # ###############################################################################
 # Subfolders
 # ###############################################################################
+add_subdirectory(descent_strategies)
 add_subdirectory(line_search)