ROL_PrimalDualInteriorPointStep.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef ROL_PRIMALDUALINTERIORPOINTSTEP_H
#define ROL_PRIMALDUALINTERIORPOINTSTEP_H

#include "ROL_VectorsNorms.hpp"
#include "ROL_PartitionedVector.hpp"
#include "ROL_InteriorPointPenalty.hpp"
#include "ROL_PrimalDualInteriorPointResidual.hpp"
#include "ROL_Krylov.hpp"


namespace ROL {
namespace InteriorPoint {
template <class Real> 
class PrimalDualInteriorPointStep : public Step<Real> {

  typedef Vector<Real>                                V;
  typedef PartitionedVector<Real>                     PV;
  typedef Objective<Real>                             OBJ;
  typedef BoundConstraint<Real>                       BND;
  typedef Krylov<Real>                                KRYLOV;
  typedef LinearOperator<Real>                        LINOP;
  typedef LinearOperatorFromConstraint<Real>          LOPEC;
  typedef Constraint<Real>                            EQCON;
  typedef StepState<Real>                             STATE;
  typedef InteriorPointPenalty<Real>                  PENALTY;
  typedef PrimalDualInteriorPointResidual<Real>       RESIDUAL;

private:

  ROL::Ptr<KRYLOV> krylov_;      // Krylov solver for the Primal Dual system 
  ROL::Ptr<LINOP>  precond_;     // Preconditioner for the Primal Dual system

  ROL::Ptr<BND> pbnd_;           // bound constraint for projecting x sufficiently away from the given bounds

  ROL::Ptr<V> x_;                // Optimization vector
  ROL::Ptr<V> g_;                // Gradient of the Lagrangian
  ROL::Ptr<V> l_;                // Lagrange multiplier

  ROL::Ptr<V> xl_;               // Lower bound vector
  ROL::Ptr<V> xu_;               // Upper bound vector

  ROL::Ptr<V> zl_;               // Lagrange multiplier for lower bound
  ROL::Ptr<V> zu_;               // Lagrange multiplier for upper bound

  ROL::Ptr<V> xscratch_;         // Scratch vector (size of x)
  ROL::Ptr<V> lscratch_;         // Scratch vector (size of l)

  ROL::Ptr<V> zlscratch_;        // Scratch vector (size of x)
  ROL::Ptr<V> zuscratch_;        // Scratch vector (size of x)

  ROL::Ptr<V> maskL_;            // Elements are 1 when xl>-INF, zero for xl = -INF
  ROL::Ptr<V> maskU_;            // Elements are 1 when xu< INF, zero for xu =  INF

  int iterKrylov_;
  int flagKrylov_;       

  bool symmetrize_;        // Symmetrize the Primal Dual system if true

  Elementwise::Multiply<Real> mult_; 

  // Parameters used by the Primal-Dual Interior Point Method

  Real mu_;                // Barrier penalty parameter
 
  Real eps_tol_;           // Error function tolerance   
  Real tau_;               // Current fraction-to-the-boundary parameter
  Real tau_min_;           // Minimum fraction-to-the-boundary parameter
  Real kappa_eps_;
  Real kappa_mu_;
  Real kappa1_;            // Feasibility projection parameter
  Real kappa2_;            // Feasibility projection parameter
  Real kappa_eps_;         // 
  Real lambda_max_;        //  multiplier maximum value
  Real theta_mu_;
  Real gamma_theta_;
  Real gamma_phi_
  Real sd_;                // Lagragian gradient scaling parameter
  Real scl_;               // Lower bound complementarity scaling parameter
  Real scu_;               // Upper bound complementarity scaling parameter
  Real smax_;              // Maximum scaling parameter

  Real diml_;              // Dimension of constaint
  Real dimx_;              // Dimension of optimization vector



