Rescaled tolerances

nonlinear-solver-spec.json

@@ -4,9 +4,13 @@
         "type": "object",
         "optional": [
             "solver",
-            "x_delta",
-            "grad_norm",
+            "x_delta_tol",
+            "grad_norm_tol",
+            "rel_grad_norm_tol",
+            "newton_decrement_tol",
+            "rel_x_delta_tol",
             "first_grad_norm_tol",
+            "norm_type",
             "max_iterations",
             "iterations_per_strategy",
             "line_search",
@@ -41,14 +45,35 @@
         "doc": "Nonlinear solver type"
     },
     {
-        "pointer": "/x_delta",
+        "pointer": "/x_delta_tol",
         "default": 0,
         "type": "float",
         "min": 0,
         "doc": "Stopping criterion: minimal change of the variables x for the iterations to continue. Computed as the L2 norm of x divide by the time step."
     },
     {
-        "pointer": "/grad_norm",
+        "pointer": "/rel_x_delta_tol",
+        "default": 0,
+        "type": "float",
+        "min": 0,
+        "doc": "Stopping criterion: minimal change of the variables x for the iterations to continue relative to first step in nonlinear solve."
+    },
+    {
+        "pointer": "/rel_grad_norm_tol",
+        "default": 0,
+        "type": "float",
+        "min": 0,
+        "doc": "Stopping criterion: minimal gradient for the iterations to continue relative to first step in nonlinear solve."
+    },
+    {
+        "pointer": "/newton_decrement_tol",
+        "default": 0,
+        "type": "float",
+        "min": 0,
+        "doc": "Stopping criterion: minimal change of energy (as estimated by Newton decrement) for the iterations to continue."
+    },
+    {
+        "pointer": "/grad_norm_tol",
         "default": 1e-08,
         "type": "float",
         "min": 0,
@@ -60,6 +85,17 @@
         "type": "float",
         "doc": "Minimal gradient norm for the iterations to not start, assume we already are at a minimum."
     },
+    {
+        "pointer": "/norm_type",
+        "default": "L2",
+        "type": "string",
+        "doc": "Norm to use when computing stopping criteria in nonlinear solve.",
+        "options": [
+            "Euclidean",
+            "L2",
+            "Linf"
+        ]
+    },
     {
         "pointer": "/max_iterations",
         "default": 500,
@@ -692,7 +728,7 @@
         "default": null,
         "type": "object",
         "optional": [
-            "f_delta",
+            "f_delta_tol",
             "f_delta_step_tol",
             "derivative_along_delta_x_tol",
             "apply_gradient_fd",
@@ -701,7 +737,7 @@
         "doc": "Nonlinear solver advanced options"
     },
     {
-        "pointer": "/advanced/f_delta",
+        "pointer": "/advanced/f_delta_tol",
         "default": 0,
         "min": 0,
         "type": "float",

src/polysolve/linear/AMGCL.hpp

@@ -76,6 +76,8 @@ namespace polysolve::linear
         // Solve the linear system Ax = b
         virtual void solve(const Ref<const VectorXd> b, Ref<VectorXd> x) override;
 
+        virtual void set_tolerance(const double tol) override {residual_error_ = tol;}
+
         // Name of the solver type (for debugging purposes)
         virtual std::string name() const override { return "AMGCL_Block" + std::to_string(BLOCK_SIZE); }
 

src/polysolve/linear/EigenSolver.hpp

@@ -64,6 +64,8 @@ namespace polysolve::linear
         // Constructor requires a solver name used for finding parameters in the json file passed to set_parameters
         EigenIterative(const std::string &name) { m_Name = name; }
 
+        virtual void set_tolerance(const double tol) override {m_Solver.setTolerance(tol);}
+
     public:
         // Set solver parameters
         virtual void set_parameters(const json &params) override;

src/polysolve/linear/HypreSolver.hpp

@@ -53,6 +53,8 @@ namespace polysolve::linear
         // Name of the solver type (for debugging purposes)
         virtual std::string name() const override { return "Hypre"; }
 
+        virtual void set_tolerance(const double tol) override {conv_tol_ = tol;}
+
     protected:
         int dimension_ = 1; // 1 = scalar (Laplace), 2 or 3 = vector (Elasticity)
         int max_iter_ = 1000;

src/polysolve/linear/Solver.hpp

@@ -113,6 +113,9 @@ namespace polysolve::linear
         /// If the problem is nullspace for multigrid solvers
         virtual void set_is_nullspace(const VectorXd &x) {}
 
+        /// Set solver tolerance 
+        virtual void set_tolerance(const double tol) {}
+
         ///
         /// @brief         { Solve the linear system Ax = b }
         ///

src/polysolve/nonlinear/BoxConstraintSolver.cpp

@@ -99,7 +99,8 @@ namespace polysolve::nonlinear
             direction);
     }
 
