optimizerBFGS.h

Go to the documentation of this file.

#ifndef _DLR_OPTIMIZERBFGS_H_
#define _DLR_OPTIMIZERBFGS_H_

#include <vector>
#include <dlrCommon/triple.h>
#include <dlrOptimization/optimizer.h>
#include <dlrOptimization/optimizerLineSearch.h>

namespace dlr {

  namespace optimization {

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

      OptimizerBFGS();

      explicit OptimizerBFGS(const Functor& functor);

      OptimizerBFGS(const OptimizerBFGS& source);

      virtual
      ~OptimizerBFGS();

      virtual std::vector<size_t>
      getNumberOfFunctionCalls() {return this->m_functionCallCount;}

      virtual std::vector<size_t>
      getNumberOfGradientCalls() {return this->m_gradientCallCount;}
    
      virtual std::vector<size_t>
      getNumberOfIterations() {return this->m_iterationCount;}
                  
      virtual void
      setParameters(size_t iterationLimit = 500,
                    size_t numberOfRestarts = 1,
                    double argumentTolerance = 1.2E-7, // Generally 4*(num. Eps.)
                    double gradientTolerance = 0.00001,
                    double lineSearchAlpha = 1.0E-4,
                    double lineSearchArgumentTolerance = 1.0e-7,
                    double numericEpsilon = 3.0E-8,
                    double maximumStepMagnitudeFactor = 100.0);
    
    
      virtual void
      setIterationLimit(size_t iterationLimit) {
        this->m_iterationLimit = iterationLimit;
      }

    
      virtual void
      setNumberOfRestarts(size_t numberOfRestarts) {
        this->m_numberOfRestarts = numberOfRestarts;
      }

    
      virtual void
      setStartPoint(const typename Functor::argument_type& startPoint);


      virtual void
      setVerbosity(int verbosity) {}

    
      virtual OptimizerBFGS&
      operator=(const OptimizerBFGS& source);
    
    protected:

      double
      gradientConvergenceMetric(const argument_type& theta,
                                const result_type& value,
                                const argument_type& gradient);

      virtual
      std::pair<typename Functor::argument_type, typename Functor::result_type>
      run();

      Triple<typename Functor::argument_type, typename Functor::result_type,
             typename Functor::argument_type>
      doBfgs(const argument_type& theta,
             const result_type& startValue,
             const argument_type& startGradient,
             size_t& numberOfFunctionCalls,
             size_t& numberOfGradientCalls,
             size_t& numberOfIterations);

      // Data members.
      double m_argumentTolerance;
      size_t m_iterationLimit;
      double m_gradientTolerance;
      double m_lineSearchAlpha;
      double m_lineSearchArgumentTolerance;
      double m_maximumStepMagnitudeFactor;
      size_t m_numberOfRestarts;
      double m_numericEpsilon;
      OptimizerLineSearch<Functor> m_optimizerLineSearch;
      argument_type m_startPoint;

      // Data members used for bookkeeping.
      std::vector<size_t> m_functionCallCount;
      std::vector<size_t> m_gradientCallCount;
      std::vector<size_t> m_iterationCount;
    
    }; // class OptimizerBFGS

  } // namespace optimization

} // namespace dlr


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

namespace dlr {

  using optimization::OptimizerBFGS;

} // namespace dlr


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

#include <math.h>
#include <dlrNumeric/array1D.h>
#include <dlrNumeric/array2D.h>
#include <dlrNumeric/utilities.h>
#include <dlrOptimization/optimizerCommon.h>

namespace dlr {

  namespace optimization {
  
    template <class Functor>
    OptimizerBFGS<Functor>::
    OptimizerBFGS()
      : Optimizer<Functor>(),
        m_argumentTolerance(),
        m_iterationLimit(),
        m_gradientTolerance(),
        m_lineSearchAlpha(),
        m_lineSearchArgumentTolerance(),
        m_maximumStepMagnitudeFactor(),
        m_numberOfRestarts(),
        m_numericEpsilon(),
        m_optimizerLineSearch(),
        m_startPoint(),
        m_functionCallCount(),
        m_gradientCallCount(),
        m_iterationCount()
    {
      this->setParameters();
    }