  void updateState( const V& x, const V &l, OBJ &obj, 
                    EQCON &con, BND &bnd, ALGO &algo_state ) {
  
    Real tol = std::sqrt(ROL_EPSILON<Real>());
    ROL::Ptr<STATE> state = Step<Real>::getState();

    obj.update(x,true,algo_state.iter);
    con.update(x,true,algo_state.iter);

    algo_state.value = obj.value(tol);
    con.value(*(state->constraintVec),x,tol);

    obj.gradient(*(state->gradientVec),x,tol);   
    con.applyAdjointJacobian(*g_,l,x,tol);

    state->gradientVec->plus(*g_); // \f$ \nabla f(x)-\nabla c(x)\lambda \f$

    state->gradientVec->axpy(-1.0,*zl_); 
    state->gradientVec->axpy(-1.0,*zu_);

    // Scaled Lagrangian gradient sup norm
    algo_state.gnorm = normLinf(*(state->gradientVec))/sd_;

    // Constraint sup norm 
    algo_state.cnorm = normLinf(*(state->constraintVec));

    Elementwise::Multiply<Real> mult;

    Real lowerViolation;
    Real upperViolation;
     
    // Deviation from complementarity
    xscratch_->set(x);
    xscratch_->applyBinary(mult,*zl_);
    
    exactLowerViolation = normLinf(*xscratch_)/scl_;
 
    xscratch_->set(x);
    xscratch_->applyBunary(mult,*zu_);
    
    exactUpperBound = normLinf(*xscratch_)/scu_;

    // Measure ||xz||  
    algo_state.aggregateModelError = std::max(exactLowerViolation,
                                              exactUpperViolation);

  }

  /* When the constraint Jacobians are ill-conditioned, we can compute 
     multiplier vectors with very large norms, making it difficult to 
     satisfy the primal-dual equations to a small tolerance. These 
     parameters allow us to rescale the constraint qualification component
     of the convergence criteria 
  */
  void updateScalingParameters(void) {

    Real nl  = normL1(*l_);
    Real nzl = normL1(*zl_);
    Real nzu = normL1(*zu_);

    sd_  = (nl+nzl+nzu)/(diml_+2*dimx);
    sd_  = std::max(smax_,sd_)

    sd_ /= smax_;

    scl_  = nzl/dimx_;
    scl_  = std::max(smax_,scl_);
    scl_ /= smax_;

    scu_  = nzu/dimx_;
    scu_  = std::max(smax_,scu_);
    scu_ /= smax_;
  }

public:

  using Step<Real>::initialize;
  using Step<Real>::compute;
  using Step<Real>::update;

  PrimalDualInteriorPointStep( ROL::ParameterList &parlist, 
                               const ROL::Ptr<Krylov<Real> > &krylov = ROL::nullPtr,
                               const ROL::Ptr<LinearOperator<Real> > &precond = ROL::nullPtr ) : 
    Step<Real>(), krylov_(krylov), precond_(precond), iterKrylov_(0), flagKrylov_(0) {

    typedef ROL::ParameterList PL;

    PL &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");

    kappa1_     = iplist.get("Bound Perturbation Coefficient 1",        1.e-2);
    kappa2_     = iplist.get("Bound Perturbation Coefficient 2",        1.e-2);
    lambda_max_ = iplist.get(" Multiplier Maximum Value",       1.e3 ); 
    smax_       = iplist.get("Maximum Scaling Parameter",               1.e2 );
    tau_min_    = iplist.get("Minimum Fraction-to-Boundary Parameter",  0.99 );
    kappa_mu_   = iplist.get("Multiplicative Penalty Reduction Factor", 0.2  );
    theta_mu_   = iplist.get("Penalty Update Power",                    1.5  );
    eps_tol_    = iplist.get("Error Tolerance",                         1.e-8);
    symmetrize_ = iplist.get("Symmetrize Primal Dual System",           true );

    PL &filter  = iplist.sublist("Filter Parameters");

    if(krylov_ == ROL::nullPtr) {
      krylov_ = KrylovFactory<Real>(parlist);
    }

    if( precond_ == ROL::nullPtr) {
      class IdentityOperator : public LINOP {
      public: 
        apply( V& Hv, const V &v, Real tol ) const {
          Hv.set(v);
        }
      }; // class IdentityOperator

      precond_ = ROL::makePtr<IdentityOperator<Real>>();
    }

  } 



