ROL_Fletcher.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_FLETCHER_H
#define ROL_FLETCHER_H

#include "ROL_FletcherBase.hpp"
#include "ROL_Objective.hpp"
#include "ROL_Constraint.hpp"
#include "ROL_Vector.hpp"
#include "ROL_Types.hpp"
#include "ROL_Ptr.hpp"
#include "ROL_Krylov.hpp"
#include "ROL_PartitionedVector.hpp"
#include <iostream>

namespace ROL {

template <class Real>
class Fletcher : public FletcherBase<Real> {
private:
  // Required for Fletcher penalty function definition
  using FletcherBase<Real>::obj_;
  using FletcherBase<Real>::con_;

  using FletcherBase<Real>::penaltyParameter_;
  using FletcherBase<Real>::quadPenaltyParameter_;

  // Evaluation counters
  using FletcherBase<Real>::nfval_;
  using FletcherBase<Real>::ngval_;
  using FletcherBase<Real>::ncval_;

  using FletcherBase<Real>::fPhi_;     // value of penalty function
  using FletcherBase<Real>::gPhi_;     // gradient of penalty function

  using FletcherBase<Real>::y_;        // multiplier estimate

  using FletcherBase<Real>::fval_;     // value of objective function
  using FletcherBase<Real>::g_;        // gradient of objective value
  using FletcherBase<Real>::c_;        // constraint value
  using FletcherBase<Real>::scaledc_;  // penaltyParameter_ * c_
  using FletcherBase<Real>::gL_;       // gradient of Lagrangian (g - A*y)

  using FletcherBase<Real>::cnorm_;    // norm of constraint violation

  using FletcherBase<Real>::isValueComputed_;
  using FletcherBase<Real>::isGradientComputed_;
  using FletcherBase<Real>::isMultiplierComputed_;
  using FletcherBase<Real>::isObjValueComputed_;
  using FletcherBase<Real>::isObjGradComputed_;
  using FletcherBase<Real>::isConValueComputed_;

  using FletcherBase<Real>::delta_;    // regularization parameter

  int HessianApprox_;

  // Temporaries
  Ptr<Vector<Real> > Tv_;       // Temporary for matvecs
  Ptr<Vector<Real> > w_;        // first component of augmented system solve solution
  Ptr<Vector<Real> > v_;        // second component of augmented system solve solution
  Ptr<Vector<Real> > wg_;       // first component of augmented system solve solution for gradient
  Ptr<Vector<Real> > vg_;       // second component of augmented system solve solution for gradient

  Ptr<Vector<Real> > xzeros_;   // zero vector
  Ptr<Vector<Real> > czeros_;   // zero vector

  using FletcherBase<Real>::useInexact_;

  using FletcherBase<Real>::multSolverError_;         // Error from augmented system solve in value()
  using FletcherBase<Real>::gradSolveError_;          // Error from augmented system solve in gradient()

  // For Augmented system solves
  using FletcherBase<Real>::krylov_;
  using FletcherBase<Real>::iterKrylov_;
  using FletcherBase<Real>::flagKrylov_;
  using FletcherBase<Real>::v1_;
  using FletcherBase<Real>::v2_;
  using FletcherBase<Real>::vv_;
  using FletcherBase<Real>::w1_;
  using FletcherBase<Real>::w2_;
  using FletcherBase<Real>::ww_;
  using FletcherBase<Real>::b1_;
  using FletcherBase<Real>::b2_;
  using FletcherBase<Real>::bb_;

  class AugSystem : public LinearOperator<Real> {
  private:
    const Ptr<Constraint<Real> > con_;
    const Ptr<const Vector<Real> > x_;
    const Real delta_;
  public:
    AugSystem(const Ptr<Constraint<Real> > &con,
              const Ptr<const Vector<Real> > &x,
              const Real delta) : con_(con), x_(x), delta_(delta) {}

    void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
      PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
      const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);

      con_->applyAdjointJacobian(*(Hvp.get(0)), *(vp.get(1)), *x_, tol);
      Hvp.get(0)->plus(*(vp.get(0)));

      con_->applyJacobian(*(Hvp.get(1)), *(vp.get(0)), *x_, tol);
      Hvp.get(1)->axpy(-delta_*delta_, *(vp.get(1)));
    }
  };

  class AugSystemPrecond : public LinearOperator<Real> {
  private:
    const Ptr<Constraint<Real> > con_;
    const Ptr<const Vector<Real> > x_;
  public:
    AugSystemPrecond(const Ptr<Constraint<Real> > con,
                     const Ptr<const Vector<Real> > x) : con_(con), x_(x) {}

    void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
      Hv.set(v.dual());
    }
    void applyInverse(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
      Real zero(0);
      PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
      const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);

      Hvp.set(0, *(vp.get(0)));
      // Second x should be dual, but unused?
      con_->applyPreconditioner(*(Hvp.get(1)),*(vp.get(1)),*x_,*x_, zero); 
    }
  };