-    double BoxConstraintSolver::compute_grad_norm(const Eigen::VectorXd &x,
+    double BoxConstraintSolver::compute_grad_norm(const Problem &objFunc,
+                                                  const Eigen::VectorXd &x,
                                                   const Eigen::VectorXd &grad) const
     {
         auto min = get_lower_bound(x, false);

src/polysolve/nonlinear/BoxConstraintSolver.hpp

@@ -31,7 +31,7 @@ namespace polysolve::nonlinear
                             const double characteristic_length,
                             spdlog::logger &logger);
 
-        double compute_grad_norm(const Eigen::VectorXd &x, const Eigen::VectorXd &grad) const override;
+        double compute_grad_norm(const Problem &objFunc, const Eigen::VectorXd &x, const Eigen::VectorXd &grad) const override;
         Eigen::VectorXd get_lower_bound(const Eigen::VectorXd &x, bool consider_max_change = true) const;
         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;

src/polysolve/nonlinear/Criteria.cpp

@@ -9,7 +9,7 @@ namespace polysolve::nonlinear
 {
     bool is_converged_status(const Status s)
     {
-        return s == Status::XDeltaTolerance || s == Status::FDeltaTolerance || s == Status::GradNormTolerance;
+        return s == Status::XDeltaTolerance || s == Status::FDeltaTolerance || s == Status::GradNormTolerance || s == Status::NewtonDecrementTolerance || s == Status::RelGradNormTolerance || s == Status::RelXDeltaTolerance;
     }
 
     Criteria::Criteria()
@@ -26,6 +26,9 @@ namespace polysolve::nonlinear
         firstGradNorm = 0;
         fDeltaCount = 0;
         xDeltaDotGrad = 0;
+        relGradNorm = 0;
+        relXDelta = 0;
+        newtonDecrement = 0;
     }
 
     void Criteria::print(std::ostream &os) const
@@ -34,12 +37,16 @@ namespace polysolve::nonlinear
     }
     std::string Criteria::print_message() const {
         return fmt::format(
-            "iters={:d} {}={:g} {}={:g} {}={:g} {}={:g}",
+            "iters={:d} {}={:g} {}={:g} {}_rel={:g} {}={:g} {}_rel={:g} {}={:g} {}={:g}",
             iterations,
             log::delta("f"), fDelta,
             log::norm(log::grad("f")), gradNorm,
+            log::norm(log::grad("f")), relGradNorm,
             log::norm(log::delta("x")), xDelta,
-            log::delta("x") + log::dot() + log::grad("f(x)"), xDeltaDotGrad);
+            log::norm(log::delta("x")), relXDelta,
+            log::delta("x") + log::dot() + log::grad("f(x)"), xDeltaDotGrad,
+            "1/2" + log::delta("x") + "^TH" + log::delta("x"), newtonDecrement
+        );
     }
 
     Status checkConvergence(const Criteria &stop, const Criteria &current)
@@ -53,6 +60,18 @@ namespace polysolve::nonlinear
         {
             return Status::GradNormTolerance;
         }
+        if (stop.relXDelta > 0 && current.relXDelta < stop.relXDelta)
+        {
+            return Status::RelXDeltaTolerance;
+        }
+        if (stop.relGradNorm > 0 && current.relGradNorm < stop.relGradNorm)
+        {
+            return Status::RelGradNormTolerance;
+        }
+        if (stop.newtonDecrement > 0 && current.newtonDecrement < stop.newtonDecrement)
+        {
+            return Status::NewtonDecrementTolerance;
+        }
         if (stop.xDelta > 0 && current.xDelta < stop.xDelta)
         {
             return Status::XDeltaTolerance;
@@ -84,6 +103,12 @@ namespace polysolve::nonlinear
             return "Change in cost function value too small";
         case Status::GradNormTolerance:
             return "Gradient vector norm too small";
+        case Status::RelGradNormTolerance:
+            return "Relative gradient vector too small";
+        case Status::RelXDeltaTolerance:
+            return "Relative change in parameter vector too small";
+        case Status::NewtonDecrementTolerance:
+            return "Newton decrement too small";
         case Status::ObjectiveCustomStop:
             return "Objective function specified to stop";
         case Status::NanEncountered:

src/polysolve/nonlinear/Criteria.hpp

@@ -16,8 +16,11 @@ namespace polysolve::nonlinear
         // Success cases
         IterationLimit,      ///< The maximum number of iterations has been reached
         XDeltaTolerance,     ///< The change in the parameter vector is below the tolerance
+        RelXDeltaTolerance,
         FDeltaTolerance,     ///< The change in the cost function is below the tolerance
         GradNormTolerance,   ///< The norm of the gradient vector is below the tolerance
+        RelGradNormTolerance,
+        NewtonDecrementTolerance,
         ObjectiveCustomStop, ///< The objective function specified to stop
         // Failure cases
         NanEncountered,        ///< The objective function returned NaN
@@ -37,6 +40,9 @@ namespace polysolve::nonlinear
         double gradNorm;      ///< Minimum norm of gradient vector
         double firstGradNorm; ///< Initial norm of gradient vector
         double xDeltaDotGrad; ///< Dot product of parameter vector and gradient vector
+        double relXDelta;
+        double relGradNorm;
+        double newtonDecrement;
         unsigned fDeltaCount; ///< Number of steps where fDelta is satisfied
 
         Criteria();

src/polysolve/nonlinear/Problem.hpp

@@ -10,6 +10,14 @@
 
 namespace polysolve::nonlinear
 {
+
+    enum class NormType
+    {
+        EUCLIDEAN = 0,
+        L2 = 1,
+        Linf = 2
+    };
+
     /// @brief Class defining optimization problem to be solved. To be defined by user code
     class Problem
     {
@@ -103,6 +111,13 @@ namespace polysolve::nonlinear
         /// @return True if the solver should stop, false otherwise.
         virtual bool stop(const TVector &x) { return false; }
 
+        virtual double grad_norm_rescaling(const NormType norm_type) const {return 1;}
+        virtual double step_norm_rescaling(const NormType norm_type) const {return 1;}
+        virtual double energy_norm_rescaling(const NormType norm_type) const {return 1;}
+
+        virtual double grad_norm(const TVector &grad, const NormType norm_type) const {return grad.norm();}
+        virtual double step_norm(const TVector &x, const NormType norm_type) const {return x.norm();}
+
         /// --- Misc ----------------------------------------------------------
 
         /// @brief Sample the function along a direction.