  ~PrimalDualInteriorPointStep() {


  void initialize( Vector<Real> &x, const Vector<Real> &g, Vector<Real> &l, const Vector<Real> &c,
                   Objective<Real> &obj, Constraint<Real> &con, BoundConstraint<Real> &bnd,
                   AlgorithmState<Real> &algo_state ) {

     
    using Elementwise::ValueSet;   

    ROL::Ptr<PENALTY> &ipPen = dynamic_cast<PENALTY&>(obj);

    // Initialize step state
    ROL::Ptr<STATE> state = Step<Real>::getState();    
    state->descentVec   = x.clone();
    state->gradientVec  = g.clone();
    state->constaintVec = c.clone();

    diml_ = l.dimension();
    dimx_ = x.dimension();

    Real one(1.0); 
    Real zero(0.0);
    Real tol = std::sqrt(ROL_EPSILON<Real>());

    x_ = x.clone();
    g_ = g.clone();
    l_ = l.clone();
    c_ = c.clone();

    xscratch_ = x.clone(); 
    lscratch_ = l.clone(); 

    zlscratch_ = x.clone(); 
    zuscratch_ = x.clone(); 

    // Multipliers for lower and upper bounds 
    zl_ = x.clone();
    zu_ = x.clone();

    /*******************************************************************************************/
    /* Identify (implicitly) the index sets of upper and lower bounds by creating mask vectors */
    /*******************************************************************************************/
  
    xl_ = bnd.getLowerBound();
    xu_ = bnd.getUpperBound();

    maskl_ = ipPen.getLowerMask();
    masku_ = ipPen.getUpperMask();
    
    // Initialize bound constraint multipliers to 1 one where the corresponding bounds are finite
    zl_->set(maskl_);
    zu_->set(masku_);

    /*******************************************************************************************/
    /* Create a new bound constraint with perturbed bounds                                     */
    /*******************************************************************************************/
    
    ROL::Ptr<V> xdiff = xu_->clone();
    xdiff->set(*xu_); 
    xdiff->axpy(-1.0,*xl_);
    xdiff->scale(kappa2_);    

    class Max1X : public Elementwise::UnaryFunction<Real> {
    public:
      Real apply( const Real &x ) const {
        return std::max(1.0,x);
      }
    };

    Max1X                            max1x;
    Elementwise::AbsoluteValue<Real> absval;
    Elementwise::Min                 min;

    // Lower perturbation vector
    ROL::Ptr<V> pl = xl_->clone();
    pl->applyUnary(absval);
    pl->applyUnary(max1x);               // pl_i = max(1,|xl_i|)
    pl->scale(kappa1_);
    pl->applyBinary(min,xdiff);          // pl_i = min(kappa1*max(1,|xl_i|),kappa2*(xu_i-xl_i))

    // Upper perturbation vector
    ROL::Ptr<V> pu = xu_->clone();
    pu->applyUnary(absval);
    pu->applyUnary(max1x);              // pu_i = max(1,|xu_i|)
    pu->scale(kappa_1); 
    pu->applyBinary(min,xdiff);         // pu_u = min(kappa1*max(1,|xu_i|,kappa2*(xu_i-xl_i)))

    // Modified lower and upper bounds so that x in [xl+pl,xu-pu] using the above perturbation vectors
    pl->plus(*xl_);
    pu->scale(-1.0);
    pu->plus(*xu_);

    pbnd_ = ROL::makePtr<BoundConstraint<Real>>(pl,pu)

    // Project the initial guess onto the perturbed bounds
    pbnd_->project(x);


    /*******************************************************************************************/
    /* Solve least-squares problem for initial equality multiplier                             */
    /*******************************************************************************************/
   
    // \f$-(\nabla f-z_l+z_u) \f$
    g_->set(*zl_);
    g_->axpy(-1.0,g);
    g_->axpy(-1.0,*zu_)

    // We only need the updated multiplier
    lscratch_->zero();
    con.solveAugmentedSystem(*xscratch_,*l_,*g_,*lscratch_,x,tol);

    // If the multiplier supremum is too large, zero the vector as this appears to avoid poor
    // initial guesses of the multiplier in the case where the constraint Jacobian is 
    // ill-conditioned 
    if( normInf(l_) > lambda_max_ ) {
      l_->zero();
    }
 
    // Initialize the algorithm state
    algo_state.nfval = 0;
    algo_state.ncval = 0;
    algo_state.ngrad = 0;

    updateState(x,l,obj,con,bnd,algo_state);

  } // initialize()


  void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l, 
                Objective<Real> &obj, Constraint<Real> &con, 
                BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {


      

    Elementwise::Fill<Real>     minus_mu(-mu_); 
    Elementwise::Divide<Real>   div;
    Elementwise::Multiply<Real> mult;

    ROL::Ptr<STATE> state = Step<Real>::getState();    

    ROL::Ptr<OBJ>   obj_ptr = ROL::makePtrFromRef(obj);
    ROL::Ptr<EQCON> con_ptr = ROL::makePtrFromRef(con); 
    ROL::Ptr<BND>   bnd_ptr = ROL::makePtrFromRef(bnd);


    /*******************************************************************************************/
    /* Form Primal-Dual system residual and operator then solve for direction vector           */
    /*******************************************************************************************/
    
          
    ROL::Ptr<V> rhs = CreatePartitionedVector(state->gradientVec,
                                         state->constraintVec,
                                         resL_,
                                         resU_);
    
    ROL::Ptr<V> sysvec = CreatePartitionedVector( x_, l_, zl_, zu_ );


    ROL::Ptr<RESIDUAL> residual = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL_,maskU_,w_,mu_,symmetrize_); 

    residual->value(*rhs,*sysvec,tol);

    ROL::Ptr<V> sol = CreatePartitionedVector( xscratch_, lscratch_, zlscratch_, zuscratch_ );

    LOPEC jacobian( sysvec, residual ); 

    
    krylov_->run(*sol,jacobian,*residual,*precond_,iterKrylov_,flagKrylov_); 
    
    /*******************************************************************************************/
    /* Perform line search                                                                     */
    /*******************************************************************************************/



  } // compute() 
  


