Commit ecd800a4 authored by David Bommes's avatar David Bommes
Browse files

improved infeasible start Newton

parent 2d3e97b5
Pipeline #17711 failed with stages
in 11 minutes and 30 seconds
......@@ -296,6 +296,7 @@ int NewtonSolver::solve_infeasible_start(NProblemInterface* _problem, const SMat
{
DEB_time_func_def;
converged_ = false;
mu_merit_ = 0.0; // reset
int max_regularization_iters = 10;
double max_regularization_value = 1e6;
......@@ -399,11 +400,15 @@ int NewtonSolver::solve_infeasible_start(NProblemInterface* _problem, const SMat
if (feasible_mode)
t = backtracking_line_search(_problem, x, g, dz, newton_decrement, fx, t_max);
else
t = backtracking_line_search_infeasible(_problem, A, b, x, nue, g, dz, fx, res_primal2, res_dual2, t_max,
30);
t = backtracking_line_search_infeasible_merit_l1(_problem, H, A, b, x, g, dz, fx, t_max,
100);
// old infeasible line search works well only for convex problems
// t = backtracking_line_search_infeasible(_problem, A, b, x, nue, g, dz, fx, res_primal2, res_dual2, t_max,
// 100);
// perform update
if (feasible_mode) // infeasible mode updates already during line search
// if (feasible_mode) // infeasible mode updates already during line search
{
// update primal variables
x += dz.head(n) * t;
......@@ -457,7 +462,11 @@ int NewtonSolver::solve_infeasible_start(NProblemInterface* _problem, const SMat
if (t < 0.1) // increase regularization if step length small
regularize_hessian = std::min(std::max(1e-3, 10.0 * regularize_hessian), max_regularization_value);
if (t == 1.0) // decrease regularization after full step
{
regularize_hessian = std::max(prescribed_hessian_regularization_relative_, 0.5 * regularize_hessian);
if(regularize_hessian < 1e-3) // set to lowest possible value when below 1e-3
regularize_hessian = prescribed_hessian_regularization_relative_;
}
}
++iter;
}
......@@ -717,8 +726,10 @@ double NewtonSolver::backtracking_line_search_infeasible(NProblemInterface* _pro
rd2 = (g_ls_ + _A.transpose() * nue_ls_).squaredNorm();
double r = rp2 + rd2;
std::cerr << " line search r=" << r << " vs r0=" << r0 << " r0-r=" << r0-r << std::endl;
// sufficient decrease in residual?
if (r < (1.0 - alpha_ls_*t) * r0)
if (r < std::pow(1.0 - alpha_ls_*t,2) * r0)
{
// succesfull line search
_res_primal2 = rp2;
......@@ -745,6 +756,70 @@ double NewtonSolver::backtracking_line_search_infeasible(NProblemInterface* _pro
//-----------------------------------------------------------------------------
double NewtonSolver::backtracking_line_search_infeasible_merit_l1(NProblemInterface* _problem, const SMatrixD& _H,
const SMatrixD& _A, const VectorD& _b,
const VectorD& _x, const VectorD& _g, VectorD& _dz, double& _fx,
const double _t_start, const int _max_ls_iters)
{
DEB_enter_func;
size_t n = _x.size();
size_t m = _b.size();
// update mu
double res_primal_1 = (_A*_x-_b).template lpNorm<1>();
double gdx = _g.transpose()*_dz.head(n);
double dxHdx = _dz.head(n).transpose()*_H*_dz.head(n);
double mu_new = 1.2*(gdx+0.5*std::max(0.0,dxHdx))/((1.0-rho_merit_)*res_primal_1);
mu_merit_ = std::max(mu_new, mu_merit_);
// current step size
double t = _t_start;
// merit function and directional derivative for t=0
// double merit_0 = _fx + mu_merit_*res_primal_1;
double merit_0 = _problem->eval_f(_x.data()) + mu_merit_*res_primal_1;
double D_merit_0 = gdx - mu_merit_*res_primal_1;
std::cerr << "mu_merit=" << mu_merit_ << std::endl;
std::cerr << "D_merit_0=" << D_merit_0 << std::endl;
double fx(0.0);
// backtracking (stable in case of NAN and with max iterations)
for(int i=0; i<_max_ls_iters; ++i)
{
// current update of x, nue and g
x_ls_ = _x + _dz.head(n)*t;
fx = _problem->eval_f(x_ls_.data());
// check if update is inside domain
if(std::isfinite(fx))
{
double merit_t = fx + mu_merit_*(_A*x_ls_-_b).template lpNorm<1>();
std::cerr << "t=" << t << ", merit(t)=" << merit_t << ", merit(0)+alpha*t*Dmerit(0)=" << merit_0 + alpha_ls_*t*D_merit_0 << std::endl;
// sufficient decrease in residual?
if (merit_t <= merit_0 + alpha_ls_*t*D_merit_0)
{
// succesfull line search
_fx = fx;
return t;
}
}
else std::cerr << "t=" << t << ", merit(t)=" << std::numeric_limits<double>::infinity() << std::endl;
// shrink with factor beta
t *= beta_ls_;
}
DEB_warning(1, "line search could not find a valid step");
return 0.0;
}
//-----------------------------------------------------------------------------
void NewtonSolver::analyze_pattern(SMatrixD& _KKT)
{
DEB_enter_func;
......
......@@ -66,6 +66,7 @@ public:
const double _beta_ls = 0.6)
: eps_(_eps), eps_ls_(_eps_line_search), max_iters_(_max_iters),
alpha_ls_(_alpha_ls), beta_ls_(_beta_ls),
rho_merit_(0.5), mu_merit_(0.0),
prescribed_constraint_regularization_absolute_(0.0), prescribed_hessian_regularization_relative_(0.0),
use_trust_region_regularization_(true),
solver_type_(LS_EigenLU), constant_hessian_structure_(false)
......@@ -163,6 +164,12 @@ protected:
VectorD& _x, VectorD& _nue, VectorD& _g, VectorD& _dz, double& _fx, double& _res_primal2, double& _res_dual2,
const double _t_start = 1.0, const int _max_ls_iters = 20);
double backtracking_line_search_infeasible_merit_l1(NProblemInterface* _problem, const SMatrixD& _H,
const SMatrixD& _A, const VectorD& _b,
const VectorD& _x, const VectorD& _g, VectorD& _dz, double& _fx,
const double _t_start=1.0, const int _max_ls_iters=20);
void analyze_pattern(SMatrixD& _KKT);
......@@ -207,6 +214,9 @@ private:
double alpha_ls_;
double beta_ls_;
double rho_merit_;
double mu_merit_;
double prescribed_constraint_regularization_absolute_;
double prescribed_hessian_regularization_relative_;
......
......@@ -64,8 +64,8 @@ public:
converged_ = false;
// number of unknowns
size_t n = _problem->n_unknowns();
DEB_line(2, "optimize via TruncatedNewtonPCG with " << n << " unknowns");
Eigen::Index n = _problem->n_unknowns();
DEB_line(2, "optimize via TruncatedNewtonPCG with " << (int)n << " unknowns");
// Newton parameters
const int max_iters = max_iters_;
......@@ -257,9 +257,9 @@ public:
converged_ = false;
// number of unknowns
size_t n = _problem->n_unknowns();
size_t m = _A.rows();
DEB_line(2, "optimize via TruncatedNewtonPCG with " << n << " unknowns and " << m << " linear constraints");
Eigen::Index n = _problem->n_unknowns();
Eigen::Index m = _A.rows();
DEB_line(2, "optimize via TruncatedNewtonPCG with " << (int)n << " unknowns and " << (int)m << " linear constraints");
// Newton parameters
const int max_iters = max_iters_;
......@@ -298,7 +298,7 @@ public:
double reg = 1e-8 * AAt.diagonal().sum() / double(m);
DEB_line(2,
"Warning: LDLT factorization failed --> regularize AAt (could lead to reduced accuracy), reg=" << reg);
for (size_t j = 0; j < AAt.rows(); ++j)
for (Eigen::Index j = 0; j < AAt.rows(); ++j)
AAt.coeffRef(j, j) += reg;
ldlt.compute(AAt);
......@@ -487,9 +487,9 @@ public:
}
// number of unknowns
size_t n = _problem->n_unknowns();
size_t m = _A.rows();
DEB_line(2, "optimize via TruncatedNewtonPCG with " << n << " unknowns and " << m << " linear constraints");
Eigen::Index n = _problem->n_unknowns();
Eigen::Index m = _A.rows();
DEB_line(2, "optimize via TruncatedNewtonPCG with " << (int)n << " unknowns and " << (int)m << " linear constraints");
// Newton parameters
const int max_iters = max_iters_;
......@@ -571,7 +571,7 @@ public:
double reg = 1e-8 * AAt.diagonal().sum() / double(m);
DEB_line(2,
"Warning: LDLT factorization failed --> regularize AAt (could lead to reduced accuracy), reg=" << reg);
for (size_t j = 0; j < AAt.rows(); ++j)
for (Eigen::Index j = 0; j < AAt.rows(); ++j)
AAt.coeffRef(j, j) += reg;
ldlt.compute(AAt);
......@@ -757,8 +757,8 @@ public:
converged_ = false;
// number of unknowns
size_t n = _problem->n_unknowns();
size_t m = _A.rows();
Eigen::Index n = _problem->n_unknowns();
Eigen::Index m = _A.rows();
std::cerr << "compute QR decomposition of A^T...\n";
SMatrixD At = _A.transpose();
......@@ -776,8 +776,8 @@ public:
Z = Q.middleCols(mq,nz); //ToDo: Is the performance of this operation ok?
std::cerr << "||A*Z|| = " << (_A*Z).norm() << std::endl;
DEB_line(2, "optimize via ReducedSystemNewtonPCG with " << n << " unknowns and " << m << " linear constraints (" << mq << " linearly independent)");
DEB_line(2, "----> reduced (unconstrained) problem has " << nz << " unknowns");
DEB_line(2, "optimize via ReducedSystemNewtonPCG with " << (int)n << " unknowns and " << (int)m << " linear constraints (" << (int)mq << " linearly independent)");
DEB_line(2, "----> reduced (unconstrained) problem has " << (int)nz << " unknowns");
// Newton parameters
const int max_iters = max_iters_;
......@@ -977,8 +977,8 @@ public:
converged_ = false;
// number of unknowns
size_t n = _problem->n_unknowns();
size_t m = _A.rows();
Eigen::Index n = _problem->n_unknowns();
Eigen::Index m = _A.rows();
std::cerr << "compute QR decomposition of A^T...\n";
StopWatch swq; swq.start();
......@@ -998,8 +998,8 @@ public:
// std::cerr << "||Z^T Z||^2 = " << std::pow((Z.transpose()*Z).norm(), 2) << std::endl;
// std::cerr << "||Z Z^T||^2 = " << std::pow((Z*Z.transpose()).norm(), 2) << std::endl;
DEB_line(2, "optimize via ReducedSystemNewtonPCGSimple with " << n << " unknowns and " << m << " linear constraints (" << mq << " linearly independent)");
DEB_line(2, "----> reduced (unconstrained) problem has " << nz << " unknowns");
DEB_line(2, "optimize via ReducedSystemNewtonPCGSimple with " << (int)n << " unknowns and " << (int)m << " linear constraints (" << (int)mq << " linearly independent)");
DEB_line(2, "----> reduced (unconstrained) problem has " << (int)nz << " unknowns");
// Newton parameters
const int max_iters = max_iters_;
......@@ -1151,10 +1151,10 @@ public:
converged_ = false;
// number of unknowns
size_t n = _problem->n_unknowns();
size_t m = _A.rows();
Eigen::Index n = _problem->n_unknowns();
Eigen::Index m = _A.rows();
DEB_line_if(!silent_, 2, "optimize via TruncatedNewtonProjectedNormalEquationsPCG with " << n << " unknowns and " << m << " linear constraints");
DEB_line_if(!silent_, 2, "optimize via TruncatedNewtonProjectedNormalEquationsPCG with " << (int)n << " unknowns and " << (int)m << " linear constraints");
// Newton parameters
const int max_iters = max_iters_;
......@@ -1236,7 +1236,7 @@ public:
double reg = 1e-8 * AWiAt.diagonal().sum() / double(m);
DEB_line_if(!silent_, 2,
"Warning: LDLT factorization failed --> regularize AAt (could lead to reduced accuracy), reg=" << reg);
for (size_t j = 0; j < AWiAt.rows(); ++j)
for (Eigen::Index j = 0; j < AWiAt.rows(); ++j)
AWiAt.coeffRef(j, j) += reg;
ldlt.compute(AWiAt);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment