diff --git a/nonlinear-auxsolver-spec.json b/nonlinear-auxsolver-spec.json new file mode 100644 index 00000000..3bb2a4ed --- /dev/null +++ b/nonlinear-auxsolver-spec.json @@ -0,0 +1,741 @@ +[ + { + "pointer": "/", + "default": null, + "type": "object", + "optional": [ + "solver", + "x_delta", + "grad_norm", + "first_grad_norm_tol", + "max_iterations", + "iterations_per_strategy", + "line_search", + "allow_out_of_iterations", + "L-BFGS", + "L-BFGS-B", + "Newton", + "ADAM", + "StochasticADAM", + "StochasticGradientDescent", + "box_constraints", + "advanced" + ], + "doc": "Settings for nonlinear solver. Interior-loop linear solver settings are defined in the solver/linear section." + }, + { + "pointer": "/solver", + "default": "Newton", + "type": "string", + "options": [ + "Newton", + "DenseNewton", + "GradientDescent", + "ADAM", + "StochasticADAM", + "StochasticGradientDescent", + "L-BFGS", + "BFGS", + "L-BFGS-B", + "MMA" + ], + "doc": "Nonlinear solver type" + }, + { + "pointer": "/x_delta", + "default": 1e-5, + "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", + "default": 1e-05, + "type": "float", + "min": 0, + "doc": "Stopping criterion: Minimal gradient norm for the iterations to continue." + }, + { + "pointer": "/first_grad_norm_tol", + "default": 1e-10, + "type": "float", + "doc": "Minimal gradient norm for the iterations to not start, assume we already are at a minimum." + }, + { + "pointer": "/max_iterations", + "default": 50, + "type": "int", + "doc": "Maximum number of iterations for a nonlinear solve." + }, + { + "pointer": "/iterations_per_strategy", + "default": 5, + "type": "int", + "doc": "Number of iterations for every substrategy before reset." + }, + { + "pointer": "/iterations_per_strategy", + "type": "list", + "doc": "Number of iterations for every substrategy before reset." + }, + { + "pointer": "/iterations_per_strategy/*", + "default": 5, + "type": "int", + "doc": "Number of iterations for every substrategy before reset." + }, + { + "pointer": "/allow_out_of_iterations", + "default": true, + "type": "bool", + "doc": "If false, an exception will be thrown when the nonlinear solver reaches the maximum number of iterations." + }, + { + "pointer": "/L-BFGS", + "default": null, + "type": "object", + "optional": [ + "history_size" + ], + "doc": "Options for LBFGS." + }, + { + "pointer": "/L-BFGS/history_size", + "default": 6, + "type": "int", + "doc": "The number of corrections to approximate the inverse Hessian matrix." + }, + { + "pointer": "/L-BFGS-B", + "default": null, + "type": "object", + "optional": [ + "history_size" + ], + "doc": "Options for the boxed L-BFGS." + }, + { + "pointer": "/L-BFGS-B/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", + "force_psd_projection", + "use_psd_projection", + "use_psd_projection_in_regularized" + ], + "doc": "Options for Newton." + }, + { + "pointer": "/Newton/residual_tolerance", + "default": -1, + "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, + "type": "float", + "doc": "Minimum regulariztion weight." + }, + { + "pointer": "/Newton/reg_weight_max", + "default": 1e8, + "type": "float", + "doc": "Maximum regulariztion weight." + }, + { + "pointer": "/Newton/reg_weight_inc", + "default": 10, + "type": "float", + "doc": "Regulariztion weight increment." + }, + { + "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": "/Newton/use_psd_projection_in_regularized", + "default": true, + "type": "bool", + "doc": "Use PSD in regularized Newton." + }, + { + "pointer": "/ADAM", + "default": null, + "type": "object", + "optional": [ + "alpha", + "beta_1", + "beta_2", + "epsilon" + ], + "doc": "Options for ADAM." + }, + { + "pointer": "/ADAM/alpha", + "default": 0.001, + "type": "float", + "doc": "Parameter alpha for ADAM." + }, + { + "pointer": "/ADAM/beta_1", + "default": 0.9, + "type": "float", + "doc": "Parameter beta_1 for ADAM." + }, + { + "pointer": "/ADAM/beta_2", + "default": 0.999, + "type": "float", + "doc": "Parameter beta_2 for ADAM." + }, + { + "pointer": "/ADAM/epsilon", + "default": 1e-8, + "type": "float", + "doc": "Parameter epsilon for ADAM." + }, + { + "pointer": "/StochasticADAM", + "default": null, + "type": "object", + "optional": [ + "alpha", + "beta_1", + "beta_2", + "epsilon", + "erase_component_probability" + ], + "doc": "Options for ADAM." + }, + { + "pointer": "/StochasticADAM/alpha", + "default": 0.001, + "type": "float", + "doc": "Parameter alpha for ADAM." + }, + { + "pointer": "/StochasticADAM/beta_1", + "default": 0.9, + "type": "float", + "doc": "Parameter beta_1 for ADAM." + }, + { + "pointer": "/StochasticADAM/beta_2", + "default": 0.999, + "type": "float", + "doc": "Parameter beta_2 for ADAM." + }, + { + "pointer": "/StochasticADAM/epsilon", + "default": 1e-8, + "type": "float", + "doc": "Parameter epsilon for ADAM." + }, + { + "pointer": "/StochasticADAM/erase_component_probability", + "default": 0.3, + "type": "float", + "doc": "Probability of erasing a component on the gradient for ADAM." + }, + { + "pointer": "/StochasticGradientDescent", + "default": null, + "type": "object", + "optional": [ + "erase_component_probability" + ], + "doc": "Options for Stochastic Gradient Descent." + }, + { + "pointer": "/StochasticGradientDescent/erase_component_probability", + "default": 0.3, + "type": "float", + "doc": "Probability of erasing a component on the gradient for StochasticGradientDescent." + }, + { + "pointer": "/solver", + "type": "list", + "doc": "List of solvers for ballback. Eg, [{'type':'Newton'}, {'type':'L-BFGS'}, {'type':'GradientDescent'}] will solve using Newton, in case of failure will fallback to L-BFGS and eventually to GradientDescent" + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "Newton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance" + ], + "doc": "Options for Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "ProjectedNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance" + ], + "doc": "Options for projected Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "RegularizedNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance", + "reg_weight_min", + "reg_weight_max", + "reg_weight_inc" + ], + "doc": "Options for regularized Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "RegularizedProjectedNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance", + "reg_weight_min", + "reg_weight_max", + "reg_weight_inc" + ], + "doc": "Options for regularized projected Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "DenseNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance" + ], + "doc": "Options for Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "DenseProjectedNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance" + ], + "doc": "Options for projected Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "DenseRegularizedNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance", + "reg_weight_min", + "reg_weight_max", + "reg_weight_inc" + ], + "doc": "Options for regularized Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "DenseRegularizedProjectedNewton", + "required": [ + "type" + ], + "optional": [ + "residual_tolerance", + "reg_weight_min", + "reg_weight_max", + "reg_weight_inc" + ], + "doc": "Options for projected regularized Newton." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "GradientDescent", + "required": [ + "type" + ], + "doc": "Options for Gradient Descent." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "StochasticGradientDescent", + "required": [ + "type" + ], + "optional": [ + "erase_component_probability" + ], + "doc": "Options for Stochastic Gradient Descent." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "L-BFGS", + "required": [ + "type" + ], + "optional": [ + "history_size" + ], + "doc": "Options for L-BFGS." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "BFGS", + "required": [ + "type" + ], + "doc": "Options for BFGS." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "ADAM", + "required": [ + "type" + ], + "optional": [ + "alpha", + "beta_1", + "beta_2", + "epsilon" + ], + "doc": "Options for ADAM." + }, + { + "pointer": "/solver/*", + "type": "object", + "type_name": "StochasticADAM", + "required": [ + "type" + ], + "optional": [ + "alpha", + "beta_1", + "beta_2", + "epsilon", + "erase_component_probability" + ], + "doc": "Options for ADAM." + }, + { + "pointer": "/solver/*/type", + "type": "string", + "options": [ + "Newton", + "DenseNewton", + "ProjectedNewton", + "DenseProjectedNewton", + "RegularizedNewton", + "DenseRegularizedNewton", + "RegularizedProjectedNewton", + "DenseRegularizedProjectedNewton", + "GradientDescent", + "StochasticGradientDescent", + "ADAM", + "StochasticADAM", + "L-BFGS", + "BFGS" + ], + "doc": "Nonlinear solver type" + }, + { + "pointer": "/solver/*/residual_tolerance", + "default": -1, + "type": "float", + "doc": "Tolerance of the linear system residual. If residual is above, the direction is rejected." + }, + { + "pointer": "/solver/*/reg_weight_min", + "default": 1e-8, + "type": "float", + "doc": "Minimum regulariztion weight." + }, + { + "pointer": "/solver/*/reg_weight_max", + "default": 1e8, + "type": "float", + "doc": "Maximum regulariztion weight." + }, + { + "pointer": "/solver/*/reg_weight_inc", + "default": 10, + "type": "float", + "doc": "Regulariztion weight increment." + }, + { + "pointer": "/solver/*/erase_component_probability", + "default": 0.3, + "type": "float", + "doc": "Probability of erasing a component on the gradient for stochastic solvers." + }, + { + "pointer": "/solver/*/history_size", + "default": 6, + "type": "int", + "doc": "The number of corrections to approximate the inverse Hessian matrix." + }, + { + "pointer": "/solver/*/alpha", + "default": 0.001, + "type": "float", + "doc": "Parameter alpha for ADAM." + }, + { + "pointer": "/solver/*/beta_1", + "default": 0.9, + "type": "float", + "doc": "Parameter beta_1 for ADAM." + }, + { + "pointer": "/solver/*/beta_2", + "default": 0.999, + "type": "float", + "doc": "Parameter beta_2 for ADAM." + }, + { + "pointer": "/solver/*/epsilon", + "default": 1e-8, + "type": "float", + "doc": "Parameter epsilon for ADAM." + }, + { + "pointer": "/line_search", + "default": null, + "type": "object", + "optional": [ + "method", + "use_grad_norm_tol", + "min_step_size", + "max_step_size_iter", + "min_step_size_final", + "max_step_size_iter_final", + "default_init_step_size", + "step_ratio", + "Armijo", + "RobustArmijo" + ], + "doc": "Settings for line-search in the nonlinear solver" + }, + { + "pointer": "/line_search/method", + "default": "RobustArmijo", + "type": "string", + "options": [ + "Armijo", + "RobustArmijo", + "Backtracking", + "None" + ], + "doc": "Line-search type" + }, + { + "pointer": "/line_search/use_grad_norm_tol", + "default": 1e-6, + "type": "float", + "doc": "When the energy is smaller than use_grad_norm_tol, line-search uses norm of gradient instead of energy" + }, + { + "pointer": "/line_search/min_step_size", + "default": 1e-8, + "type": "float", + "doc": "Mimimum step size" + }, + { + "pointer": "/line_search/max_step_size_iter", + "default": 20, + "type": "int", + "doc": "Number of iterations" + }, + { + "pointer": "/line_search/min_step_size_final", + "default": 1e-10, + "type": "float", + "doc": "Mimimum step size for last descent strategy" + }, + { + "pointer": "/line_search/max_step_size_iter_final", + "default": 50, + "type": "int", + "doc": "Number of iterations for last descent strategy" + }, + { + "pointer": "/line_search/default_init_step_size", + "default": 1, + "type": "float", + "doc": "Initial step size" + }, + { + "pointer": "/line_search/step_ratio", + "default": 0.5, + "type": "float", + "doc": "Ratio used to decrease the step" + }, + { + "pointer": "/line_search/Armijo", + "default": null, + "type": "object", + "optional": [ + "c" + ], + "doc": "Options for Armijo." + }, + { + "pointer": "/line_search/Armijo/c", + "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", + "optional": [ + "bounds", + "max_change" + ], + "default": null + }, + { + "pointer": "/box_constraints/bounds", + "default": [], + "type": "list", + "doc": "Box constraints on optimization variables." + }, + { + "pointer": "/box_constraints/bounds/*", + "type": "list", + "doc": "Box constraint values on optimization variables." + }, + { + "pointer": "/box_constraints/bounds/*/*", + "type": "float", + "doc": "Box constraint values on optimization variables." + }, + { + "pointer": "/box_constraints/bounds/*", + "type": "float", + "doc": "Box constraint values on optimization variables." + }, + { + "pointer": "/box_constraints/max_change", + "default": -1, + "type": "float", + "doc": "Maximum change of optimization variables in one iteration, only for solvers with box constraints. Negative value to disable this constraint." + }, + { + "pointer": "/box_constraints/max_change", + "type": "list", + "doc": "Maximum change of optimization variables in one iteration, only for solvers with box constraints." + }, + { + "pointer": "/box_constraints/max_change/*", + "type": "float", + "doc": "Maximum change of every optimization variable in one iteration, only for solvers with box constraints." + }, + { + "pointer": "/advanced", + "default": null, + "type": "object", + "optional": [ + "f_delta", + "f_delta_step_tol", + "derivative_along_delta_x_tol", + "apply_gradient_fd", + "gradient_fd_eps" + ], + "doc": "Nonlinear solver advanced options" + }, + { + "pointer": "/advanced/f_delta", + "default": 0, + "min": 0, + "type": "float", + "doc": "Dangerous Option: Quit the optimization if the solver reduces the energy by less than f_delta for consecutive f_delta_step_tol steps." + }, + { + "pointer": "/advanced/f_delta_step_tol", + "default": 100, + "type": "int", + "doc": "Dangerous Option: Quit the optimization if the solver reduces the energy by less than f_delta for consecutive f_delta_step_tol steps." + }, + { + "pointer": "/advanced/derivative_along_delta_x_tol", + "default": 0, + "min": 0, + "type": "float", + "doc": "Quit the optimization if the directional derivative along the descent direction is smaller than this tolerance." + }, + { + "pointer": "/advanced/apply_gradient_fd", + "default": "None", + "type": "string", + "options": [ + "None", + "DirectionalDerivative", + "FullFiniteDiff" + ], + "doc": "Expensive Option: For every iteration of the nonlinear solver, run finite difference to verify gradient of energy." + }, + { + "pointer": "/advanced/gradient_fd_eps", + "default": 1e-7, + "type": "float", + "doc": "Expensive Option: Eps for finite difference to verify gradient of energy." + } +] \ No newline at end of file diff --git a/nonlinear-solver-spec.json b/nonlinear-solver-spec.json index fc1ff34b..9208f80b 100644 --- a/nonlinear-solver-spec.json +++ b/nonlinear-solver-spec.json @@ -1,4 +1,5 @@ -[{ +[ + { "pointer": "/", "default": null, "type": "object", @@ -136,7 +137,7 @@ }, { "pointer": "/Newton/residual_tolerance", - "default": 1e-5, + "default": 1e10, "type": "float", "doc": "Tolerance of the linear system residual. If residual is above, the direction is rejected." }, @@ -479,7 +480,7 @@ }, { "pointer": "/solver/*/residual_tolerance", - "default": 1e-5, + "default": 1e10, "type": "float", "doc": "Tolerance of the linear system residual. If residual is above, the direction is rejected." }, diff --git a/src/polysolve/nonlinear/Criteria.cpp b/src/polysolve/nonlinear/Criteria.cpp index 3dd6e112..23e771fc 100644 --- a/src/polysolve/nonlinear/Criteria.cpp +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -18,6 +18,8 @@ namespace polysolve::nonlinear void Criteria::reset() { + constexpr double NaN = std::numeric_limits::quiet_NaN(); + iterations = 0; xDelta = 0; fDelta = 0; @@ -25,13 +27,18 @@ namespace polysolve::nonlinear firstGradNorm = 0; fDeltaCount = 0; xDeltaDotGrad = 0; + + energy = NaN; + alpha = NaN; + step = NaN; } void Criteria::print(std::ostream &os) const { os << print_message(); } - std::string Criteria::print_message() const { + std::string Criteria::print_message() const + { return fmt::format( "iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} Δx⋅∇f(x)={:g}", iterations, fDelta, gradNorm, xDelta, xDeltaDotGrad); @@ -64,7 +71,8 @@ namespace polysolve::nonlinear return Status::Continue; } - std::string_view status_message(Status s) { + std::string_view status_message(Status s) + { switch (s) { case Status::NotStarted: diff --git a/src/polysolve/nonlinear/Criteria.hpp b/src/polysolve/nonlinear/Criteria.hpp index d4e770c5..61407145 100644 --- a/src/polysolve/nonlinear/Criteria.hpp +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -39,6 +39,9 @@ namespace polysolve::nonlinear double xDeltaDotGrad; ///< Dot product of parameter vector and gradient vector unsigned fDeltaCount; ///< Number of steps where fDelta is satisfied + double alpha; ///< LS alpha + double energy; ///< Energy at the current step + double step; ///< alpha * grad.norm() Criteria(); void reset(); @@ -50,7 +53,7 @@ namespace polysolve::nonlinear Status checkConvergence(const Criteria &stop, const Criteria ¤t); std::string_view status_message(Status s); - std::string criteria_message(const Criteria& s); + std::string criteria_message(const Criteria &s); std::ostream &operator<<(std::ostream &os, const Status &s); diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index 1b98faa1..799d4d91 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -299,6 +299,7 @@ namespace polysolve::nonlinear POLYSOLVE_SCOPED_STOPWATCH("compute objective function", obj_fun_time, m_logger); energy = objFunc(x); } + m_current.energy = energy; if (!std::isfinite(energy)) { @@ -413,6 +414,7 @@ namespace polysolve::nonlinear POLYSOLVE_SCOPED_STOPWATCH("line search", line_search_time, m_logger); rate = m_line_search->line_search(x, delta_x, objFunc); } + m_current.alpha = rate; if (std::isnan(rate)) { @@ -466,6 +468,7 @@ namespace polysolve::nonlinear // Post update // ----------- const double step = (rate * delta_x).norm(); + m_current.step = step; // m_logger.debug("[{}][{}] rate={:g} ‖step‖={:g}", // descent_strategy_name(), m_line_search->name(), rate, step); @@ -481,6 +484,12 @@ namespace polysolve::nonlinear m_current.fDeltaCount = (m_current.fDelta < m_stop.fDelta) ? (m_current.fDeltaCount + 1) : 0; + if (m_iteration_callback && m_iteration_callback(m_current)) + { + m_status = Status::ObjectiveCustomStop; + m_logger.debug("[{}][{}] Iteration callback decided to stop", descent_strategy_name(), m_line_search->name()); + } + m_logger.debug( "[{}][{}] {} (stopping criteria: {})", descent_strategy_name(), m_line_search->name(), m_current.print_message(), m_stop.print_message()); diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index 211ea981..9f291d04 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -79,15 +79,16 @@ namespace polysolve::nonlinear void set_line_search(const json ¶ms); const json &info() const { return solver_info; } + void set_iteration_callback(std::function callback) { m_iteration_callback = callback; } + /// @brief If true the solver will not throw an error if the maximum number of iterations is reached bool allow_out_of_iterations = false; - /// @brief Get the line search object const std::shared_ptr &line_search() const { return m_line_search; }; protected: - /// @brief Compute direction in which the argument should be updated + /// @brief Compute direction in which the argument should be updated /// @param objFunc Problem to be minimized /// @param x Current input (n x 1) /// @param grad Gradient at current step (n x 1) @@ -106,8 +107,8 @@ namespace polysolve::nonlinear Criteria m_stop; /// @brief Current criteria - Criteria m_current; + Criteria m_current; /// @brief Current status Status m_status = Status::NotStarted; @@ -139,7 +140,7 @@ namespace polysolve::nonlinear // ==================================================================== // Solver state // ==================================================================== - + /// @brief Reset the solver at the start of a minimization /// @param ndof number of degrees of freedom void reset(const int ndof); @@ -151,12 +152,14 @@ namespace polysolve::nonlinear std::vector m_iter_per_strategy; + std::function m_iteration_callback = nullptr; + // ==================================================================== // Solver info // ==================================================================== /// @brief Update solver info JSON object - /// @param energy + /// @param energy void update_solver_info(const double energy); /// @brief Reset timing members to 0 diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.cpp b/src/polysolve/nonlinear/descent_strategies/Newton.cpp index 6ddea33b..64c1c0b6 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.cpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.cpp @@ -68,9 +68,6 @@ namespace polysolve::nonlinear linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger); if (linear_solver->is_dense() == sparse) log_and_throw_error(logger, "Newton linear solver must be {}, instead got {}", sparse ? "sparse" : "dense", linear_solver->name()); - - if (residual_tolerance <= 0) - log_and_throw_error(logger, "Newton residual_tolerance must be > 0, instead got {}", residual_tolerance); } Newton::Newton( @@ -145,7 +142,7 @@ namespace polysolve::nonlinear is_sparse ? solve_sparse_linear_system(objFunc, x, grad, direction) : solve_dense_linear_system(objFunc, x, grad, direction); - if (std::isnan(residual) || residual > residual_tolerance * characteristic_length) + if (std::isnan(residual) || (residual_tolerance > 0 && residual > residual_tolerance * characteristic_length)) { m_logger.debug("[{}] large (or nan) linear solve residual {}>{} (‖∇f‖={})", name(), residual, residual_tolerance * characteristic_length, grad.norm());