Commit 49e4334c authored by Max Lyon's avatar Max Lyon

Merge branch 'NewtonWithNumericalRegularization' into 'master'

Newton with numerical regularization

See merge request !23
parents 416f4819 32df97e8
Pipeline #5342 passed with stages
in 6 minutes and 44 seconds
......@@ -8,27 +8,11 @@ clang-c++11:
tags:
- Linux
gcc-c++98:
script: "CI/ci-linux.sh gcc C++98"
tags:
- Linux
clang-c++98:
script: "CI/ci-linux.sh clang C++98"
tags:
- Linux
macos-c++11:
script: "CI/ci-mac.sh C++11"
tags:
- Apple
macos-c++98:
script: "CI/ci-mac.sh C++98"
tags:
- Apple
CoMISo-VS2013-Qt-5.5.1-x64:
variables:
BUILD_PLATFORM: "VS2013"
......
......@@ -121,6 +121,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
{
DEB_time_func_def;
const double KKT_res_eps = 1e-6;
const int max_KKT_regularization_iters = 40;
double regularize_constraints_limit = 1e-6;
const double max_allowed_constraint_violation2 = 1e-12;
// number of unknowns
size_t n = _problem->n_unknowns();
// number of constraints
......@@ -133,6 +138,8 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
VectorD x(n);
_problem->initial_x(x.data());
double initial_constraint_violation2 = (_A*x-_b).squaredNorm();
// storage of update vector dx and rhs of KKT system
VectorD dx(n+m), rhs(n+m), g(n);
rhs.setZero();
......@@ -144,21 +151,67 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
lu_solver_.isSymmetric(true);
// start with no regularization
double regularize(0.0);
double regularize_hessian(0.0);
double regularize_constraints(0.0);
int iter=0;
bool first_factorization = true;
while( iter < max_iters_)
{
// get Newton search direction by solving LSE
bool first_factorization = (iter ==0);
factorize(_problem, _A, _b, x, regularize, first_factorization);
double kkt_res2(0.0);
double constraint_res2(0.0);
int reg_iters(0);
do
{
// get Newton search direction by solving LSE
bool fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
first_factorization = false;
if(fact_ok)
{
// get rhs
_problem->eval_gradient(x.data(), g.data());
rhs.head(n) = -g;
rhs.tail(m) = _b - _A*x;
// get rhs
_problem->eval_gradient(x.data(), g.data());
rhs.head(n) = -g;
rhs.tail(m) = _b - _A*x;
// solve KKT system
solve_kkt_system(rhs, dx);
// check numerical stability of KKT system and regularize if necessary
kkt_res2 = (KKT_*dx-rhs).squaredNorm();
constraint_res2 = (_A*dx.head(n)-rhs.tail(m)).squaredNorm();
}
if(!fact_ok || kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2)
{
// alternatingly regularize hessian and constraints
if(reg_iters % 2 == 0 || regularize_constraints >= regularize_constraints_limit)
{
DEB_line(2, "Warning: numerical issues in KKT system with residual^2 " << kkt_res2 << " (" << constraint_res2 << ") -> regularize hessian");
if(regularize_hessian == 0.0)
regularize_hessian = 1e-6;
else
regularize_hessian *= 2.0;
}
else
{
DEB_line(2, "Warning: numerical issues in KKT system with residual^2 " << kkt_res2 << " (" << constraint_res2 << ") -> regularize constraints");
if(regularize_constraints == 0.0)
regularize_constraints = 1e-8;
else
regularize_constraints *= 2.0;
}
}
++reg_iters;
}
while( (kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
// solve KKT system
solve_kkt_system(rhs, dx);
// no valid step could be found?
if(kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2 || reg_iters >= max_KKT_regularization_iters)
{
DEB_line(2, "Warning: numerical issues in KKT system could not be resolved -> terminating NewtonSolver with current solution");
_problem->store_result(x.data());
return 0;
}
// get maximal reasonable step
double t_max = std::min(1.0,
......@@ -167,44 +220,31 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
// perform line-search
double newton_decrement(0.0);
double fx(0.0);
//double t = backtracking_line_search(_problem, x, g, dx, newton_decrement, fx, t_max);
const VectorD x0 = x;
const double fx0 = _problem->eval_f(x0.data());
double t = t_max;
bool x_ok = false;
for (int i = 0; i < 10; ++i < 3 ? t *= 0.95 : t /= 2)
{// clamp back t very mildly the first 3 steps and then aggressively (/ 2)
x = x0 + dx.head(n) * t;
fx = _problem->eval_f(x.data());
if (fx > fx0) // function value is larger
continue;
newton_decrement = fx0 - fx;
//check constraint violation
const VectorD r = _b - _A*x;
const double cnstr_vltn = r.norm();
if (cnstr_vltn > 1e-8)
continue;
x_ok = true;
break; // exit here to avoid changing t
}
double t = backtracking_line_search(_problem, x, g, dx, newton_decrement, fx, t_max);
if (!x_ok)
// perform update
x += dx.head(n)*t;
double constraint_violation2 = (_A*x-_b).squaredNorm();
if(constraint_violation2 > 2*initial_constraint_violation2 && constraint_violation2 > max_allowed_constraint_violation2)
{
DEB_line(2, "Line search failed, pulling back to the intial solution");
t = 0;
x = x0;
fx = fx0;
DEB_line(2, "Warning: numerical issues in KKT system leads to constraint violation -> recovery phase");
// restore old solution
x -= dx.head(n)*t;
regularize_constraints *= 0.5;
regularize_constraints_limit = regularize_constraints;
}
DEB_line(2, "iter: " << iter
<< ", f(x) = " << fx << ", t = " << t << " (tmax=" << t_max << ")"
<< (t < t_max ? " _clamped_" : "")
<< ", eps = [Newton decrement] = " << newton_decrement
<< ", constraint violation prior = " << rhs.tail(m).norm()
<< ", after = " << (_b - _A*x).norm());
<< ", after = " << (_b - _A*x).norm()
<< ", KKT residual^2 = " << kkt_res2);
// converged?
if(newton_decrement < eps_ || std::abs(t) < eps_ls_)
......@@ -224,8 +264,8 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
//-----------------------------------------------------------------------------
void NewtonSolver::factorize(NProblemInterface* _problem,
const SMatrixD& _A, const VectorD& _b, const VectorD& _x, double& _regularize,
bool NewtonSolver::factorize(NProblemInterface* _problem,
const SMatrixD& _A, const VectorD& _b, const VectorD& _x, double& _regularize_hessian, double& _regularize_constraints,
const bool _first_factorization)
{
DEB_enter_func;
......@@ -258,48 +298,31 @@ void NewtonSolver::factorize(NProblemInterface* _problem,
trips.push_back(Triplet(it.col(),it.row()+n,it.value()));
}
if(_regularize != 0.0)
// regularize constraints
// if(_regularize_constraints != 0.0)
for( int i=0; i<m; ++i)
trips.push_back(Triplet(n+i,n+i,_regularize));
trips.push_back(Triplet(n+i,n+i,_regularize_constraints));
// regularize Hessian
// if(_regularize_hessian != 0.0)
{
double ad(0.0);
for( int i=0; i<n; ++i)
ad += H.coeffRef(i,i);
ad *= _regularize_hessian/double(n);
for( int i=0; i<n; ++i)
trips.push_back(Triplet(i,i,ad));
}
// create KKT matrix
SMatrixD KKT(nf,nf);
KKT.setFromTriplets(trips.begin(), trips.end());
KKT_.resize(nf,nf);
KKT_.setFromTriplets(trips.begin(), trips.end());
// compute LU factorization
if(_first_factorization)
analyze_pattern(KKT);
analyze_pattern(KKT_);
bool success = numerical_factorization(KKT);
if(!success)
{
// add more regularization
if(_regularize == 0.0)
_regularize = 1e-8;
else
_regularize *= 2.0;
// print information
DEB_line(2, "Linear Solver reported problem while factoring KKT system: ");
DEB_line_if(solver_type_ == LS_EigenLU, 2, lu_solver_.lastErrorMessage());
// for( int i=0; i<m; ++i)
// trips.push_back(Triplet(n+i,n+i,_regularize));
// regularize full system
for( int i=0; i<n+m; ++i)
trips.push_back(Triplet(i,i,_regularize));
// create KKT matrix
KKT.setFromTriplets(trips.begin(), trips.end());
// compute LU factorization
analyze_pattern(KKT);
numerical_factorization(KKT);
// IGM_THROW_if(lu_solver_.info() != Eigen::Success, QGP_BOUNDED_DISTORTION_FAILURE);
}
return numerical_factorization(KKT_);
}
......
......@@ -89,9 +89,9 @@ public:
protected:
void factorize(NProblemInterface* _problem, const SMatrixD& _A,
const VectorD& _b, const VectorD& _x, double& _regularize,
const bool _first_factorization);
bool factorize(NProblemInterface* _problem, const SMatrixD& _A,
const VectorD& _b, const VectorD& _x, double& _regularize_hessian,
double& _regularize_constraints, const bool _first_factorization);
double backtracking_line_search(NProblemInterface* _problem, VectorD& _x,
VectorD& _g, VectorD& _dx, double& _newton_decrement, double& _fx,
......@@ -133,7 +133,6 @@ private:
double eps_;
double eps_ls_;
int max_iters_;
// double accelerate_;
double alpha_ls_;
double beta_ls_;
......@@ -141,6 +140,9 @@ private:
LinearSolver solver_type_;
// cache KKT Matrix
SMatrixD KKT_;
// Sparse LU decomposition
Eigen::SparseLU<SMatrixD> lu_solver_;
......
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