prox-solver.hpp

Go to the documentation of this file.

#pragma once
 
#include "proxsuite-nlp/fwd.hpp"
#include "proxsuite-nlp/problem-base.hpp"
#include "proxsuite-nlp/pdal.hpp"
#include "proxsuite-nlp/workspace.hpp"
#include "proxsuite-nlp/results.hpp"
#include "proxsuite-nlp/helpers-base.hpp"
#include "proxsuite-nlp/logger.hpp"
#include "proxsuite-nlp/bcl-params.hpp"
#include "proxsuite-nlp/bfgs-strategy.hpp"
 
#include <boost/mpl/bool.hpp>
 
#include "proxsuite-nlp/modelling/costs/squared-distance.hpp"
 
#include "proxsuite-nlp/linesearch-base.hpp"
 
namespace proxsuite {
namespace nlp {
 
enum class MultiplierUpdateMode { NEWTON, PRIMAL, PRIMAL_DUAL };
 
enum class HessianApprox {
  EXACT,
  GAUSS_NEWTON,
  BFGS,
  IDENTITY
};
 
enum InertiaFlag { INERTIA_OK = 0, INERTIA_BAD = 1, INERTIA_HAS_ZEROS = 2 };
 
enum KktSystem { KKT_CLASSIC, KKT_PRIMAL_DUAL };

template <typename _Scalar> class ProxNLPSolverTpl {
public:
  using Scalar = _Scalar;
  PROXSUITE_NLP_DYNAMIC_TYPEDEFS(Scalar);
  using Problem = ProblemTpl<Scalar>;
  using Workspace = WorkspaceTpl<Scalar>;
  using Results = ResultsTpl<Scalar>;
  using Manifold = ManifoldAbstractTpl<Scalar>;
  using LinesearchOptions = typename Linesearch<Scalar>::Options;
  using CallbackPtr = shared_ptr<helpers::base_callback<Scalar>>;
  using ConstraintSet = ConstraintSetTpl<Scalar>;
  using ConstraintObject = ConstraintObjectTpl<Scalar>;
  using BFGS = BFGSStrategy<Scalar, BFGSType::Hessian>;
 
protected:
  Problem *problem_;
 
public:
  ALMeritFunctionTpl<Scalar> merit_fun;
  QuadraticDistanceCostTpl<Scalar> prox_penalty;

  VerboseLevel verbose = QUIET;
  HessianApprox hess_approx = HessianApprox::GAUSS_NEWTON;
  LinesearchStrategy ls_strat = LinesearchStrategy::ARMIJO;
  MultiplierUpdateMode mul_update_mode = MultiplierUpdateMode::NEWTON;
 
  mutable BFGS bfgs_strategy_;

  std::size_t max_refinement_steps_ = 5;
  Scalar kkt_tolerance_ = 1e-13;
  LDLTChoice ldlt_choice_;
  KktSystem kkt_system_ = KKT_CLASSIC;

 
  Scalar inner_tol0 = 1.;
  Scalar prim_tol0 = 1.;
  Scalar inner_tol_ = inner_tol0;
  Scalar prim_tol_ = prim_tol0;
  Scalar rho_init_; //< Initial primal proximal penalty parameter.
  Scalar mu_init_;  //< Initial penalty parameter.
private:
  Scalar rho_ = rho_init_;   //< Primal proximal penalty parameter.
  Scalar mu_ = mu_init_;     //< Penalty parameter.
  Scalar mu_inv_ = 1. / mu_; //< Inverse penalty parameter.
public:
  Scalar inner_tol_min = 1e-9; //< Lower safeguard for the subproblem tolerance.
  Scalar mu_upper_ = 1.;       //< Upper safeguard for the penalty parameter.
  Scalar mu_lower_ = 1e-9;     //< Lower safeguard for the penalty parameter.
  Scalar pdal_beta_ = 0.5;     //< primal-dual weight for the dual variables.

  BCLParamsTpl<Scalar> bcl_params;

  LinesearchOptions ls_options;

  Scalar target_tol;

  BaseLogger logger{};

 
  Scalar del_inc_k = 8.;      //< Inertia corrector increase factor
  Scalar del_inc_big = 100.;  //< Inertia corrector increase factor (big)
  Scalar del_dec_k = 1. / 3.; //< Inertia corrector decrease factor
  Scalar DELTA_MIN = 1e-14;   //< Minimum nonzero regularization strength.
  Scalar DELTA_MAX = 1e6;     //< Maximum regularization strength.
  Scalar DELTA_NONZERO_INIT = 1e-4;
  Scalar DELTA_INIT = 0.;
 
  std::size_t max_iters = 100;     //< Solver maximum number of iterations.
  std::size_t max_al_iters = 1000; //< Maximum outer (AL) iterations.

  std::vector<CallbackPtr> callbacks_;
 
  unique_ptr<Workspace> workspace_;
  unique_ptr<Results> results_;
 
  ProxNLPSolverTpl(Problem &prob, const Scalar tol = 1e-6,
                   const Scalar mu_eq_init = 1e-2, const Scalar rho_init = 0.,
                   const VerboseLevel verbose = QUIET,
                   const Scalar mu_lower = 1e-9, const Scalar prim_alpha = 0.1,
                   const Scalar prim_beta = 0.9, const Scalar dual_alpha = 1.,
                   const Scalar dual_beta = 1.,
                   LDLTChoice ldlt_blocked = LDLTChoice::BUNCHKAUFMAN,
                   const LinesearchOptions ls_options = LinesearchOptions());
 
  const Problem &problem() const { return *problem_; }
 
  const Manifold &manifold() const { return *problem_->manifold_; }
 
  void setup() {
    workspace_ = std::make_unique<Workspace>(*problem_, ldlt_choice_);
    results_ = std::make_unique<Results>(*problem_);
 
    if (hess_approx == HessianApprox::BFGS) {
      // TODO: work on implementation of BFGS on manifolds
      if (problem_->ndx() == problem_->nx()) {
        BFGS valid_bfgs_strategy_ = BFGS(problem_->ndx());
        bfgs_strategy_ = std::move(valid_bfgs_strategy_);
        //  set workspace hessian to identity matrix (init to zero in workspace)
        workspace_->objective_hessian.setIdentity();
      } else {
        throw std::runtime_error("BFGS for hessian approximation currently "
                                 "only works for Euclidean manifolds.");
      }
    }
  }

  ConvergenceFlag solve(const ConstVectorRef &x0,
                        const std::vector<VectorRef> &lams0);
 
  PROXSUITE_NLP_DEPRECATED
  const Workspace &getWorkspace() const { return *workspace_; }
  PROXSUITE_NLP_DEPRECATED
  const Results &getResults() const { return *results_; }

  ConvergenceFlag solve(const ConstVectorRef &x0,
                        const ConstVectorRef &lams0 = VectorXs(0));
 
  void innerLoop(Workspace &workspace, Results &results);
 
  void assembleKktMatrix(Workspace &workspace);

  PROXSUITE_NLP_INLINE bool iterativeRefinement(Workspace &workspace) const;

  inline void updatePenalty();

  void setDualPenalty(const Scalar beta) { pdal_beta_ = beta; }

  void setPenalty(const Scalar &new_mu) noexcept;

  void setProxParameter(const Scalar &new_rho) noexcept;

  inline void registerCallback(const CallbackPtr &cb) noexcept {
    callbacks_.push_back(cb);
  }

  void clearCallbacks() noexcept { callbacks_.clear(); }

  void updateToleranceFailure() noexcept;

  void updateToleranceSuccess() noexcept;
 
  inline void tolerancePostUpdate() noexcept;

  void acceptMultipliers(Results &results, Workspace &workspace) const;

  void computeMultipliers(const ConstVectorRef &inner_lams_data,
                          Workspace &workspace) const;

  void computeProblemDerivatives(const ConstVectorRef &x, Workspace &workspace,
                                 boost::mpl::false_) const;
  void computeProblemDerivatives(const ConstVectorRef &x, Workspace &workspace,
                                 boost::mpl::true_) const;

  void computePrimalResiduals(Workspace &workspace, Results &results) const;

  void tryStep(Workspace &workspace, const Results &results, Scalar alpha);
 
  void invokeCallbacks(Workspace &workspace, Results &results) {
    for (auto cb : callbacks_) {
      cb->call(workspace, results);
    }
  }
};

inline InertiaFlag checkInertia(const int ndx, const int nc,
                                const Eigen::VectorXi &signature);
 
} // namespace nlp
} // namespace proxsuite
 
#include "proxsuite-nlp/prox-solver.hxx"
 
#ifdef PROXSUITE_NLP_ENABLE_TEMPLATE_INSTANTIATION
#include "proxsuite-nlp/prox-solver.txx"
#endif