optimizerLineSearch.h

Go to the documentation of this file.

#ifndef _DLR_OPTIMIZERLINESEARCH_H_
#define _DLR_OPTIMIZERLINESEARCH_H_

#include <dlrOptimization/optimizer.h>
#include <dlrOptimization/optimizerCommon.h>

namespace dlr {

  namespace optimization {

    template <class Functor>
    class OptimizerLineSearch
      : public Optimizer<Functor>
    {
    public:
      // Typedefs for convenience
      typedef typename Functor::argument_type argument_type;
      typedef typename Functor::result_type result_type;

      OptimizerLineSearch();

      explicit OptimizerLineSearch(const Functor& functor);

      OptimizerLineSearch(const Functor& functor, 
                          const argument_type& startPoint,
                          const argument_type& startGradient,
                          const result_type& startValue);

      OptimizerLineSearch(const OptimizerLineSearch& source);

      virtual
      ~OptimizerLineSearch();

      size_t
      getNumberOfFunctionCalls() {return this->m_functionCallCount;}
    
      virtual void
      setStartPoint(const argument_type& startPoint);

      virtual void
      setStartPoint(const argument_type& startPoint,
                    const result_type& startValue,
                    const argument_type& startGradient);

      virtual void
      setInitialStep(const argument_type& initialStep);

      void
      setParameters(double argumentTolerance=1.0e-7,
                    double alpha=1.0e-4,
                    double maximumStepMagnitude=100.0);
    
      OptimizerLineSearch&
      operator=(const OptimizerLineSearch& source);
    
    protected:

      void
      checkState() {
        // Warning(xxx): currently no state checking is implemented.
      }
    
      std::pair<typename Functor::argument_type, typename Functor::result_type>
      run();
    
      // Data members
      double m_alpha;
      double m_argumentTolerance;
      size_t m_functionCallCount;
      argument_type m_initialStep;
      double m_initialStepMagnitude;
      double m_maximumStepMagnitude;
      argument_type m_startGradient;
      argument_type m_startPoint;
      result_type m_startValue;

    }; // class OptimizerLineSearch

  } // namespace optimization

} // namespace dlr


/* ======= Declarations to maintain compatibility with legacy code. ======= */

namespace dlr {

  using optimization::OptimizerLineSearch;

} // namespace dlr


/*******************************************************************
 * Member function definitions follow.  This would be a .cpp file
 * if it weren't templated.
 *******************************************************************/

#include <cmath>
#include <sstream>
#include <dlrCommon/exception.h>
#include <dlrOptimization/optimizerCommon.h>

namespace dlr {

  namespace optimization {

    template <class Functor>
    OptimizerLineSearch<Functor>::
    OptimizerLineSearch()
      : Optimizer<Functor>(),
        m_alpha(0.0),
        m_argumentTolerance(0.0),
        m_functionCallCount(0),
        m_initialStep(),
        m_initialStepMagnitude(0.0),
        m_maximumStepMagnitude(0.0),
        m_startGradient(),
        m_startPoint(),
        m_startValue(0.0)
    {
      this->setParameters();
    }

    template <class Functor>
    OptimizerLineSearch<Functor>::
    OptimizerLineSearch(const Functor& functor)
      : Optimizer<Functor>(functor),
        m_alpha(0.0),
        m_argumentTolerance(0.0),
        m_functionCallCount(0),
        m_initialStep(),
        m_initialStepMagnitude(0.0),
        m_maximumStepMagnitude(0.0),
        m_startGradient(),
        m_startPoint(),
        m_startValue(0.0)
    {
      this->setParameters();
    }
  
    template<class Functor>
    OptimizerLineSearch<Functor>::
    OptimizerLineSearch(const OptimizerLineSearch& source)
      : Optimizer<Functor>(source),
        m_alpha(source.m_alpha),
        m_argumentTolerance(source.argumentTolerance),
        m_functionCallCount(source.m_functionCallCount),
        m_initialStep(source.m_initialStep.size()),
        m_initialStepMagnitude(source.initialStepMagnitude),
        m_maximumStepMagnitude(source.m_maximumStepMagnitude),
        m_startGradient(source.m_startGradient.size()),
        m_startPoint(source.m_startPoint.size()),
        m_startValue(source.m_startValue)
    {
      copyArgumentType(source.m_initialStep, this->m_initialStep);
      copyArgumentType(source.m_startGradient, this->m_startGradient);
      copyArgumentType(source.m_startPoint, this->m_startPoint);
    }

    template <class Functor>
    OptimizerLineSearch<Functor>::
    ~OptimizerLineSearch()
    {
      // Empty
    }