  void update( Vector<Real> &x, Vector<Real> &l, const Vector<Real> &s,
               Objective<Real> &obj, Constraint<Real> &con, 
               BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {

    // Check deviation from shifted complementarity
    Elementwise::Shift<Real>    minus_mu(-mu_);

    xscratch_->set(x);
    xscratch_->applyBinary(mult,*zl_);
    xscratch_->applyUnary(minus_mu);

    lowerViolation = normLinf(*xscratch_)/scl_; // \f$ \max_i xz_l^i - \mu \f$
 
    xscratch_->set(x);
    xscratch_->applyBinary(mult,*zu_);
    xscratch_->applyUnary(minus_mu);            
    
    upperBound = normLinf(*xscratch_)/scu_;

    // Evaluate \f$E_\mu(x,\lambda,z_l,z_u)\f$
    Real Emu = algo_state.gnorm;
    Emu = std::max(Emu,algo_state.cnorm);
    Emu = std::max(Emu,upperBound);
    Emu = std::max(Emu,lowerBound);

    // If sufficiently converged for the current mu, update it
    if(Emu < (kappa_epsilon_*mu_) ) {
      Real mu_old = mu_;
 
      /* \mu_{j+1} = \max\left{ \frac{\epsilon_\text{tol}}{10},
                                \min\{ \kappa_{\mu} \mu_j, 
                                       \mu_j^{\theta_\mu}\} \right\} */
      mu_ = std::min(kappa_mu_*mu_old,std::pow(mu_old,theta_mu_));
      mu_ = std::max(eps_tol_/10.0,mu_);
 
     // Update fraction-to-boundary parameter
     tau_ = std::max(tau_min_,1.0-mu_);     

          
 
    } 

  } // update() 




  // TODO: Implement header print out
  std::string printHeader( void ) const {
    std::string head("");
    return head;
  }

  // TODO: Implement name print out
  std::string printName( void ) const {
    std::string name("");
    return name;
  }


}; // class PrimalDualInteriorPointStep

} // namespace InteriorPoint
} // namespace ROL


#endif //  ROL_PRIMALDUALINTERIORPOINTSTEP_H