smooth_nesterov.hxx

Go to the documentation of this file.

/*
 * smooth_nesterov.hxx
 *
 *  Created on: Dec 23, 2013
 *      Author: bsavchyn
 */

#ifndef SMOOTH_NESTEROV_HXX_
#define SMOOTH_NESTEROV_HXX_
#include <opengm/inference/inference.hxx>
#include <opengm/inference/trws/trws_base.hxx>
#include <opengm/inference/auxiliary/primal_lpbound.hxx>
#include <opengm/inference/trws/smoothing_strategy.hxx>

namespace opengm{

template<class GM>
struct Nesterov_Parameter : public trws_base::SmoothingBasedInference_Parameter<GM>

{
   typedef typename  GM::ValueType ValueType;
   typedef trws_base::DecompositionStorage<GM> Storage;
   typedef typename trws_base::SmoothingBasedInference_Parameter<GM> parent;
   typedef typename parent::SmoothingParametersType SmoothingParametersType;
   typedef typename parent::SumProdSolverParametersType SumProdSolverParametersType;
   typedef typename parent::MaxSumSolverParametersType MaxSumSolverParametersType;
   typedef typename parent::PrimalLPEstimatorParametersType PrimalLPEstimatorParametersType;
   typedef typename parent::SmoothingStrategyType SmoothingStrategyType;
   typedef enum {ADAPTIVE_STEP,WC_STEP,JOJIC_STEP} GradientStepType;

   static GradientStepType getGradientStepType(const std::string& name)
   {
         if (name.compare("WC_STEP")==0) return WC_STEP;
         else if (name.compare("JOJIC_STEP")==0)  return JOJIC_STEP;
         else return ADAPTIVE_STEP;
   }

   static std::string getString(GradientStepType steptype)
   {
      switch (steptype)
      {
      case ADAPTIVE_STEP: return std::string("ADAPTIVE_STEP");
      case WC_STEP   : return std::string("WC_STEP");
      case JOJIC_STEP   : return std::string("JOJIC_STEP");
      default: return std::string("UNKNOWN");
      }
   }

   Nesterov_Parameter(size_t numOfExternalIterations=0,
             ValueType precision=1.0,
             bool absolutePrecision=true,
             size_t numOfInternalIterations=3,
             typename Storage::StructureType decompositionType=Storage::GENERALSTRUCTURE,
             ValueType smoothingGapRatio=4,
             ValueType startSmoothingValue=0.0,
             ValueType primalBoundPrecision=std::numeric_limits<ValueType>::epsilon(),
             size_t maxPrimalBoundIterationNumber=100,
             size_t presolveMaxIterNumber=100,
             bool canonicalNormalization=true,
             ValueType presolveMinRelativeDualImprovement=0.01,
             bool lazyLPPrimalBoundComputation=true,
             ValueType smoothingDecayMultiplier=-1.0,
             SmoothingStrategyType smoothingStrategy=SmoothingParametersType::ADAPTIVE_DIMINISHING,
             bool fastComputations=true,
             bool verbose=false,
             GradientStepType gradientStep=ADAPTIVE_STEP,
             ValueType gamma0=1e6,
             bool plainGradient=false
             )
    :parent(numOfExternalIterations,
          precision,
          absolutePrecision,
          numOfInternalIterations,
          decompositionType,
          smoothingGapRatio,
          startSmoothingValue,
          primalBoundPrecision,
          maxPrimalBoundIterationNumber,
          presolveMaxIterNumber,
          canonicalNormalization,
          presolveMinRelativeDualImprovement,
          lazyLPPrimalBoundComputation,
          smoothingDecayMultiplier,
          smoothingStrategy,
          fastComputations,
          verbose
          ),
          gradientStep_(gradientStep),
          gamma0_(gamma0),
          plainGradient_(plainGradient)
     {};

   GradientStepType gradientStep_;
   ValueType gamma0_;
   bool plainGradient_;

#ifdef TRWS_DEBUG_OUTPUT
     void print(std::ostream& fout)const
     {
        parent::print(fout);
        fout << "gradientStep_="<<getString(gradientStep_)<<std::endl;
        fout << "gamma0_=" << gamma0_ <<std::endl;
        fout << "plainGradient_=" << plainGradient_ <<std::endl;
     }
#endif

};


template<class GM, class ACC>
class NesterovAcceleratedGradient : public trws_base::SmoothingBasedInference<GM, ACC> //Inference<GM, ACC>
{
public:
   typedef trws_base::SmoothingBasedInference<GM, ACC> parent;
   typedef ACC AccumulationType;
   typedef GM GraphicalModelType;
   OPENGM_GM_TYPE_TYPEDEFS;