    template <class Functor>
    OptimizerBFGS<Functor>::
    OptimizerBFGS(const Functor& functor)
      : Optimizer<Functor>(functor),
        m_argumentTolerance(),
        m_iterationLimit(),
        m_gradientTolerance(),
        m_lineSearchAlpha(),
        m_lineSearchArgumentTolerance(),
        m_maximumStepMagnitudeFactor(),
        m_numberOfRestarts(),
        m_numericEpsilon(),
        m_optimizerLineSearch(functor),
        m_startPoint(),
        m_functionCallCount(),
        m_gradientCallCount(),
        m_iterationCount()    
    {
      this->setParameters();
    }
  
    template<class Functor>
    OptimizerBFGS<Functor>::
    OptimizerBFGS(const OptimizerBFGS& source)
      : Optimizer<Functor>(source),
        m_argumentTolerance(source.m_argumentTolerance),
        m_iterationLimit(source.m_iterationLimit),
        m_gradientTolerance(source.m_gradientTolerance),
        m_lineSearchAlpha(source.m_lineSearchAlpha),
        m_lineSearchArgumentTolerance(source.m_lineSearchArgumentTolerance),
        m_maximumStepMagnitudeFactor(source.m_maximumStepMagnitudeFactor),
        m_numberOfRestarts(source.m_numberOfRestarts),
        m_numericEpsilon(source.m_numericEpsilon),
        m_optimizerLineSearch(source.m_optimizerLineSearch),
        m_startPoint(source.m_startPoint.size()),
        m_functionCallCount(source.m_functionCallCount),
        m_gradientCallCount(source.m_gradientCallCount),
        m_iterationCount(source.m_iterationCount)
    {
      copyArgumentType(source.m_startPoint, this->m_startPoint);
    }

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

    template<class Functor>
    void
    OptimizerBFGS<Functor>::
    setParameters(size_t iterationLimit,
                  size_t numberOfRestarts,
                  double argumentTolerance,
                  double gradientTolerance,
                  double lineSearchAlpha,
                  double lineSearchArgumentTolerance,
                  double numericEpsilon,
                  double maximumStepMagnitudeFactor)
    {
      // Copy input arguments.
      this->m_iterationLimit = iterationLimit;
      this->m_numberOfRestarts = numberOfRestarts;
      this->m_argumentTolerance = argumentTolerance;
      this->m_gradientTolerance = gradientTolerance;
      this->m_lineSearchAlpha = lineSearchAlpha;
      this->m_lineSearchArgumentTolerance = lineSearchArgumentTolerance;
      this->m_numericEpsilon = numericEpsilon;
      this->m_maximumStepMagnitudeFactor = maximumStepMagnitudeFactor;

      // Reset memory of previous minimization.
      this->m_functionCallCount.clear();
      this->m_gradientCallCount.clear();
      this->m_iterationCount.clear();
    
      // We've changed the parameters, so we'll have to rerun the 
      // optimization.  Indicate this by setting the inherited member
      // m_needsOptimization.    
      Optimizer<Functor>::m_needsOptimization = true;
    }    

    template<class Functor>
    void
    OptimizerBFGS<Functor>::
    setStartPoint(const typename Functor::argument_type& startPoint)
    {
      copyArgumentType(startPoint, this->m_startPoint);

      // Reset memory of previous minimization.
      this->m_functionCallCount.clear();
      this->m_gradientCallCount.clear();
      this->m_iterationCount.clear();

      // We've changed the parameters, so we'll have to rerun the 
      // optimization.  Indicate this by setting the inherited member
      // m_needsOptimization.    
      Optimizer<Functor>::m_needsOptimization = true;
    }

    template<class Functor>
    OptimizerBFGS<Functor>&
    OptimizerBFGS<Functor>::
    operator=(const OptimizerBFGS<Functor>& source)
    {
      Optimizer<Functor>::operator=(source);
      this->m_argumentTolerance = source.m_argumentTolerance;
      this->m_iterationLimit = source.m_iterationLimit;
      this->m_gradientTolerance = source.m_gradientTolerance;
      this->m_lineSearchAlpha = source.m_lineSearchAlpha;
      this->m_lineSearchArgumentTolerance = source.m_lineSearchArgumentTolerance;
      this->m_maximumStepMagnitudeFactor = source.m_maximumStepMagnitudeFactor;
      this->m_numberOfRestarts = source.m_numberOfRestarts;
      this->m_numericEpsilon = source.m_numericEpsilon;
      this->m_optimizerLineSearch = source.m_optimizerLineSearch;
      copyArgumentType(source.m_startPoint, this->m_startPoint);

      this->m_functionCallCount = source.m_functionCallCount;
      this->m_gradientCallCount = source.m_gradientCallCount;
      this->m_iterationCount = source.m_iterationCount;
    
      return *this;
    }


