C++ Objective::update方法代码示例

 void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
           const Real &gs, const Vector<Real> &s, const Vector<Real> &x, 
           Objective<Real> &obj, BoundConstraint<Real> &con ) {
   Real tol = std::sqrt(ROL_EPSILON<Real>());
   ls_neval = 0;
   ls_ngrad = 0;
   // Get initial line search parameter
   alpha = LineSearch<Real>::getInitialAlpha(ls_neval,ls_ngrad,fval,gs,x,s,obj,con);
   // Update iterate
   LineSearch<Real>::updateIterate(*xnew_,x,s,alpha,con);
   // Get objective value at xnew
   Real fold = fval;
   obj.update(*xnew_);
   fval = obj.value(*xnew_,tol);
   ls_neval++;
   // Perform backtracking
   while ( !LineSearch<Real>::status(LINESEARCH_BACKTRACKING,ls_neval,ls_ngrad,alpha,fold,gs,fval,*xnew_,s,obj,con) ) {
     alpha *= rho_;
     // Update iterate
     LineSearch<Real>::updateIterate(*xnew_,x,s,alpha,con);
     // Get objective value at xnew
     obj.update(*xnew_);
     fval = obj.value(*xnew_,tol);
     ls_neval++;
   }
 }

  void update( Vector<Real> &x, const Vector<Real> &s,
               Objective<Real> &obj, BoundConstraint<Real> &bnd,
               AlgorithmState<Real> &algo_state ) {
    Real tol = std::sqrt(ROL_EPSILON<Real>());
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();

    // Update iterate
    algo_state.iter++;
    x.plus(s);
    (step_state->descentVec)->set(s);
    algo_state.snorm = s.norm();

    // Compute new gradient
    if ( useSecantPrecond_ ) {
      gp_->set(*(step_state->gradientVec));
    }
    obj.update(x,true,algo_state.iter);
    if ( computeObj_ ) {
      algo_state.value = obj.value(x,tol);
      algo_state.nfval++;
    }
    obj.gradient(*(step_state->gradientVec),x,tol);
    algo_state.ngrad++;

    // Update Secant Information
    if ( useSecantPrecond_ ) {
      secant_->updateStorage(x,*(step_state->gradientVec),*gp_,s,algo_state.snorm,algo_state.iter+1);
    }

    // Update algorithm state
    (algo_state.iterateVec)->set(x);
    algo_state.gnorm = step_state->gradientVec->norm();
  }

  virtual Real getInitialAlpha(int &ls_neval, int &ls_ngrad, const Real fval, const Real gs, 
                               const Vector<Real> &x, const Vector<Real> &s, 
                               Objective<Real> &obj, BoundConstraint<Real> &con) {
    Real val = 1.0;
    if (useralpha_) {
      val = alpha0_;
    }
    else {
      if (edesc_ == DESCENT_STEEPEST || edesc_ == DESCENT_NONLINEARCG) {
        Real tol = std::sqrt(ROL_EPSILON);
        Real alpha = 1.0;
        // Evaluate objective at x + s
        updateIterate(*d_,x,s,alpha,con);
        obj.update(*d_);
        Real fnew = obj.value(*d_,tol);
        ls_neval++;
        // Minimize quadratic interpolate to compute new alpha
        alpha = -gs/(2.0*(fnew-fval-gs));
        val = ((std::abs(alpha) > std::sqrt(ROL_EPSILON)) ? std::abs(alpha) : 1.0);

        alpha0_ = val;
        useralpha_ = true;
      }
      else {
        val = 1.0;
      }
    }
    return val;
  }

 // Run Iteration scaled line search
 void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
           const Real &gs, const Vector<Real> &s, const Vector<Real> &x, 
           Objective<Real> &obj, BoundConstraint<Real> &con ) {
   Real tol = std::sqrt(ROL_EPSILON<Real>());
   ls_neval = 0;
   ls_ngrad = 0;
   // Update target objective value
   if ( fval < min_value_ ) {
     min_value_ = fval;
   }
   target_ = rec_value_ - 0.5*delta_;
   if ( fval < target_ ) {
     rec_value_ = min_value_; 
     sigma_ = 0.0;
   }
   else {
     if ( sigma_ > bound_ ) {
       rec_value_ = min_value_;
       sigma_ = 0.0;
       delta_ *= 0.5;
     }
   }
   target_ = rec_value_ - delta_;
   // Get line-search parameter
   alpha = (fval - target_)/std::abs(gs);
   // Update iterate
   LineSearch<Real>::updateIterate(*xnew_,x,s,alpha,con);
   // Compute objective function value
   obj.update(*xnew_);
   fval = obj.value(*xnew_,tol);
   ls_neval++;
   // Update sigma 
   sigma_ += alpha*std::sqrt(std::abs(gs));
 }

 virtual Real getInitialAlpha(int &ls_neval, int &ls_ngrad, const Real fval, const Real gs, 
                              const Vector<Real> &x, const Vector<Real> &s, 
                              Objective<Real> &obj, BoundConstraint<Real> &con) {
   Real val = 1.0;
   if (useralpha_) {
     val = alpha0_;
   }
   else {
     if (edesc_ == DESCENT_STEEPEST || edesc_ == DESCENT_NONLINEARCG) {
       Real tol = std::sqrt(ROL_EPSILON);
       Real alpha = 1.0;
       // Evaluate objective at x + s
       updateIterate(*d_,x,s,alpha,con);
       obj.update(*d_);
       Real fnew = obj.value(*d_,tol);
       ls_neval++;
       // Minimize quadratic interpolate to compute new alpha
       alpha = -gs/(2.0*(fnew-fval-gs));
       // Evaluate objective at x + alpha s 
       updateIterate(*d_,x,s,alpha,con);
       obj.update(*d_);
       fnew = obj.value(*d_,tol);
       ls_neval++;
       // Ensure that sufficient decrease and curvature conditions are satisfied
       bool stat = status(LINESEARCH_BISECTION,ls_neval,ls_ngrad,alpha,fval,gs,fnew,x,s,obj,con);
       if ( !stat ) {
         alpha = 1.0;
       }
       val = alpha;
     }
     else {
       val = 1.0;
     }
   }
   return val;
 }

 // Run Iteration scaled line search
 void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
           const Real &gs, const Vector<Real> &s, const Vector<Real> &x, 
           Objective<Real> &obj, BoundConstraint<Real> &con ) {
   Real tol = std::sqrt(ROL_EPSILON);
   ls_neval = 0;
   ls_ngrad = 0;
   // Get line search parameter
   algo_iter_++;
   alpha = LineSearch<Real>::getInitialAlpha(ls_neval,ls_ngrad,fval,gs,x,s,obj,con)/algo_iter_;
   // Update iterate
   LineSearch<Real>::updateIterate(*xnew_,x,s,alpha,con);
   // Compute objective function value
   obj.update(*xnew_);
   fval = obj.value(*xnew_,tol);
   ls_neval++;
 }

  void update( Vector<Real> &x, const Vector<Real> &s,
               Objective<Real> &obj, BoundConstraint<Real> &bnd,
               AlgorithmState<Real> &algo_state ) {
    Real tol = std::sqrt(ROL_EPSILON<Real>()), one(1);
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();

    // Update iterate and store previous step
    algo_state.iter++;
    d_->set(x);
    x.plus(s);
    bnd.project(x);
    (step_state->descentVec)->set(x);
    (step_state->descentVec)->axpy(-one,*d_);
    algo_state.snorm = s.norm();

    // Compute new gradient
    gp_->set(*(step_state->gradientVec));
    obj.update(x,true,algo_state.iter);
    if ( computeObj_ ) {
      algo_state.value = obj.value(x,tol);
      algo_state.nfval++;
    }
    obj.gradient(*(step_state->gradientVec),x,tol);
    algo_state.ngrad++;

    // Update Secant Information
    secant_->updateStorage(x,*(step_state->gradientVec),*gp_,s,algo_state.snorm,algo_state.iter+1);

    // Update algorithm state
    (algo_state.iterateVec)->set(x);
    if ( useProjectedGrad_ ) {
      gp_->set(*(step_state->gradientVec));
      bnd.computeProjectedGradient( *gp_, x );
      algo_state.gnorm = gp_->norm();
    }
    else {
      d_->set(x);
      d_->axpy(-one,(step_state->gradientVec)->dual());
      bnd.project(*d_);
      d_->axpy(-one,x);
      algo_state.gnorm = d_->norm();
    }
  }

  /** \brief Initialize step.  

             This includes projecting the initial guess onto the constraints, 
             computing the initial objective function value and gradient, 
             and initializing the dual variables.

             @param[in,out]    x           is the initial guess 
             @param[in]        obj         is the objective function
             @param[in]        con         are the bound constraints
             @param[in]        algo_state  is the current state of the algorithm
  */
  void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g, 
                   Objective<Real> &obj, BoundConstraint<Real> &con, 
                   AlgorithmState<Real> &algo_state ) {
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();
    // Initialize state descent direction and gradient storage
    step_state->descentVec  = s.clone();
    step_state->gradientVec = g.clone();
    step_state->searchSize  = 0.0;
    // Initialize additional storage
    xlam_ = x.clone(); 
    x0_   = x.clone();
    xbnd_ = x.clone();
    As_   = s.clone(); 
    xtmp_ = x.clone(); 
    res_  = g.clone();
    Ag_   = g.clone(); 
    rtmp_ = g.clone(); 
    gtmp_ = g.clone(); 
    // Project x onto constraint set
    con.project(x);
    // Update objective function, get value, and get gradient
    Real tol = std::sqrt(ROL_EPSILON);
    obj.update(x,true,algo_state.iter);
    algo_state.value = obj.value(x,tol);
    algo_state.nfval++;
    algo_state.gnorm = computeCriticalityMeasure(x,obj,con,tol);
    algo_state.ngrad++;
    // Initialize dual variable
    lambda_ = s.clone(); 
    lambda_->set((step_state->gradientVec)->dual());
    lambda_->scale(-1.0);
    //con.setVectorToLowerBound(*lambda_);
    // Initialize Hessian and preconditioner
    Teuchos::RCP<Objective<Real> > obj_ptr = Teuchos::rcp(&obj, false);
    Teuchos::RCP<BoundConstraint<Real> > con_ptr = Teuchos::rcp(&con, false);
    hessian_ = Teuchos::rcp( 
      new PrimalDualHessian<Real>(secant_,obj_ptr,con_ptr,algo_state.iterateVec,xlam_,useSecantHessVec_) );
    precond_ = Teuchos::rcp( 
      new PrimalDualPreconditioner<Real>(secant_,obj_ptr,con_ptr,algo_state.iterateVec,xlam_,
                                         useSecantPrecond_) );
  }

  /** \brief Update step, if successful.

             This function returns \f$x_{k+1} = x_k + s_k\f$.
             It also updates secant information if being used.

             @param[in]        x           is the new iterate
             @param[out]       s           is the step computed via PDAS
             @param[in]        obj         is the objective function
             @param[in]        con         are the bound constraints
             @param[in]        algo_state  is the current state of the algorithm
  */
  void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj, BoundConstraint<Real> &con,
               AlgorithmState<Real> &algo_state ) {
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();

    x.plus(s);
    feasible_ = con.isFeasible(x);
    algo_state.snorm = s.norm();
    algo_state.iter++;
    Real tol = std::sqrt(ROL_EPSILON);
    obj.update(x,true,algo_state.iter);
    algo_state.value = obj.value(x,tol);
    algo_state.nfval++;
    
    if ( secant_ != Teuchos::null ) {
      gtmp_->set(*(step_state->gradientVec));
    }
    algo_state.gnorm = computeCriticalityMeasure(x,obj,con,tol);
    algo_state.ngrad++;

    if ( secant_ != Teuchos::null ) {
      secant_->update(*(step_state->gradientVec),*gtmp_,s,algo_state.snorm,algo_state.iter+1);
    }
    (algo_state.iterateVec)->set(x);
  }

    void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj, BoundConstraint<Real> &con,
                 AlgorithmState<Real> &algo_state ) {
        Real tol = std::sqrt(ROL_EPSILON<Real>());
        Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();

        // Update iterate and store step
        algo_state.iter++;
        x.plus(s);
        (step_state->descentVec)->set(s);
        algo_state.snorm = s.norm();

        // Compute new gradient
        obj.update(x,true,algo_state.iter);
        if ( computeObj_ ) {
            algo_state.value = obj.value(x,tol);
            algo_state.nfval++;
        }
        obj.gradient(*(step_state->gradientVec),x,tol);
        algo_state.ngrad++;

        // Update algorithm state
        (algo_state.iterateVec)->set(x);
        algo_state.gnorm = (step_state->gradientVec)->norm();
    }

 void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
           const Real &gs, const Vector<Real> &s, const Vector<Real> &x, 
           Objective<Real> &obj, BoundConstraint<Real> &con ) {
   Real tol = std::sqrt(ROL_EPSILON);
   ls_neval = 0;
   ls_ngrad = 0;
   // Get initial line search parameter
   alpha = LineSearch<Real>::getInitialAlpha(ls_neval,ls_ngrad,fval,gs,x,s,obj,con);
   // Update iterate
   LineSearch<Real>::updateIterate(*xnew_,x,s,alpha,con);
   // Get objective value at xnew
   Real fold = fval;
   obj.update(*xnew_);
   fval = obj.value(*xnew_,tol);
   ls_neval++;
   // Initialize
   Real fvalp  = 0.0;
   Real alpha1 = 0.0;
   Real alpha2 = 0.0;
   Real a      = 0.0;
   Real b      = 0.0;
   Real x1     = 0.0;
   Real x2     = 0.0;
   bool first_iter = true;
   // Perform cubic interpolation back tracking
   while ( !LineSearch<Real>::status(LINESEARCH_CUBICINTERP,ls_neval,ls_ngrad,alpha,fold,gs,fval,x,s,obj,con) ) {
     if ( first_iter ) { // Minimize quadratic interpolate
       alpha1 = -gs*alpha*alpha/(2.0*(fval-fold-gs*alpha));
       first_iter = false;
     }
     else {              // Minimize cubic interpolate
       x1 = fval-fold-alpha*gs;
       x2 = fvalp-fval-alpha2*gs;
       a = (1.0/(alpha - alpha2))*( x1/(alpha*alpha) - x2/(alpha2*alpha2));
       b = (1.0/(alpha - alpha2))*(-x1*alpha2/(alpha*alpha) + x2*alpha/(alpha2*alpha2));
       if ( std::abs(a) < ROL_EPSILON ) {
         alpha1 = -gs/(2.0*b);
       }
       else {
         alpha1 = (-b+sqrt(b*b-3.0*a*gs))/(3.0*a);
       }
       if ( alpha1 > 0.5*alpha ) {
         alpha1 = 0.5*alpha;
       }
     }
     alpha2 = alpha;
     fvalp  = fval;
     // Back track if necessary
     if ( alpha1 <= 0.1*alpha ) {
       alpha *= 0.1;
     }
     else if ( alpha1 >= 0.5*alpha ) {
       alpha *= 0.5;
     }
     else {
       alpha = alpha1;
     }
     // Update iterate
     LineSearch<Real>::updateIterate(*xnew_,x,s,alpha,con);
     // Get objective value at xnew
     obj.update(*xnew_);
     fval = obj.value(*xnew_,tol);
     ls_neval++;
   }
 }

  void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
            const Real &gs, const Vector<Real> &s, const Vector<Real> &x, 
            Objective<Real> &obj, BoundConstraint<Real> &con ) {
    Real tol = std::sqrt(ROL_EPSILON);
    ls_neval = 0;
    ls_ngrad = 0;
    // Get initial line search parameter
    alpha = LineSearch<Real>::getInitialAlpha(ls_neval,ls_ngrad,fval,gs,x,s,obj,con);

    // Compute value phi(0)
    Real tl = 0.0;         // Left interval point
    Real val_tl = fval;

    // Initialize value phi(t)
    Real tc = 0.0;        // Center interval point
    Real val_tc = 0.0;

    // Compute value phi(alpha)
    Real tr = alpha;      // Right interval point
    LineSearch<Real>::updateIterate(*xnew_,x,s,tr,con);
    obj.update(*xnew_);
    Real val_tr = obj.value(*xnew_,tol);
    ls_neval++;

    // Check if phi(alpha) < phi(0)
    if ( val_tr < val_tl ) {
      if ( LineSearch<Real>::status(LINESEARCH_BRENTS,ls_neval,ls_ngrad,tr,fval,gs,val_tr,x,s,obj,con) ) {
        alpha = tr;
        fval  = val_tr;
        return;
      }
    }

    // Compute min( phi(0), phi(alpha) )
    Real t     = 0.0;
    Real val_t = 0.0;
    if ( val_tl < val_tr ) {
      t     = tl;
      val_t = val_tl;
    }
    else {
      t     = tr;
      val_t = val_tr;
    }

    // Determine bracketing triple
    const Real gr                = (1.0+sqrt(5.0))/2.0;
    const Real inv_gr2           = 1.0/(gr*gr);
    const Real goldinv           = 1.0/(1.0+gr);
    const Real tiny              = sqrt(ROL_EPSILON);
    const Real max_extrap_factor = 100.0;
    Real tmp    = 0.0;
    Real q      = 0.0;
    Real r      = 0.0; 
    Real tm     = 0.0;
    Real tlim   = 0.0; 
    Real val_tm = 0.0;

    int itbt = 0;
    while ( val_tr > val_tl && itbt < 8 ) {
      tc     = tr;
      val_tc = val_tr;

      tr     = goldinv * (tc + gr*tl);
      LineSearch<Real>::updateIterate(*xnew_,x,s,tr,con);
      obj.update(*xnew_);
      val_tr = obj.value(*xnew_,tol);
      ls_neval++;

      itbt++;
    }
    if ( val_tr > val_tl ) {
      tmp    = tl;
      tl     = tr;
      tr     = tmp;
      tmp    = val_tr;
      val_tr = val_tl;
      val_tl = tmp;
      tc     = 0.0;
    }

    if ( std::abs(tc) < ROL_EPSILON ) {
      tc = tl + (gr-1.0)*(tr-tl);
      LineSearch<Real>::updateIterate(*xnew_,x,s,tc,con);
      obj.update(*xnew_);
      val_tc = obj.value(*xnew_,tol);
      ls_neval++;
    }

    if ( val_tl <= val_tr && val_tl <= val_tc ) {
      t     = tl;
      val_t = val_tl;
    }
    else if ( val_tc <= val_tr && val_tc <= val_tl ) {
      t     = tc;
      val_t = val_tc;
    }
    else {
      t     = tr;
      val_t = val_tr;
//.........这里部分代码省略.........

  virtual bool status( const ELineSearch type, int &ls_neval, int &ls_ngrad, const Real alpha, 
                       const Real fold, const Real sgold, const Real fnew, 
                       const Vector<Real> &x, const Vector<Real> &s, 
                       Objective<Real> &obj, BoundConstraint<Real> &con ) { 
    Real tol = std::sqrt(ROL_EPSILON);

    // Check Armijo Condition
    bool armijo = false;
    if ( con.isActivated() ) {
      Real gs = 0.0;
      if ( edesc_ == DESCENT_STEEPEST ) {
        updateIterate(*d_,x,s,alpha,con);
        d_->scale(-1.0);
        d_->plus(x);
        gs = -s.dot(*d_);
      }
      else {
        d_->set(s);
        d_->scale(-1.0);
        con.pruneActive(*d_,*(grad_),x,eps_);
        gs = alpha*(grad_)->dot(*d_);
        d_->zero();
        updateIterate(*d_,x,s,alpha,con);
        d_->scale(-1.0);
        d_->plus(x);
        con.pruneInactive(*d_,*(grad_),x,eps_);
        gs += d_->dot(grad_->dual());
      }
      if ( fnew <= fold - c1_*gs ) {
        armijo = true;
      }
    }
    else {
      if ( fnew <= fold + c1_*alpha*sgold ) {
        armijo = true;
      }
    }

    // Check Maximum Iteration
    bool itcond = false;
    if ( ls_neval >= maxit_ ) { 
      itcond = true;
    }

    // Check Curvature Condition
    bool curvcond = false;
    if ( armijo && ((type != LINESEARCH_BACKTRACKING && type != LINESEARCH_CUBICINTERP) ||
                    (edesc_ == DESCENT_NONLINEARCG)) ) {
      if (econd_ == CURVATURECONDITION_GOLDSTEIN) {
        if (fnew >= fold + (1.0-c1_)*alpha*sgold) {
          curvcond = true;
        }
      }
      else if (econd_ == CURVATURECONDITION_NULL) {
        curvcond = true;
      }
      else { 
        updateIterate(*xtst_,x,s,alpha,con);
        obj.update(*xtst_);
        obj.gradient(*g_,*xtst_,tol);
        Real sgnew = 0.0;
        if ( con.isActivated() ) {
          d_->set(s);
          d_->scale(-alpha);
          con.pruneActive(*d_,s,x);
          sgnew = -d_->dot(g_->dual());
        }
        else {
          sgnew = s.dot(g_->dual());
        }
        ls_ngrad++;
   
        if (    ((econd_ == CURVATURECONDITION_WOLFE)       
                     && (sgnew >= c2_*sgold))
             || ((econd_ == CURVATURECONDITION_STRONGWOLFE) 
                     && (std::abs(sgnew) <= c2_*std::abs(sgold)))
             || ((econd_ == CURVATURECONDITION_GENERALIZEDWOLFE) 
                     && (c2_*sgold <= sgnew && sgnew <= -c3_*sgold))
             || ((econd_ == CURVATURECONDITION_APPROXIMATEWOLFE) 
                     && (c2_*sgold <= sgnew && sgnew <= (2.0*c1_ - 1.0)*sgold)) ) {
          curvcond = true;
        }
      }
    }

    if (type == LINESEARCH_BACKTRACKING || type == LINESEARCH_CUBICINTERP) {
      if (edesc_ == DESCENT_NONLINEARCG) {
        return ((armijo && curvcond) || itcond);
      }
      else {
        return (armijo || itcond);
      }
    }
    else {
      return ((armijo && curvcond) || itcond);
    }
  }