    template <class Functor>
    void OptimizerLineSearch<Functor>::
    setStartPoint(const argument_type& startPoint)
    {
      // Copy input argument.
      copyArgumentType(startPoint, this->m_startPoint);

      // Compute other starting info.  Note that we don't need to use
      // copyArgumentType(), because we have sole ownership of the
      // return value from m_functor.gradient().
      this->m_startValue = this->m_functor(startPoint);
      this->m_startGradient = this->m_functor.gradient(startPoint);
    
      // Reset the count of function calls.  We ignore the calls used 
      // to initialize m_startValue, etc. above.
      this->m_functionCallCount = 0;

      // We've changed the starting point, so we'll have to rerun the 
      // optimization.  Indicate this by setting the inherited member
      // m_needsOptimization.
      this->m_needsOptimization = true;
    }

    template<class Functor>
    void
    OptimizerLineSearch<Functor>::
    setStartPoint(const argument_type& startPoint,
                  const result_type& startValue,
                  const argument_type& startGradient)
    {
      // Copy input arguments.
      copyArgumentType(startPoint, this->m_startPoint);
      this->m_startValue = startValue;
      copyArgumentType(startGradient, this->m_startGradient);
    
      // Reset the count of function calls.
      this->m_functionCallCount = 0;

      // We've changed the starting point, so we'll have to rerun the 
      // optimization.  Indicate this by setting the inherited member
      // m_needsOptimization.
      this->m_needsOptimization = true;
    }

    template<class Functor>
    void OptimizerLineSearch<Functor>::
    setInitialStep(const argument_type& initialStep)
    {
      // First, make sure this is a valid initialStep.  
      double initialStepMagnitude =
        std::sqrt(dotArgumentType(initialStep, initialStep));
      if(initialStepMagnitude == 0.0) {
        DLR_THROW3(ValueException, "OptimizerLineSearch::setInitialStep()",
                   "initialStep magnitude is equal to 0.0.");
      }

      // OK, valid initialStep, save it for later.  Since we'll also
      // be needing the magnitude of initialStep, save it too.
      this->m_initialStepMagnitude = initialStepMagnitude;
      copyArgumentType(initialStep, this->m_initialStep);
    
      // Reset the count of function calls.
      this->m_functionCallCount = 0;

      // We've changed the line direction, so we'll have to rerun the 
      // optimization.  Indicate this by setting the inherited member
      // m_needsOptimization.
      this->m_needsOptimization = true;
    }
  
    template<class Functor>
    void OptimizerLineSearch<Functor>::
    setParameters(double argumentTolerance,
                  double alpha,
                  double maximumStepMagnitude)
    {
      this->m_argumentTolerance = argumentTolerance;
      this->m_alpha = alpha;
      this->m_maximumStepMagnitude = maximumStepMagnitude;

      // Reset the count of function calls.
      this->m_functionCallCount = 0;

      // We've changed the parameters, so we'll have to rerun the 
      // optimization.  Indicate this by setting the inherited member
      // m_needsOptimization.    
      this->m_needsOptimization = true;
    }    
  
    template <class Functor>
    OptimizerLineSearch<Functor>& OptimizerLineSearch<Functor>::
    operator=(const OptimizerLineSearch& source)
    {
      // First, call the assignment operator of the parent class.
      Optimizer<Functor>::operator=(source);

      // Then copy all of the variables which are specific to
      // OptimizerLineSearch.
      this->m_alpha = source.m_alpha;
      this->m_argumentTolerance = source.m_argumentTolerance;
      this->m_functionCallCount = source.m_functionCallCount;
      this->m_initialStepMagnitude = source.m_initialStepMagnitude;
      this->m_maximumStepMagnitude = source.m_maximumStepMagnitude;
      this->m_startValue = source.m_startValue;
      copyArgumentType(source.m_initialStep, this->m_initialStep);
      copyArgumentType(source.m_startGradient, this->m_startGradient);
      copyArgumentType(source.m_startPoint, this->m_startPoint);
      return *this;
    }

