Skip to content

Commit

Permalink
newton clenup
Browse files Browse the repository at this point in the history
  • Loading branch information
@teseoch
teseoch committed Nov 10, 2023
1 parent de5e9d1 commit 74612f4
Show file tree
Hide file tree
Showing 5 changed files with 52 additions and 36 deletions.
23 changes: 15 additions & 8 deletions non-linear-solver-spec.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,14 +102,21 @@
"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,
Expand All @@ -128,18 +135,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,
Expand Down
10 changes: 6 additions & 4 deletions src/polysolve/nonlinear/Solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -247,7 +247,8 @@ namespace polysolve::nonlinear

if (!ok || std::isnan(grad_norm) || (m_strategies[m_descent_strategy]->is_direction_descent() && grad_norm != 0 && delta_x.dot(grad) >= 0))
{
++m_descent_strategy;
if (!m_strategies[m_descent_strategy]->handle_error())
++m_descent_strategy;
if (m_descent_strategy >= m_strategies.size())
{
this->m_status = cppoptlib::Status::UserDefined;
Expand All @@ -267,7 +268,8 @@ namespace polysolve::nonlinear
const double delta_x_norm = delta_x.norm();
if (std::isnan(delta_x_norm))
{
++m_descent_strategy;
if (!m_strategies[m_descent_strategy]->handle_error())
++m_descent_strategy;

if (m_descent_strategy >= m_strategies.size())
{
Expand Down Expand Up @@ -297,8 +299,8 @@ namespace polysolve::nonlinear
if (std::isnan(rate))
{
assert(this->m_status == cppoptlib::Status::Continue);

++m_descent_strategy;
if (!m_strategies[m_descent_strategy]->handle_error())
++m_descent_strategy;
if (m_descent_strategy >= m_strategies.size())
{
this->m_status = cppoptlib::Status::UserDefined; // Line search failed on gradient descent, so quit!
Expand Down
1 change: 1 addition & 0 deletions src/polysolve/nonlinear/descent_strategies/DescentStrategy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ namespace polysolve::nonlinear
virtual void update_solver_info(json &solver_info, const double per_iteration) {}

virtual bool is_direction_descent() { return true; }
virtual bool handle_error() { return false; }

virtual bool compute_update_direction(
Problem &objFunc,
Expand Down
49 changes: 27 additions & 22 deletions src/polysolve/nonlinear/descent_strategies/Newton.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,22 @@ namespace polysolve::nonlinear
solver_params, linear_solver_params,
characteristic_length, logger));

// TODO disable regularization?
res.push_back(std::make_unique<RegularizedNewton>(
sparse,
solver_params, linear_solver_params,
characteristic_length, logger));

// TODO disable projection?
res.push_back(std::make_unique<ProjectedNewton>(
sparse,
solver_params, linear_solver_params,
characteristic_length, logger));
const double reg_weight_min = solver_params["Newton"]["reg_weight_min"];
if (reg_weight_min > 0)
res.push_back(std::make_unique<RegularizedNewton>(
sparse,
solver_params, linear_solver_params,
characteristic_length, logger));

const bool use_psd_projection = solver_params["Newton"]["use_psd_projection"];
if (!use_psd_projection)
res.push_back(std::make_unique<ProjectedNewton>(
sparse,
solver_params, linear_solver_params,
characteristic_length, logger));

if (res.empty())
log_and_throw_error(logger, "Newton needs to have at least one of force_psd_projection=false, reg_weight_min>0, or use_psd_projection=true");

return res;
}
Expand All @@ -42,6 +47,8 @@ namespace polysolve::nonlinear
: Superclass(solver_params, characteristic_length, logger),
is_sparse(sparse), m_characteristic_length(characteristic_length)
{
m_residual_tolerance = solver_params["Newton"]["residual_tolerance"];

linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger);
assert(linear_solver->is_dense() == !sparse);
}
Expand All @@ -67,7 +74,8 @@ namespace polysolve::nonlinear
reg_weight_min = solver_params["Newton"]["reg_weight_min"];
reg_weight_max = solver_params["Newton"]["reg_weight_max"];
reg_weight_inc = solver_params["Newton"]["reg_weight_inc"];
reg_weight_dec = solver_params["Newton"]["reg_weight_dec"];

reg_weight = reg_weight_min;
}

// =======================================================================
Expand All @@ -81,7 +89,7 @@ namespace polysolve::nonlinear
void RegularizedNewton::reset(const int ndof)
{
Superclass::reset(ndof);
reg_weight = 0;
reg_weight = reg_weight_min;
}

// =======================================================================
Expand Down Expand Up @@ -259,13 +267,18 @@ namespace polysolve::nonlinear
}
// =======================================================================

bool RegularizedNewton::handle_error()
{
reg_weight *= reg_weight_inc;
return reg_weight < reg_weight_max;
}
// =======================================================================

bool Newton::check_direction(
const double residual, const TVector &grad, const TVector &direction)
{
// gradient descent, check descent direction
if (std::isnan(residual) || residual > std::max(1e-8 * grad.norm(), 1e-5) * m_characteristic_length)
if (std::isnan(residual) || residual > m_residual_tolerance * m_characteristic_length)
{
m_logger.debug("[{}] large (or nan) linear solve residual {} (||∇f||={})",
name(), residual, grad.norm());
Expand All @@ -277,14 +290,6 @@ namespace polysolve::nonlinear
m_logger.trace("linear solve residual {}", residual);
}

// do this check here because we need to repeat the solve without resetting reg_weight
if (grad.norm() != 0 && grad.dot(direction) >= 0)
{
m_logger.debug("[{}] direction is not a descent direction (‖g‖={:g}; ‖Δx‖={:g}; Δx⋅g={:g}≥0)",
name(), grad.norm(), direction.norm(), direction.dot(grad));
return false;
}

return true;
}

Expand Down
5 changes: 3 additions & 2 deletions src/polysolve/nonlinear/descent_strategies/Newton.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ namespace polysolve::nonlinear
json internal_solver_info = json::array();

const double m_characteristic_length;
double m_residual_tolerance;
const bool is_sparse;

std::unique_ptr<polysolve::linear::Solver> linear_solver; ///< Linear solver used to solve the linear system
Expand Down Expand Up @@ -111,14 +112,14 @@ namespace polysolve::nonlinear
std::string name() const override { return internal_name() + "RegularizedNewton"; }

void reset(const int ndof) override;
bool handle_error() override;

private:
double reg_weight_min; // needs to be greater than zero
double reg_weight_max;
double reg_weight_inc;
double reg_weight_dec;

double reg_weight = 0; ///< Regularization Coefficients
double reg_weight; ///< Regularization Coefficients
protected:
void compute_hessian(Problem &objFunc,
const TVector &x,
Expand Down

0 comments on commit 74612f4

Please sign in to comment.