    // =============== Protected member functions below =============== //

    template <class Functor>
    double
    OptimizerBFGS<Functor>::
    gradientConvergenceMetric(const argument_type& theta,
                              const result_type& value,
                              const argument_type& gradient)
    {
      double returnValue = 0.0;
      double denominator = std::max(static_cast<double>(value), 1.0);
      for(size_t index = 0; index < theta.size(); ++index) {
        double thetaAbsValue = std::fabs(theta[index]);
        double gradientAbsValue = std::fabs(gradient[index]);
        double candidate =
          gradientAbsValue * std::max(thetaAbsValue, 1.0) / denominator;
        if(candidate > returnValue) {
          returnValue = candidate;
        }
      }
      return returnValue;
    }

    template <class Functor>
    std::pair<typename Functor::argument_type, typename Functor::result_type>
    OptimizerBFGS<Functor>::
    run()
    {
      // Check that we have a valid startPoint.
      if(this->m_startPoint.size() == 0) {
        DLR_THROW3(StateException, "OptimizerBFGS<Functor>::run()",
                   "startPoint has not been initialized.");
      }
    
      // Initialize working location so that we start at the right place.
      argument_type theta(this->m_startPoint.size());
      copyArgumentType(this->m_startPoint, theta);

      // Compute initial values of function and its gradient.
      result_type startValue = this->m_functor(theta);
      argument_type startGradient = this->m_functor.gradient(theta);

      // Now run the optimization.
      this->m_functionCallCount.clear();
      this->m_gradientCallCount.clear();
      this->m_iterationCount.clear();
      size_t functionCallCount;
      size_t gradientCallCount;
      size_t iterationCount;
      Triple<argument_type, result_type, argument_type>
        optimum_optimalValue_gradient =
        this->doBfgs(theta, startValue, startGradient, functionCallCount,
                     gradientCallCount, iterationCount);
      // When accounting function calls, don't forget to add the initial
      // function and gradient evaluations above.
      this->m_functionCallCount.push_back(functionCallCount + 1);
      this->m_gradientCallCount.push_back(gradientCallCount + 1);
      this->m_iterationCount.push_back(iterationCount);

      // Restart as many times as requested.
      for(size_t index = 0; index < this->m_numberOfRestarts; ++index) {
        copyArgumentType(optimum_optimalValue_gradient.first, theta);
        startValue = optimum_optimalValue_gradient.second;
        copyArgumentType(optimum_optimalValue_gradient.third,
                         startGradient);
        optimum_optimalValue_gradient =
          this->doBfgs(theta, startValue, startGradient, functionCallCount,
                       gradientCallCount, iterationCount);
        this->m_functionCallCount.push_back(functionCallCount);
        this->m_gradientCallCount.push_back(gradientCallCount);
        this->m_iterationCount.push_back(iterationCount);
      }

      return std::make_pair(optimum_optimalValue_gradient.first,
                            optimum_optimalValue_gradient.second);
    }
    