    template <class Functor>
    std::pair<typename Functor::argument_type, typename Functor::result_type>
    OptimizerLineSearch<Functor>::
    run()
    {
      // Make sure all required starting values have been specified.
      this->checkState();

      // Reset the count of function calls.
      this->m_functionCallCount = 0;
    
      // Initialize a few variables.
      argument_type initialStep;
      copyArgumentType(this->m_initialStep, initialStep);
      double initialStepMagnitude = this->m_initialStepMagnitude;

      // Restrict length of initial step.
      if(this->m_initialStepMagnitude > this->m_maximumStepMagnitude) {
        initialStepMagnitude = this->m_maximumStepMagnitude;
        for(size_t index = 0; index < initialStep.size(); ++index) {
          initialStep[index] *=
            this->m_maximumStepMagnitude / this->m_initialStepMagnitude;
        }
      }

      // Compute degree of similarity between startGradient and initialStep.
      double slope = dotArgumentType(this->m_startGradient, initialStep);
      if(slope >= 0.0) {
        std::ostringstream message;
        message << "Initial search direction: " << initialStep << " "
                << "is not downhill with respect to initial gradient: "
                << this->m_startGradient << ".";
        DLR_THROW3(StateException, "OptimizerLineSearch::run()",
                   message.str().c_str());
                   
      }

      // Compute smallest allowable step, allowing for numerical issues.
      double scale = contextSensitiveScale(initialStep, this->m_startPoint);
      if(scale == 0.0) {
        DLR_THROW3(RunTimeException, "OptimizerLineSearch::run()",
                   "Invalid initial scale.  "
                   "Perhaps initialStep is the zero vector.");
      }
      double minimumLambda = this->m_argumentTolerance / scale;

      // Now do the search.
      double lambdaValue = 1.0;
      argument_type nextSamplePoint(initialStep.size());
      double previousLambdaValue = 1.0;
      double previousNextFunctionValue = this->m_startValue;
      while(1) {
        // Check first termination criterion.
        // Note: this differs from NRC in that NRC first computes
        // nextFunctionValue before doing this check, and consequently
        // returns it to the calling context if the following conditional
        // passes.
        if(lambdaValue < minimumLambda) {
          // In this case, root finding programs should double check
          // convergence!
          copyArgumentType(this->m_startPoint, nextSamplePoint);
          // Note: This return value differs from NRC, in that we return
          // startValue instead of the most recently computed value.
          // return std::make_pair(nextSamplePoint, nextFunctionValue);
          return std::make_pair(nextSamplePoint, this->m_startValue);
        }

        // Compute next place to evaluate function.
        for(size_t index = 0; index < initialStep.size(); ++index) {
          nextSamplePoint[index] =
            this->m_startPoint[index] + (lambdaValue * initialStep[index]);
        }

        // Do the function evaluation.
        result_type nextFunctionValue = this->m_functor(nextSamplePoint);
        ++(this->m_functionCallCount);

        // Check second termination criterion.
        if(nextFunctionValue
           < this->m_startValue + (this->m_alpha * lambdaValue * slope)) {
          return std::make_pair(nextSamplePoint, nextFunctionValue);
        }

        // Compute next value of lambda.
        double lambdaHat;
        if(lambdaValue == 1.0) {
          lambdaHat =
            (-slope / (2.0 * (nextFunctionValue - this->m_startValue
                              - slope)));
        } else {
          // Buld vector.
          double b0 = (nextFunctionValue - this->m_startValue
                       - (lambdaValue * slope));
          double b1 = (previousNextFunctionValue - this->m_startValue
                       - (previousLambdaValue * slope));
          // Build matrix.
          double deltaLambda = lambdaValue - previousLambdaValue;
          double a00 = 1.0 / (lambdaValue * lambdaValue * deltaLambda);
          double a01 = -1.0 / (previousLambdaValue * previousLambdaValue
                               * deltaLambda);
          double a10 = -previousLambdaValue / (lambdaValue * lambdaValue
                                               * deltaLambda);
          double a11 = lambdaValue / (previousLambdaValue * previousLambdaValue
                                      * deltaLambda);
          // Do matrix multiplication.
          double ab0 = a00 * b0 + a01 * b1;
          double ab1 = a10 * b0 + a11 * b1;

          // Use result to choose next value of lambda.
          if(ab0 == 0.0) {
            if(ab1 == 0.0) {
              DLR_THROW3(RunTimeException, "OptimizerLineSearch::run()",
                         "Invalid value for internal variable ab1.");
            }
            lambdaHat = -slope / (2.0 * ab1);
          } else {
            double discriminant = ab1 * ab1 - (3.0 * ab0 * slope);
            if(discriminant < 0.0) {
              DLR_THROW3(RunTimeException, "OptimizerLineSearch::run()",
                         "Roundoff error, discriminant < 0.0");
            }
            lambdaHat = (-ab1 + std::sqrt(discriminant)) / (3.0 * ab0);
          }
          if(lambdaHat > (0.5 * lambdaValue)) {
            lambdaHat = 0.5 * lambdaValue;
          }
        }
        previousLambdaValue = lambdaValue;
        previousNextFunctionValue = nextFunctionValue;
        lambdaValue = std::max(lambdaHat, 0.1 * lambdaValue);
      }
    }

  } // namespace optimization

} // namespace dlr

#endif // #ifndef _DLR_OPTIMIZERLINESEARCH_H_