TRGNAlg.hh

Go to the documentation of this file.

#ifndef RVL_UMIN_TRALG_
#define RVL_UMIN_TRALG_

#include "table.hh"
#include "op.hh"
#include "alg.hh"
#include "ioterm.hh"
#include "boolterm.hh"

namespace RVLUmin {

  using namespace RVL;
  using namespace RVLAlg;

  template<typename Scalar, typename Policy >
  class TRGNStep: public Algorithm {
    
    typedef typename ScalarFieldTraits<Scalar>::AbsType atype;
    
  public:
    
    void run() {
      //cerr<<"TRGNStep::run\n";
      // run then delete solver - first check whether trust region was
      // transgressed
      solver->run();  
      bool active = false;
      {
    Terminator & tsolver = dynamic_cast<Terminator &>(*solver);
    active = tsolver.query();
      }
      delete solver;    
      solver=NULL;

      // test for nontrivial step - failure is fatal
      if (p.norm()<minstep*pol.Delta) {
    nostep=true;
    return;
      }

      // here rnorm is final residual from TRLS solver
      atype jvalsave = jval;
      predred = jvalsave - 0.5*rnorm*rnorm;
      str<<"** predicted reduction = "<<predred<<" active = "<<active<<"\n";
      str<<"** length of step = "<<p.norm()<<"\n";
      //cerr<<"TRGNStep: predred="<<predred<<endl;
      // provisionally update eval point - save current point in case 
      // need to revert
      Vector<Scalar> & x = opeval.getPoint();
      xsave.copy(x);
      x.linComb(-ScalarFieldTraits<Scalar>::One(),p);
      // reset (reconstruct) solver
      //cerr<<"TRGNStep::run - main callg before call to pol.build Delta="<<pol.Delta<<endl;
      //      solver = pol.build(p,opeval.getDeriv(),opeval.getValue(),rnorm,nrnorm,str);   
      // by convention this should have recalculated the objective
      // function, gradient norm
      //      jval = 0.5*rnorm*rnorm;
      jval = 0.5*opeval.getValue().normsq();

      // compute actual reduction
      actred = jvalsave - jval;
      // Case 1: fail to satisfy G-A : rescind eval point update, shrink Delta
      str<<"** actual reduction: old jval = "<<jvalsave<<" new jval = "<<jval<<" actred = "<<actred<<"\n";
      if (actred < eta1*predred) {
    pol.Delta *= gamma1;
    x.copy(xsave);
      }
      // Case 2: otherwise keep update
      // Case 2.1: only grow Delta if step hit TR boundary
      if ((actred > eta2*predred) && active) { 
    int dcount = 0;
    Vector<Scalar> plast(opeval.getDomain());
    while (active && dcount<5) {
      str<<"active trust region constraint and good step - attempt to increase TR\n";
      atype deltalast=pol.Delta;
      atype jvallast = jval;
      
      // save last update vector
      plast.copy(p);
      // stretch it
      p.scale(gamma2);
      // stretch TR
      pol.Delta *= gamma2; 
      str<<"last TR = "<<deltalast<<" new TR = "<<pol.Delta<<"\n";
      // restore x
      x.copy(xsave);
      // compute DF(x)p - x not updated
      Vector<Scalar> dp(opeval.getRange());
      opeval.getDeriv().applyOp(p,dp);
      // create quadratic residual
      dp.linComb(-ScalarFieldTraits<Scalar>::One(),opeval.getValue());
      // estimate reduction - note jvalsave is resid at xsave
      predred = jvalsave - 0.5*dp.normsq();
      str<<"new predred="<<predred<<"\n";
      // update x
      x.linComb(-ScalarFieldTraits<Scalar>::One(),p);
      // create real residual
      jval = 0.5*opeval.getValue().normsq();
      str<<"last jval = "<<jvalsave<<" new jval = "<<jval<<"\n";
      // compute actual reduction
      actred = jvalsave - jval;
      str<<"new actred = "<<actred<<"\n";
      // if G-A fails, back off
      if (actred < eta1*predred || jval > jvallast) {
        str<<"G-A failed, revert to previous estimate, unset active flag\n";
        pol.Delta = deltalast;
        jval=jvallast;
        p.copy(plast);
        x.copy(xsave);
        x.linComb(-ScalarFieldTraits<Scalar>::One(),p);
        active=false;
      }   
      else {
        str<<"G-A passed";
      }
      dcount++;
      if (dcount==5) {
        str<<", max number of TR expansions reached\n";
      }
      else {
        str<<", try it again\n";
      }
    }
      }

      // rejigger everything for the new point
      solver = pol.build(p,opeval.getDeriv(),opeval.getValue(),rnorm,nrnorm,str);         
      agnrm=nrnorm;
      rgnrm=nrnorm*gnrmrecip;

      //cerr<<"at end of TRGNStep::run - Delta = "<<pol.Delta<<endl;
    }
    