    template <class Functor>
    Triple<typename Functor::argument_type, typename Functor::result_type,
           typename Functor::argument_type>
    OptimizerBFGS<Functor>::
    doBfgs(const argument_type& theta,
           const result_type& startValue,
           const argument_type& startGradient,
           size_t& numberOfFunctionCalls,
           size_t& numberOfGradientCalls,
           size_t& numberOfIterations)
    {
      // Basic initializations.
      size_t dimensionality = theta.size();
      numberOfFunctionCalls = 0;
      numberOfGradientCalls = 0;
      numberOfIterations = 0;

      // We'll need non-const versions of the input arguments.
      argument_type thetaLocal(theta.size());
      copyArgumentType(theta, thetaLocal);
      result_type currentValue = startValue;
      // Of course, we should only copy startGradient if it's actually
      // been initialized.
      argument_type currentGradient;
      if(startGradient.size() != 0) {
        copyArgumentType(startGradient, currentGradient);
      } else {
        currentGradient = this->m_functor.gradient(thetaLocal);
        ++numberOfGradientCalls;
      }

      // If gradient magnitude is zero, we're already at an extremum or
      // a saddle point.  In either case, we don't know which way to go.
      // Instead, just terminate.
      if(dotArgumentType(currentGradient, currentGradient) == 0) {
        return makeTriple(thetaLocal, currentValue, currentGradient);
      }
    
      // Check that gradient dimension is correct.
      if(currentGradient.size() != dimensionality) {
        std::ostringstream message;
        message << "startPoint has dimensionality " << dimensionality
                << " but objective function returns gradient with "
                << "dimensionality " << currentGradient.size() << ".";
        DLR_THROW3(RunTimeException, "OptimizerBFGS<Functor>::doBfgs()",
                   message.str().c_str());
      }
    
      // Set up inverse hessian estimate.
      Array2D<double> inverseHessian(dimensionality, dimensionality);
      inverseHessian = 0.0;
      for(size_t index = 0; index < dimensionality; ++index) {
        // initialize to identity.
        inverseHessian(index, index) = 1.0;
      }

      // Set up initial search step.
      argument_type searchStep(dimensionality);
      for(size_t index = 0; index < dimensionality; ++index) {
        searchStep[index] = -(currentGradient[index]);
      }

      // Compute maximum allowable step size, allowing for numerical issues.
      double thetaMagnitudeSquared = 0.0;
      for(size_t index = 0; index < dimensionality; ++index) {
        thetaMagnitudeSquared += thetaLocal[index] * thetaLocal[index];
      }
      double thetaMagnitude = std::sqrt(thetaMagnitudeSquared);
      double maximumStepMagnitude =
        (this->m_maximumStepMagnitudeFactor
         * std::max(thetaMagnitude, static_cast<double>(dimensionality)));

      // Now set line search parameters
      // Make sure line search optimizer has the right objective function.
      // Since optimizer::setObjectiveFunction() doesn't update
      // m_optimizerLineSearch.
      this->m_optimizerLineSearch.setObjectiveFunction(this->m_functor);
      // Set line search parameters.
      this->m_optimizerLineSearch.setParameters(
        this->m_lineSearchArgumentTolerance, this->m_lineSearchAlpha,
        maximumStepMagnitude);
    
      // Perform the bfgs iteration.
      while(1) {
        // Perform line search.
        this->m_optimizerLineSearch.setStartPoint(thetaLocal, currentValue,
                                                  currentGradient);
        this->m_optimizerLineSearch.setInitialStep(searchStep);
        argument_type thetaNew = this->m_optimizerLineSearch.optimum();
        currentValue = this->m_optimizerLineSearch.optimalValue();
        numberOfFunctionCalls +=
          this->m_optimizerLineSearch.getNumberOfFunctionCalls();

        // Revise our initial step to match the data, and update current
        // position in parameter space.
        for(size_t index = 0; index < dimensionality; ++index) {
          searchStep[index] = thetaNew[index] - thetaLocal[index];
          thetaLocal[index] = thetaNew[index];
        }

        // Test for "insufficient parameter change" convergence.
        if(contextSensitiveScale(searchStep, thetaLocal)
           < this->m_argumentTolerance) {
          // Gradient at thetaLocal has not been computed yet.  Return
          // empty gradient.
          return makeTriple(thetaLocal, currentValue, argument_type());
        }

        // Temporarily save gradient value.
        Array1D<double> deltaGradient(dimensionality);
        for(size_t index = 0; index < dimensionality; ++index) {
          deltaGradient[index] = -currentGradient[index];
        }

        // Update gradient.
        currentGradient = this->m_functor.gradient(thetaLocal);
        ++numberOfGradientCalls;

        // Sanity check
        if(currentGradient.size() != dimensionality) {
          std::ostringstream message;
          message << "dimensionality of gradient changed mid-stream from "
                  << dimensionality << " to " << currentGradient.size() << ".";
          DLR_THROW3(RunTimeException, "OptimizerBFGS<Functor>::doBfgs()",
                     message.str().c_str());
        }

        // And compute change in gradient.  Remember that deltaGradient was
        // previously set to -currentGradient.
        for(size_t index = 0; index < dimensionality; ++index) {
          deltaGradient[index] += currentGradient[index];
        }

        // Test for "small gradient" convergence.
        if(this->gradientConvergenceMetric(thetaLocal, currentValue,
                                           currentGradient)
           < this->m_gradientTolerance) {
          return makeTriple(thetaLocal, currentValue, currentGradient);
        }

        // Now prepare to update estimate of the inverse hessian.
        Array1D<double> inverseHessianTimesDeltaGradient =
          matrixMultiply(inverseHessian, deltaGradient);
        double fac = dotArgumentType(deltaGradient, searchStep);
        if(fac < 0.0) {
          fac = 0.0;
        }
        double fae = dotArgumentType(
          deltaGradient, inverseHessianTimesDeltaGradient);
        double mag2DeltaGradient = dotArgumentType(deltaGradient, deltaGradient);
        double mag2SearchStep = dotArgumentType(searchStep, searchStep);

        // Are gradient change and search direction sufficiently aligned?
        if((fac * fac)
           > (this->m_numericEpsilon * mag2DeltaGradient * mag2SearchStep)) {
          // Some useful scalars.
          double oneOverFac = 1.0 / fac;
          double fad = 1.0 / fae;

          // Use deltaGradient as temporary storage
          for(size_t index = 0; index < dimensionality; ++index) {
            deltaGradient[index] =
              (oneOverFac * searchStep[index]
               - fad * inverseHessianTimesDeltaGradient[index]);
          }

          // Now update inverse Hession estimate.
          for(size_t row = 0; row < inverseHessian.rows(); ++row) {
            for(size_t column = 0; column < inverseHessian.columns(); ++column) {
              // First DFP term.
              // 
              // Note(xxx): redundant calculation?  this multiplication
              // is done immediately above.
              double increment0 = (oneOverFac * searchStep[row]
                                   * searchStep[column]);
              // Second DFP term.
              double decrement1 = (fad * inverseHessianTimesDeltaGradient[row]
                                   * inverseHessianTimesDeltaGradient[column]);
              // BFGS term.  Remember that deltaGradient has been preempted for
              // temporary storage at this point.
              double increment2 = (fae * deltaGradient[row]
                                   * deltaGradient[column]);
              // Now do the update.
              inverseHessian(row, column) += (increment0 - decrement1
                                              + increment2);
            }
          }

          // Array2D<double> increment0 =
          //   outerProduct(oneOverFac * searchStep, searchStep);

          // Second DFP term.
          // Array2D<double> decrement1 =
          //   outerProduct(fad * inverseHessianTimesDeltaGradient,
          //                inverseHessianTimesDeltaGradient);

          // BFGS term.  Remember that deltaGradient has been preempted for
          // temporary storage at this point.
          // Array2D<double> increment2 =
          //   outerProduct(fae * deltaGradient, deltaGradient);

          // Actually do the update.
          // inverseHessian += increment0;
          // inverseHessian -= decrement1;
          // inverseHessian += increment2;
        }

        // Use Newton's method to choose the next update.
        matrixMultiplyArgumentType(inverseHessian, currentGradient, searchStep);
        for(size_t index = 0; index < dimensionality; ++index) {
          searchStep[index] = -(searchStep[index]);
        }

        // Check for convergence failure.
        ++numberOfIterations;
        if(numberOfIterations > this->m_iterationLimit) {
          // We're going to bail out, but save the result anyway, just
          // in case someone cares.
          this->setOptimum(thetaLocal, currentValue, false);

          // Now throw the exception.
          std::ostringstream message;
          message << "Iteration limit of " << this->m_iterationLimit
                  << " exceeded.";
          DLR_THROW3(RunTimeException, "OptimizerBFGS<Functor>::doBfgs()",
                     message.str().c_str());
        }
      }
    }

  } // namespace optimization

} // namespace dlr

#endif // #ifndef _DLR_OPTIMIZERBFGS_H_

---