   //typedef std::vector<typename GM::ValueType> DDVectorType;

   typedef typename parent::Storage Storage;
   typedef typename Storage::DDVectorType DDVectorType;
   typedef typename parent::SumProdSolver SumProdSolver;
   typedef typename parent::MaxSumSolver MaxSumSolver;
   typedef typename parent::PrimalBoundEstimator PrimalBoundEstimator;

   typedef Nesterov_Parameter<GM> Parameter;

// typedef visitors::ExplicitVerboseVisitor<NesterovAcceleratedGradient<GM, ACC> > VerboseVisitorType;
// typedef visitors::ExplicitTimingVisitor <NesterovAcceleratedGradient<GM, ACC> > TimingVisitorType;
// typedef visitors::ExplicitEmptyVisitor  <NesterovAcceleratedGradient<GM, ACC> > EmptyVisitorType;

   typedef visitors::VerboseVisitor<NesterovAcceleratedGradient<GM, ACC> > VerboseVisitorType;
   typedef visitors::TimingVisitor <NesterovAcceleratedGradient<GM, ACC> > TimingVisitorType;
   typedef visitors::EmptyVisitor  <NesterovAcceleratedGradient<GM, ACC> > EmptyVisitorType;


   NesterovAcceleratedGradient(const GraphicalModelType& gm,const Parameter& param
#ifdef TRWS_DEBUG_OUTPUT
         ,std::ostream& fout=std::cout
#endif
   )
   :
      parent(gm,param
#ifdef TRWS_DEBUG_OUTPUT
            ,(param.verbose_ ? fout : *OUT::nullstream::Instance())
#endif
      ),
      _parameters(param),
      _currentDualVector(_getDualVectorSize(),0.0)
   {
#ifdef TRWS_DEBUG_OUTPUT
      parent::_fout << "Parameters of the "<< name() <<" algorithm:"<<std::endl;
      param.print(parent::_fout);
#endif

      if (param.numOfExternalIterations_==0) throw std::runtime_error("NEST: a strictly positive number of iterations must be provided!");
   };

   template<class VISITOR>
   InferenceTermination infer(VISITOR & visitor);

   std::string name() const{ return "NEST"; }
   InferenceTermination infer(){EmptyVisitorType visitor; return infer(visitor);};

// InferenceTermination marginal(const IndexType varID, IndependentFactorType& out) //const
// {
//      parent::_marginalsTemp.resize(parent::_storage.numberOfLabels(varID));
//      parent::_sumprodsolver.GetMarginals(varID, parent::_marginalsTemp.begin());
//      OPENGM_ASSERT(parent::_marginalsTemp.size() == out.size());
//      out.assign(parent::_storage.masterModel(), &varID, &varID+1, ACC::template neutral<ValueType>());
//      for (LabelType i=0;i<out.size();++i)
//         out(i)=parent::_marginalsTemp[i];
// }

private:
   ValueType _evaluateGradient(const DDVectorType& point,DDVectorType* pgradient);
   ValueType _evaluateSmoothObjective(const DDVectorType& point);
   size_t    _getDualVectorSize()const{return parent::_storage.getDDVectorSize();}
   void      _SetDualVariables(const DDVectorType& lambda);
   ValueType _estimateOmega0()const{return 1;};//TODO: exchange with a reasonable value
   void _InitSmoothing();//TODO: refactor me
   void _GradientStep(const DDVectorType& gradient, const DDVectorType& startPoint, DDVectorType& endPoint, ValueType stepsize);

   ValueType getLipschitzConstant()const;