    atype getDelta() { return pol.Delta; }

    TRGNStep(Policy const & _pol,
         OperatorEvaluation<Scalar> & _opeval,
         atype _eta1,
         atype _eta2,
         atype _gamma1,
         atype _gamma2,
         atype _minstep,
         bool & _nostep,
         atype & _predred,
         atype & _actred,
         atype & _jval,
         atype & _agnrm,
         atype & _rgnrm,
         ostream & _str) 
      : pol(_pol),
    opeval(_opeval), 
    eta1(_eta1),eta2(_eta2),gamma1(_gamma1),gamma2(_gamma2),minstep(_minstep),
    nostep(_nostep),predred(_predred), actred(_actred),
    jval(_jval), agnrm(_agnrm), rgnrm(_rgnrm),
    p(_opeval.getDomain()), xsave(_opeval.getDomain()), solver(NULL), str(_str) {

      //cerr<<"TRGNStep:: construction\n";
      // sanity test params
      if ( (gamma1 <= ScalarFieldTraits<atype>::Zero()) ||
       (gamma1 > ScalarFieldTraits<atype>::One())  ||
       (gamma2 < ScalarFieldTraits<atype>::One())  ||
       (gamma1*gamma2 > ScalarFieldTraits<atype>::One()-100.0*numeric_limits<atype>::epsilon())  ||
       (eta1  < ScalarFieldTraits<atype>::Zero())  ||
       (eta1 > eta2) ) {
    RVLException e;
    e<<"Error: TRTRGNStep constructor\n";
    e<<"insane TR param inputs\n";
    e<<"gamma1= "<<gamma1<<"\n";
    e<<"gamma2= "<<gamma2<<"\n";
    e<<"eta1  = "<<eta1<<"\n";
    e<<"eta2  = "<<eta2<<"\n";
    throw e;
      }

      // no step taken yet so set reductions to zero
      actred = ScalarFieldTraits<atype>::Zero();
      predred = ScalarFieldTraits<atype>::Zero();

      //cerr<<"TRGNStep -> Solver construction via policy\n";
      // initial solver assignment
      solver = pol.build(p,opeval.getDeriv(),opeval.getValue(),rnorm,nrnorm,str);
      //cerr<<"TRGNStep -> return from Solver construction\n";
      // solver policy: least squares solution of Ax=b. 
      // pass b=opeval.getValue, A=opeval.getDeriv references - cause reinitialization if needed
      // on creation, initialize x=0, r=b, rnorm=b.norm, nrnorm=A^Tb.norm 
      jval = 0.5*rnorm*rnorm;
      agnrm = nrnorm;

      // record initial reciprocal grad norm for future ref
      if (ProtectedDivision<atype>(ScalarFieldTraits<atype>::One(),nrnorm,gnrmrecip)) {
    // reset grad norms to zero
    agnrm = ScalarFieldTraits<atype>::Zero();
    rgnrm = ScalarFieldTraits<atype>::Zero();
    gnrmrecip = ScalarFieldTraits<atype>::Zero();
    str<<"NOTE: TRGNStep constructor\n";
    str<<"  initial value of normal residual norm = "<<nrnorm<<" below division threshold\n";
    str<<"  set residual norm values to zero, return - this should stop iteration\n"; 
      }
      else {
    rgnrm = ScalarFieldTraits<atype>::One();
      }
      // don't touch initial Delta
      //cerr<<"Exit TRGNStep constructor\n";

    }

    ~TRGNStep() { if (solver) delete solver; }

  private:
    
    Policy const & pol;
    OperatorEvaluation<Scalar> & opeval;
    atype & predred;
    atype & actred;
    atype & jval;
    atype & agnrm;
    atype & rgnrm;
    bool & nostep;
    atype gnrmrecip;
    atype eta1;
    atype eta2;
    atype gamma1;
    atype gamma2;
    atype minstep;
    atype rnorm;
    atype nrnorm;
    Vector<Scalar> p;
    Vector<Scalar> xsave;
    // workspace for solver
    mutable Algorithm * solver;
    ostream & str;
  };

  template<typename Scalar, typename Policy>
  class TRGNAlg: public Algorithm, public Policy {

    typedef typename ScalarFieldTraits<Scalar>::AbsType atype;

  public:

    void run() {
      //cerr<<"TRGNAlg::run\n";
      try {
    bool nostep = false;
    //cerr<<"TRGNAlg->TRGNStep constructor\n";
    TRGNStep<Scalar, Policy> step(*this,opeval,eta1,eta2,gamma1,gamma2,minstep,
                      nostep,predred,actred,jval,agnrm,rgnrm,str);
    BoolTerminator bstop(nostep);
    VectorCountingThresholdIterationTable<atype> vstop(maxcount,names,nums,tols,str);
    vstop.init();
    atype jval0=jval;
    atype agnrm0=agnrm;
    atype delta0=step.getDelta();
    OrTerminator stop(bstop,vstop);
    //cerr<<"TRGNAlg->LoopAlg\n";
    LoopAlg doit(step,stop);
    doit.run();
    actcount=vstop.getCount();
    str<<"\n******************* TRGN Summary *******************\n";
    if (nostep && maxcount > 0) { 
      str<<"TRGN: no update from GN step at TR step "<<actcount<<"\n";
      str<<"  possibly gradient already below threshhold\n";
    }
    if (actcount > 0) {
      str<<"iteration count            = "<<actcount<<"\n";
          str<<"initial objective          = "<<jval0<<"\n";
          str<<"final objective            = "<<jval<<"\n";
          str<<"objective redn             = "<<jval/jval0<<"\n";
          str<<"initial gradient norm      = "<<agnrm0<<"\n";
          str<<"gradient norm              = "<<agnrm<<"\n";
          str<<"gradient redn              = "<<rgnrm<<"\n";
      str<<"initial trust radius       = "<<delta0<<"\n";
      str<<"final trust radius         = "<<step.getDelta()<<"\n";
    }     
    else {
      if (!nostep) {
        str<<"TRGN no initial step due to internal error\n";
      }
    }
      }
      catch (RVLException & e) {
    e<<"\ncalled from TRGNAlg::run\n";
    throw e;
      }
    }
    
    int getCount() const { return actcount; }
    
    TRGNAlg(Operator<Scalar> const & op, 
        Vector<Scalar> & x,
        int _maxcount,
        atype _jtol,
        atype _agtol,
        atype _rgtol,
        atype _eta1,
        atype _eta2,
        atype _gamma1,
        atype _gamma2,
        atype _minstep,
        ostream & _str)
      : Policy(), opeval(op,x),
    jtol(_jtol),agtol(_agtol),rgtol(_rgtol),
    eta1(_eta1),eta2(_eta2),gamma1(_gamma1),gamma2(_gamma2),
        maxcount(_maxcount),minstep(_minstep),
    names(6),nums(6),tols(6),actcount(0),
    str(_str) { this->inittable(); }
    
    TRGNAlg(Operator<Scalar> const & op, 
        Vector<Scalar> & x,
        Table & t,
        ostream & _str)
      : opeval(op,x),
    jtol(getValueFromTable<atype>(t,"ResidualTol")),
    agtol(getValueFromTable<atype>(t,"AbsGradTol")),
    rgtol(getValueFromTable<atype>(t,"RelGradTol")),
    eta1(getValueFromTable<Scalar>(t,"MinDecrease")),
    eta2(getValueFromTable<Scalar>(t,"GoodDecrease")),
    gamma1(getValueFromTable<Scalar>(t,"StepDecrFactor")),
    gamma2(getValueFromTable<Scalar>(t,"StepIncrFactor")),
    maxcount(getValueFromTable<int>(t,"MaxItn")), 
    minstep(getValueFromTable<int>(t,"MinStepTol")),
    names(6),nums(6),tols(6),actcount(0),
    str(_str) { 
      this->inittable(); }

  private:
    OperatorEvaluation<Scalar> opeval;
    mutable atype predred;
    mutable atype actred;
    mutable atype jval;
    atype jtol;
    mutable atype agnrm;
    mutable atype rgnrm;
    atype agtol;
    atype rgtol;
    atype eta1;
    atype eta2;
    atype gamma1;
    atype gamma2;
    int maxcount;
    atype minstep;
    vector<string> names;
    vector<atype *> nums;
    vector<atype> tols;
    mutable int actcount;
    ostream & str;

    void inittable() {
      names[0]="LS Objective"; nums[0]=&jval; tols[0]=jtol;
      names[1]="Gradient Norm"; nums[1]=&agnrm; tols[1]=agtol;
      names[2]="Rel Grad Norm"; nums[2]=&rgnrm; tols[2]=rgtol;
      names[3]="Actual Redn"; nums[3]=&actred; tols[3]=-numeric_limits<atype>::max();
      names[4]="Predicted Redn"; nums[4]=&predred; tols[4]=ScalarFieldTraits<atype>::Zero();
      names[5]="Trust Radius"; nums[5]=&(this->Delta); tols[5]=ScalarFieldTraits<atype>::Zero();
    }
  };
}
#endif

---