public:
  Fletcher(const ROL::Ptr<Objective<Real> > &obj,
           const ROL::Ptr<Constraint<Real> > &con,
           const Vector<Real> &optVec,
           const Vector<Real> &conVec,
           ROL::ParameterList &parlist)
      : FletcherBase<Real>(obj, con) {

      gPhi_    = optVec.dual().clone();
      y_       = conVec.dual().clone();
      g_       = optVec.dual().clone();
      gL_      = optVec.dual().clone();
      c_       = conVec.clone();
      scaledc_ = conVec.clone();

      Tv_ = optVec.dual().clone();
      w_ = optVec.dual().clone();
      v_ = conVec.dual().clone();
      wg_ = optVec.dual().clone();
      vg_ = conVec.dual().clone();

      xzeros_ = optVec.dual().clone();
      xzeros_->zero();
      czeros_ = conVec.clone();
      czeros_->zero();

      v1_ = optVec.dual().clone();
      v2_ = conVec.dual().clone();
      vv_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({v1_, v2_}));

      w1_ = optVec.dual().clone();
      w2_ = conVec.dual().clone();
      ww_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({w1_, w2_}));

      b1_ = optVec.dual().clone();
      b2_ = conVec.clone();
      bb_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({b1_, b2_}));

      ROL::ParameterList& sublist = parlist.sublist("Step").sublist("Fletcher");
      HessianApprox_ = sublist.get("Level of Hessian Approximation",  0);
      penaltyParameter_ = sublist.get("Penalty Parameter", 1.0);
      quadPenaltyParameter_ = sublist.get("Quadratic Penalty Parameter", 0.0);

      delta_ = sublist.get("Regularization Parameter", 0.0);

      useInexact_ = sublist.get("Inexact Solves", false);
      
      ROL::ParameterList krylovList;
      Real atol = static_cast<Real>(1e-12);
      Real rtol = static_cast<Real>(1e-2);
      krylovList.sublist("General").sublist("Krylov").set("Type", "GMRES");
      krylovList.sublist("General").sublist("Krylov").set("Absolute Tolerance", atol);
      krylovList.sublist("General").sublist("Krylov").set("Relative Tolerance", rtol);
      krylovList.sublist("General").sublist("Krylov").set("Iteration Limit", 200);
      krylov_ = KrylovFactory<Real>(krylovList);
  }

  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
    obj_->update(x,flag,iter);
    con_->update(x,flag,iter);
    isValueComputed_ = (flag ? false : isValueComputed_);
    isGradientComputed_ = (flag ? false : isGradientComputed_);
    isMultiplierComputed_ = (flag ? false : isMultiplierComputed_);
    isObjValueComputed_ = (flag ? false : isObjValueComputed_);
    isObjGradComputed_ = (flag ? false : isObjGradComputed_);
    isConValueComputed_ = (flag ? false : isConValueComputed_);
  }

  Real value( const Vector<Real> &x, Real &tol ) {
    if( isValueComputed_ && multSolverError_*cnorm_ <= tol) {
      tol = multSolverError_*cnorm_;
      return fPhi_;
    }

    Real zero(0);

    // Reset tolerances
    Real origTol = tol;
    Real tol2 = origTol;

    FletcherBase<Real>::objValue(x, tol2); tol2 = origTol;
    multSolverError_ = origTol / (static_cast<Real>(2) * std::max(static_cast<Real>(1), cnorm_));
    computeMultipliers(x, multSolverError_);
    tol =  multSolverError_*cnorm_;

    fPhi_ = fval_ - c_->dot(y_->dual());

    if( quadPenaltyParameter_ > zero ) {
      fPhi_ = fPhi_ + Real(0.5)*quadPenaltyParameter_*(c_->dot(c_->dual()));
    }

    isValueComputed_ = true;

    return fPhi_;
  }

  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
    if( isGradientComputed_ && gradSolveError_ <= tol) {
      tol = gradSolveError_;
      g.set(*gPhi_);
      return;
    }

    Real zero(0);

    // Reset tolerances
    Real origTol = tol;
    Real tol2 = origTol;

    gradSolveError_ = origTol / static_cast<Real>(2);
    computeMultipliers(x, gradSolveError_);

    bool refine = isGradientComputed_;

    // gPhi = sum y_i H_i w + sigma w + sum v_i H_i gL - H w + gL
    solveAugmentedSystem( *wg_, *vg_, *xzeros_, *c_, x, gradSolveError_, refine );
    gradSolveError_ += multSolverError_;
    tol = gradSolveError_;

    con_->applyAdjointHessian( *gPhi_, *y_, *wg_, x, tol2 ); tol2 = origTol;
    gPhi_->axpy( penaltyParameter_, *wg_ );

    obj_->hessVec( *Tv_, *wg_, x, tol2 ); tol2 = origTol;
    gPhi_->axpy( static_cast<Real>(-1), *Tv_ );

    con_->applyAdjointHessian( *Tv_, *vg_, *gL_, x, tol2 ); tol2 = origTol;
    gPhi_->plus( *Tv_ );

    gPhi_->plus( *gL_ );

    if( quadPenaltyParameter_ > zero ) {
      con_->applyAdjointJacobian( *Tv_, *c_, x, tol2 ); tol2 = origTol;
      gPhi_->axpy( quadPenaltyParameter_, *Tv_ );
    }

    g.set(*gPhi_);
    isGradientComputed_ = true;
  }

  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
    Real zero(0);

    // Reset tolerances
    Real origTol = tol;
    Real tol2 = origTol;

    if( !isMultiplierComputed_ || !useInexact_) {
      // hessVec tol is always set to ~1e-8. So if inexact linear system solves are used, then
      // the multipliers will always get re-evaluated to high precision, which defeats the purpose
      // of computing the objective/gradient approximately.
      // Thus if inexact linear system solves are used, we will not update the multiplier estimates
      // to high precision.
      computeMultipliers(x, tol);
    }

    obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
    con_->applyAdjointHessian( *Tv_, *y_, v, x, tol2 ); tol2 = origTol;
    hv.axpy(static_cast<Real>(-1), *Tv_ );

    tol2 = tol;
    solveAugmentedSystem( *w_, *v_, hv, *czeros_, x, tol2 ); tol2 = origTol;
    hv.scale( static_cast<Real>(-1) );
    hv.plus( *w_ );

    Tv_->set(v);
    tol2 = tol;
    solveAugmentedSystem( *w_, *v_, *Tv_, *czeros_, x, tol2 ); tol2 = origTol;
    hv.axpy(static_cast<Real>(-2)*penaltyParameter_, *w_);

    obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = origTol;
    hv.plus( *Tv_ );
    con_->applyAdjointHessian( *Tv_, *y_, *w_, x, tol2 ); tol2 = origTol;
    hv.axpy( static_cast<Real>(-1), *Tv_ );

    hv.axpy( static_cast<Real>(2)*penaltyParameter_, v );

    if( quadPenaltyParameter_ > zero ) {
      con_->applyJacobian( *b2_, v, x, tol2 ); tol2 = origTol;
      con_->applyAdjointJacobian( *Tv_, *b2_, x, tol2 ); tol2 = origTol;
      hv.axpy( quadPenaltyParameter_, *Tv_ );
      con_->applyAdjointHessian( *Tv_, *c_, v, x, tol2); tol2 = origTol;
      hv.axpy( -quadPenaltyParameter_, *Tv_ );
    }

  }

  void solveAugmentedSystem(Vector<Real> &v1,
                            Vector<Real> &v2,
                            const Vector<Real> &b1,
                            const Vector<Real> &b2,
                            const Vector<Real> &x,
                            Real &tol,
                            bool refine = false) {
    // Ignore tol for now
    ROL::Ptr<LinearOperator<Real> > K
      = ROL::makePtr<AugSystem>(con_, makePtrFromRef(x), delta_);
    ROL::Ptr<LinearOperator<Real> > P
      = ROL::makePtr<AugSystemPrecond>(con_, makePtrFromRef(x));

    b1_->set(b1);
    b2_->set(b2);

    if( refine ) {
      // TODO: Make sure this tol is actually ok...
      Real origTol = tol;
      w1_->set(v1);
      w2_->set(v2);
      K->apply(*vv_, *ww_, tol); tol = origTol;

      b1_->axpy( static_cast<Real>(-1), *v1_ );
      b2_->axpy( static_cast<Real>(-1), *v2_ );
    }

    v1_->zero();
    v2_->zero();

    // If inexact, change tolerance
    if( useInexact_ ) {
      krylov_->resetAbsoluteTolerance(tol);
    }

    flagKrylov_ = 0;
    tol = krylov_->run(*vv_,*K,*bb_,*P,iterKrylov_,flagKrylov_);

    if( refine ) {
      v1.plus(*v1_);
      v2.plus(*v2_);
    } else {
      v1.set(*v1_);
      v2.set(*v2_);
    }
  }

  void computeMultipliers(const Vector<Real>& x, const Real tol) {
    if( isMultiplierComputed_ && multSolverError_ <= tol) {
      return;
    }

    if( !isMultiplierComputed_ ) {
      Real tol2 = tol;
      FletcherBase<Real>::objGrad(x, tol2); tol2 = tol;
      FletcherBase<Real>::conValue(x, tol2);
      cnorm_ = c_->norm();
    }

    bool refine = isMultiplierComputed_;

    multSolverError_ = tol;
    solveAugmentedSystem(*gL_, *y_, *g_, *scaledc_, x, multSolverError_, refine);

    isMultiplierComputed_ = true;
  }

}; // class Fletcher

} // namespace ROL

#endif