   Parameter     _parameters;
   DDVectorType     _currentDualVector;
};

template<class GM, class ACC>
typename NesterovAcceleratedGradient<GM,ACC>::ValueType
NesterovAcceleratedGradient<GM,ACC>::getLipschitzConstant()const
{
 ValueType result=0;
 for (IndexType modelId=0;modelId<parent::_storage.numberOfModels();++modelId)
    result+=(ValueType)parent::_storage.size(modelId);

 return result/parent::_sumprodsolver.GetSmoothing();
}

template<class GM, class ACC>
typename NesterovAcceleratedGradient<GM,ACC>::ValueType
NesterovAcceleratedGradient<GM,ACC>::_evaluateGradient(const DDVectorType& point,DDVectorType* pgradient)
{
   ValueType bound=_evaluateSmoothObjective(point);
   parent::_sumprodsolver.GetMarginalsMove();
   std::vector<ValueType> buffer1st;
   std::vector<ValueType> buffer;
   //transform marginals to dual vector
   pgradient->resize(_currentDualVector.size());
   typename DDVectorType::iterator gradientIt=pgradient->begin();
   for (IndexType varId=0;varId<parent::_storage.masterModel().numberOfVariables();++varId)// all variables
   {
      const typename Storage::SubVariableListType& varList=parent::_storage.getSubVariableList(varId);

      if (varList.size()==1) continue;
      typename Storage::SubVariableListType::const_iterator modelIt=varList.begin();
      IndexType firstModelId=modelIt->subModelId_;
      IndexType firstModelVariableId=modelIt->subVariableId_;
      buffer1st.resize(parent::_storage.masterModel().numberOfLabels(varId));
      buffer.resize(parent::_storage.masterModel().numberOfLabels(varId));
      parent::_sumprodsolver.GetMarginalsForSubModel(firstModelId,firstModelVariableId,buffer1st.begin());
      ++modelIt;
      for(;modelIt!=varList.end();++modelIt) //all related models
      {
         parent::_sumprodsolver.GetMarginalsForSubModel(modelIt->subModelId_,modelIt->subVariableId_,buffer.begin());
         gradientIt=std::transform(buffer.begin(),buffer.end(),buffer1st.begin(),gradientIt,std::minus<ValueType>());
      }
   }

   return bound;
}


template<class GM, class ACC>
typename NesterovAcceleratedGradient<GM,ACC>::ValueType
NesterovAcceleratedGradient<GM,ACC>::_evaluateSmoothObjective(const DDVectorType& point)
{
   _SetDualVariables(point);
   parent::_sumprodsolver.ForwardMove();
   return parent::_sumprodsolver.bound();
}

template<class GM, class ACC>
void NesterovAcceleratedGradient<GM,ACC>::_SetDualVariables(const DDVectorType& lambda)
{
   DDVectorType delta(_currentDualVector.size());
   std::transform(lambda.begin(),lambda.end(),_currentDualVector.begin(),delta.begin(),std::minus<ValueType>());
   _currentDualVector=lambda;
   parent::_storage.addDDvector(delta);
};

template<class GM, class ACC>
void NesterovAcceleratedGradient<GM,ACC>::_InitSmoothing()
{
   if (_parameters.smoothing_ > 0.0)
      parent::_sumprodsolver.SetSmoothing(_parameters.smoothing_);
   else
      throw std::runtime_error("NesterovAcceleratedGradient::_InitSmoothing(): Error! Automatic smoothing selection is not implemented yet.");
};

template<class GM, class ACC>
void NesterovAcceleratedGradient<GM,ACC>::_GradientStep(const DDVectorType& gradient, const DDVectorType& startPoint, DDVectorType& endPoint, ValueType stepsize)
{
   std::transform(gradient.begin(),gradient.end(),endPoint.begin(),std::bind1st(std::multiplies<ValueType>(),stepsize));
   std::transform(startPoint.begin(),startPoint.end(),endPoint.begin(),endPoint.begin(),std::plus<ValueType>());
}

template<class GM, class ACC>
template<class VISITOR>
InferenceTermination NesterovAcceleratedGradient<GM,ACC>::infer(VISITOR & vis)
{
   visitors::VisitorWrapper<VISITOR,NesterovAcceleratedGradient<GM, ACC> > visitor(&vis,this);

   size_t counter=0;
   visitor.addLog("primalLPbound");
   visitor.addLog("oracleCalls");
   visitor.begin();

   if (parent::_sumprodsolver.GetSmoothing()<=0.0)
   {

   parent::_maxsumsolver.ForwardMove();
   parent::_maxsumsolver.EstimateIntegerLabelingAndBound();
   parent::_SelectOptimalBoundsAndLabeling();
    ++counter;

   if (parent::_sumprodsolver.CheckDualityGap(parent::value(),parent::bound()))
   {
   #ifdef TRWS_DEBUG_OUTPUT
      parent::_fout << "NesterovAcceleratedGradient::_CheckStoppingCondition(): Precision attained! Problem solved!"<<std::endl;
   #endif

       return NORMAL;
   }

   counter+=parent::_EstimateStartingSmoothing(visitor);

   }else
   {
      parent::_sumprodsolver.SetSmoothing(_parameters.startSmoothingValue());
   }

   DDVectorType gradient(_currentDualVector.size()),
         lambda(_currentDualVector.size()),
         y(_currentDualVector),
         v(_currentDualVector);
   DDVectorType w(_currentDualVector.size());//temp variable
   ValueType   alpha,

   gamma= _parameters.gamma0_,
   omega=_estimateOmega0();

   omega=omega/2.0;

   for (size_t i=0;i<_parameters.maxNumberOfIterations();++i)
   {
#ifdef TRWS_DEBUG_OUTPUT
      parent::_fout <<"i="<<i<<std::endl;
#endif
      //gradient step with approximate linear search:
      ValueType doubledLipschitzConstant=2*getLipschitzConstant();//depends on a smoothing value
//===================== begin of internal loop ===========================================
      for (size_t j=0;j<_parameters.numberOfInternalIterations();++j)
      {
      ValueType mul=1.0;

      counter+=2; ValueType oldObjVal=_evaluateGradient(y,&gradient);
#ifdef TRWS_DEBUG_OUTPUT
      parent::_fout <<"Dual smooth objective ="<<oldObjVal<<std::endl;
#endif

      switch (_parameters.gradientStep_)
      {
      case Parameter::ADAPTIVE_STEP:
      {
      ValueType norm2=std::inner_product(gradient.begin(),gradient.end(),gradient.begin(),(ValueType)0);//squared L2 norm
      ValueType newObjVal;

      omega/=4.0;

      do
      {
         omega*=2.0;
         ACC::iop(-1.0,1.0,mul);
         _GradientStep(gradient,y,lambda,mul/omega);

         //newObjVal=_evaluateSmoothObjective(lambda,((j+1)==_parameters.numberOfInternalIterations()));
         ++counter; newObjVal=_evaluateSmoothObjective(lambda);
      }
      while ( ACC::bop(newObjVal,(ValueType)(oldObjVal+mul*norm2/2.0/omega)) && (omega < doubledLipschitzConstant));//TODO: +/- and >/< depending on ACC
      ++counter; parent::_sumprodsolver.GetMarginalsAndDerivativeMove();

      if (omega >= doubledLipschitzConstant)
      {
#ifdef TRWS_DEBUG_OUTPUT
         parent::_fout << "Step size is smaller then the inverse Lipschitz constant. Passing to smoothing update." <<std::endl;
#endif
      }
      }
      break;

      case Parameter::WC_STEP:
         omega=doubledLipschitzConstant/2.0;
         _GradientStep(gradient,y,lambda,mul/omega);
         break;

      case Parameter::JOJIC_STEP:
         omega=1.0/parent::_sumprodsolver.GetSmoothing();
         _GradientStep(gradient,y,lambda,mul/omega);
         break;
      default:
         std::runtime_error("NesterovAcceleratedGradient::infer():Error! Unknown value of the step size selector _parameters.gradientStep_.");
         break;
      }//switch


      if (!_parameters.plainGradient_)
      {
         //updating parameters
         alpha=(sqrt(gamma*gamma+4*omega*gamma)-gamma)/omega/2.0;
         gamma*=(1-alpha);
         //v+=(alpha/gamma)*gradient;
         trws_base::transform_inplace(gradient.begin(),gradient.end(),std::bind1st(std::multiplies<ValueType>(),mul*alpha/gamma));
         std::transform(v.begin(),v.end(),gradient.begin(),v.begin(),std::plus<ValueType>());

         //y=alpha*v+(1-alpha)*lambda;
         trws_base::transform_inplace(lambda.begin(),lambda.end(),std::bind1st(std::multiplies<ValueType>(),(1-alpha)));
         std::transform(v.begin(),v.end(),w.begin(),std::bind1st(std::multiplies<ValueType>(),alpha));
         std::transform(w.begin(),w.end(),lambda.begin(),y.begin(),std::plus<ValueType>());
      }else //plain gradient algorithm
      {
         std::copy(lambda.begin(),lambda.end(),y.begin());
      }
}

//=================================== end of internal loop ===============================================
         parent::_maxsumsolver.ForwardMove();//initializes a move, makes a forward move and computes the dual bound, is used also in derivative computation in the next line
         parent::_maxsumsolver.EstimateIntegerLabelingAndBound();
         ++counter;
#ifdef TRWS_DEBUG_OUTPUT
         parent::_fout << "_maxsumsolver.bound()=" <<parent::_maxsumsolver.bound()<<", _maxsumsolver.value()=" <<parent::_maxsumsolver.value() <<std::endl;
#endif

         ValueType  derivative=parent::_EstimateRhoDerivative();
#ifdef TRWS_DEBUG_OUTPUT
            parent::_fout << "derivative="<<derivative<<std::endl;
#endif

         InferenceTermination returncode;
         if ( parent::_CheckStoppingCondition(&returncode))
         {
            visitor();
            //std::cout << "counter=" << counter <<std::endl;
            visitor.log("oracleCalls",(double)counter);
            visitor.log("primalLPbound",(double)parent::_bestPrimalLPbound);
            visitor.end();
            return NORMAL;
         }

         size_t flag=visitor();
         //std::cout << "counter=" << counter <<std::endl;
         visitor.log("oracleCalls",(double)counter);
         visitor.log("primalLPbound",(double)parent::_bestPrimalLPbound);
         if( flag != visitors::VisitorReturnFlag::ContinueInf ){
            break;
         }

         parent::_UpdateSmoothing(parent::_bestPrimalBound,parent::_maxsumsolver.bound(),parent::_sumprodsolver.bound(),derivative,i+1);


   }
   //update smoothing
   parent::_SelectOptimalBoundsAndLabeling();
      visitor();
      visitor.log("oracleCalls",(double)counter);
      visitor.log("primalLPbound",(double)parent::_bestPrimalLPbound);
   visitor.end();

   return NORMAL;
}

//template<class GM>
//class NesterovAcceleratedGradient<GM, opengm::Integrator>
//{
//public:
// typedef GM GraphicalModelType;
// OPENGM_GM_TYPE_TYPEDEFS;
// typedef NesterovAcceleratedGradient <GM, opengm::Maximizer> parent;
// typedef typename parent::AccumulationType AccumulationType;
//
// typedef typename parent::Parameter Parameter;
// typedef typename parent::Storage Storage;
// typedef typename Storage::DDVectorType DDVectorType;
// typedef typename parent::SumProdSolver SumProdSolver;
// typedef typename parent::MaxSumSolver MaxSumSolver;
// typedef typename parent::PrimalBoundEstimator PrimalBoundEstimator;
//
// typedef typename parent::VerboseVisitorType VerboseVisitorType;
// typedef typename parent::TimingVisitorType TimingVisitorType;
//    typedef typename parent::EmptyVisitorType EmptyVisitorType;
//
// NesterovAcceleratedGradient(const GraphicalModelType& gm, const Parameter& param
// #ifdef TRWS_DEBUG_OUTPUT
//            ,std::ostream& fout=std::cout
// #endif
//            )
//    {
//       Parameter param1=param;
//       param1.smoothingStrategy()   = Parameter::SmoothingParametersType::FIXED;
//       param1.setStartSmoothingValue(1.0);
//
//       _pparent = new parent(gm,param1
// #ifdef TRWS_DEBUG_OUTPUT
//             ,fout
// #endif
//       );
//    }
// virtual ~NesterovAcceleratedGradient(){delete _pparent;}
//
// std::string name() const{ return "NEST"; }
//
// InferenceTermination marginal(const IndexType varID, IndependentFactorType& out) const
// {
//    return _pparent->marginal(varID,out);
// }
//
// template<class VISITOR>
// InferenceTermination infer(VISITOR & visitor){return _pparent->infer(visitor);};
//
// InferenceTermination infer(){return _pparent->infer();}
//
//private:
// parent* _pparent;
//};


}//namespace opengm

#endif /* SMOOTH_NESTEROV_HXX_ */