optimizerLM.h

Go to the documentation of this file.

#ifndef _DLR_OPTIMIZERLM_H_
#define _DLR_OPTIMIZERLM_H_

#include <vector>
#include <dlrCommon/types.h>
#include <dlrOptimization/optimizer.h>
#include <dlrOptimization/optimizerLineSearch.h>
#include <dlrLinearAlgebra/linearAlgebra.h>
#include <dlrNumeric/utilities.h>

namespace dlr {

  namespace optimization {

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


      OptimizerLM();


      explicit OptimizerLM(const Functor& functor);


      OptimizerLM(const OptimizerLM& source);


      virtual
      ~OptimizerLM();

    
//     /** 
//      * This method returns the number of function calls required to
//      * complete the previous minimization.  If the minimization
//      * parameter "restarts" is 0, there will be only one element in
//      * the returned vector.  If restarts is greater than 0, the first
//      * element of the return value will reflect the number of function
//      * calls in the initial minimization, and subsequent numbers will
//      * reflect the number of function calls in the following restarted
//      * minimizations.  If a valid minimization has not been performed
//      * since the last update to startPoint, parameters, etc., then the
//      * return value will be an empty vector.
//      *
//      * @return a vector of function call counts.
//      */
//     virtual std::vector<size_t>
//     getNumberOfFunctionCalls() {return this->m_functionCallCount;}

    
//     /** 
//      * This method returns the number of gradient calls required to
//      * complete the previous minimization.  If the minimization
//      * parameter "restarts" is 0, there will be only one element in
//      * the returned vector.  If restarts is greater than 0, the first
//      * element of the return value will reflect the number of gradient
//      * calls in the initial minimization, and subsequent numbers will
//      * reflect the number of gradient calls in the following restarted
//      * minimizations.  If a valid minimization has not been performed
//      * since the last update to startPoint, parameters, etc., then the
//      * return value will be an empty vector.
//      *
//      * @return a vector of gradient call counts.
//      */
//     virtual std::vector<size_t>
//     getNumberOfGradientCalls() {return this->m_gradientCallCount;}
    
//     /** 
//      * This method returns the number of iterations required to
//      * complete the previous minimization.  If the minimization
//      * parameter "restarts" is 0, there will be only one element in
//      * the returned vector.  If restarts is greater than 0, the first
//      * element of the return value will reflect the number of
//      * iterations in the initial minimization, and subsequent numbers
//      * will reflect the number of iterations in the following
//      * restarted minimizations.  If a valid minimization has not been
//      * performed since the last update to startPoint, parameters,
//      * etc., then the return value will be an empty vector.
//      *
//      * @return a vector of iteration counts.
//      */
//     virtual std::vector<size_t>
//     getNumberOfIterations() {return this->m_iterationCount;}


      virtual void
      setMinimumGradientMagnitude(double minimumGradientMagnitude);

    
      virtual void
      setParameters(double initialLambda = 1.0,
                    size_t maxIterations = 40,
                    double maxLambda = 1.0E7,
                    double minLambda = 1.0E-13,
                    double minError = 0.0,
                    double minimumGradientMagnitude = 1.0E-5,
                    double minDrop = 1.0E-4,
                    size_t strikes = 3,
                    int maxBackSteps = -1,
                    int verbosity = 0);
    
    
      virtual void
      setStartPoint(const typename Functor::argument_type& startPoint);

    
      virtual void
      setVerbosity(int verbosity) {m_verbosity = verbosity;}

    
      virtual OptimizerLM&
      operator=(const OptimizerLM& 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();


      inline virtual void
      verboseWrite(const char* message, int verbosity);


      template <class Type>
      inline void
      verboseWrite(const char* intro, const Type& subject, int verbosity);
    
      // Data members.
      double m_initialLambda;
      int m_maxBackSteps;
      size_t m_maxIterations;
      double m_maxLambda;
      double m_minDrop;
      double m_minError;
      double m_minGrad;
      double m_minLambda;
      argument_type m_startPoint;
      size_t m_strikes;
      int m_verbosity;

    }; // class OptimizerLM

  } // namespace optimization

} // namespace dlr


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

namespace dlr {

  using optimization::OptimizerLM;

} // 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>
    OptimizerLM<Functor>::
    OptimizerLM()
      : Optimizer<Functor>(),
        m_initialLambda(),
        m_maxBackSteps(),
        m_maxIterations(),
        m_maxLambda(),
        m_minDrop(),
        m_minError(),
        m_minGrad(),
        m_minLambda(),
        m_startPoint(),
        m_strikes(),
        m_verbosity()
    {
      this->setParameters();
    }


    template <class Functor>
    OptimizerLM<Functor>::
    OptimizerLM(const Functor& functor)
      : Optimizer<Functor>(functor),
        m_initialLambda(),
        m_maxBackSteps(),
        m_maxIterations(),
        m_maxLambda(),
        m_minDrop(),
        m_minError(),
        m_minGrad(),
        m_minLambda(),
        m_startPoint(),
        m_strikes(),
        m_verbosity()
    {
      this->setParameters();
    }
  

    template<class Functor>
    OptimizerLM<Functor>::
    OptimizerLM(const OptimizerLM& source)
      : Optimizer<Functor>(source),
        m_initialLambda(source.m_initialLambda),
        m_maxBackSteps(source.m_maxBackSteps),
        m_maxIterations(source.m_maxIterations),
        m_maxLambda(source.m_maxLambda),
        m_minDrop(source.m_minDrop),
        m_minError(source.m_minError),
        m_minGrad(source.m_minGrad),
        m_minLambda(source.m_minLambda),
        m_startPoint(source.m_startPoint.size()),
        m_strikes(source.m_strikes),
        m_verbosity(source.m_verbosity)
    {
      copyArgumentType(source.m_startPoint, this->m_startPoint);
    }

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


    // This method sets one of the termination criteria of the
    // optimization.
    template<class Functor>
    void
    OptimizerLM<Functor>::
    setMinimumGradientMagnitude(double minimumGradientMagnitude)
    {
      this->m_minGrad = minimumGradientMagnitude * minimumGradientMagnitude;
    }

  
    template<class Functor>
    void
    OptimizerLM<Functor>::
    setParameters(double initialLambda,
                  size_t maxIterations,
                  double maxLambda,
                  double minLambda,
                  double minError,
                  double minimumGradientMagnitude,
                  double minDrop,
                  size_t strikes,
                  int maxBackSteps,
                  int verbosity)
    {
      // Copy input arguments.
      this->m_initialLambda = initialLambda;
      this->m_maxIterations = maxIterations;
      this->m_maxLambda = maxLambda;
      this->m_minLambda = minLambda;
      this->m_minError = minError;
      this->m_minGrad = minimumGradientMagnitude * minimumGradientMagnitude;
      this->m_minDrop = minDrop;
      this->m_strikes = strikes;
      this->m_maxBackSteps = maxBackSteps;
      this->m_verbosity = verbosity;

//     // 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
    OptimizerLM<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>
    OptimizerLM<Functor>&
    OptimizerLM<Functor>::
    operator=(const OptimizerLM<Functor>& source)
    {
      if(&source != this) {
        Optimizer<Functor>::operator=(source);
        m_initialLambda = source.m_initialLambda;
        m_maxBackSteps = source.m_maxBackSteps;
        m_maxIterations = source.m_maxIterations;
        m_maxLambda = source.m_maxLambda;
        m_minDrop = source.m_minDrop;
        m_minError = source.m_minError;
        m_minGrad = source.m_minGrad;
        m_minLambda = source.m_minLambda;
        copyArgumentType(source.m_startPoint, this->m_startPoint);
        m_strikes = source.m_strikes;
        m_verbosity = source.m_verbosity;
      }
      return *this;
    }


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

    template <class Functor>
    std::pair<typename Functor::argument_type, typename Functor::result_type>
    OptimizerLM<Functor>::
    run()
    {
      // Check that we have a valid startPoint.
      if(this->m_startPoint.size() == 0) {
        DLR_THROW3(StateException, "OptimizerLM<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);

      // Initialize variables relating to convergence and convergence
      // failure.
      size_t strikes = 0;
      int backtrackCount = 0;

      // Initialize variables which will take function return values and
      // derivatives.
      result_type errorValue;
      argument_type dEdX(theta.size());
      Array2D<Float64> d2EdX2(theta.size(), theta.size());

      // Initialize intermediate values used by the minimization.
      Array2D<Float64> BMatrix(theta.size(), theta.size());
      Array1D<Float64> deltaX(theta.size());
      argument_type xCond(theta.size());
      Array1D<result_type> errorHistory =
        zeros(m_maxIterations + 1, type_tag<result_type>());
      Float64 lambda = m_initialLambda;

      // Get initial value of error function.
      errorValue = this->m_functor(theta);
      errorHistory[0] = errorValue;
      this->verboseWrite("Error History:\n", errorHistory, 1);

      // Loop until termination.
      size_t iterationIndex = 0;
      for(; iterationIndex < m_maxIterations; ++iterationIndex) {

        // Compute gradient.
        this->m_functor.computeGradientAndHessian(theta, dEdX, d2EdX2);

        // Gradient almost zero?
        if(dotArgumentType(dEdX,dEdX) <= m_minGrad) {
          this->verboseWrite("Tiny gradient, terminating iteration.\n", 1);
          break;
        }

        // Adjust lambda.
        while(lambda <= m_maxLambda) {

          // Precondition the Hessian matrix.
          std::copy(d2EdX2.begin(), d2EdX2.end(), BMatrix.begin());
          for(size_t diagIndex = 0; diagIndex < theta.size(); ++diagIndex) {
            BMatrix(diagIndex, diagIndex) += lambda;
          }

          // Solve B * deltaX = dEdX'
          copyArgumentType(dEdX, deltaX);
          linearSolveInPlace(BMatrix, deltaX);
        
          // We have a new candidate location in the error space.
          for(size_t elementIndex = 0; elementIndex < theta.size();
              ++elementIndex) {
            xCond[elementIndex] = theta[elementIndex] - deltaX[elementIndex];
          }

          // Do we have a decrease in the error function at the candidate
          // location?
          errorValue = m_functor(xCond);
          if(errorValue < errorHistory[iterationIndex]) {
            // Yes. Go on to the next iteration.
            backtrackCount = 0;
            errorHistory[iterationIndex + 1] = errorValue;
            copyArgumentType(xCond, theta);
            lambda /= 10.0;
            if(lambda < m_minLambda) {
              lambda = m_minLambda;
            }
            this->verboseWrite("Lambda = ", lambda, 1);
            this->verboseWrite("Theta = ", theta, 2);
            break;
          } else {
            // Error did not decrease.  Try a bigger lambda.
            ++backtrackCount;
            lambda *= 10.0;
            if(lambda > m_maxLambda) {
              break;
            }
            this->verboseWrite("Lambda = ", lambda, 1);

            // Make sure we haven't exceeded the maxBackSteps
            // termination criterion.
            if(m_maxBackSteps >= 0
               && backtrackCount > m_maxBackSteps) {
              break;
            }

          }
        }

        if(m_verbosity >= 1 ) {
          std::cout << "Error History:\n" << errorHistory << std::endl;
        }
// Note(xxx):  Not sure if we still want this functionality.
//       if(m_iterationFunctorPtr != 0) {
//         m_iterationFunctorPtr(theta);
//       }
    
        // Test termination conditions.
        Float64 drop =
          (errorHistory[iterationIndex] - errorValue)
          / errorHistory[iterationIndex];
        if(drop < m_minDrop) {
          ++strikes;
          if(m_verbosity >= 2) {
            fprintf(stderr, "strikes = %d\n", strikes );
          }
        } else {
          strikes = 0;
        }

        if(lambda >= m_maxLambda 
           || strikes == m_strikes 
           || errorHistory[iterationIndex] <= m_minError 
           || (m_maxBackSteps >= 0 && backtrackCount >= m_maxBackSteps)) {
          if(m_verbosity >= 1) {
            std::cout << "Stopping with lambda = " << lambda
                      << " (" << m_maxLambda << ")\n"
                      << "              strikes = " << strikes
                      << " (" << m_strikes << ")\n"
                      << "              error = " << errorHistory[iterationIndex]
                      << " (" << m_minError << ")\n"
                      << "              backTrackCount = " << backtrackCount
                      << " (" << m_maxBackSteps << ")" << std::endl;
          }
          break;
        }
      }
      return std::make_pair(theta, errorHistory[iterationIndex]);
    }


    template <class Functor>
    inline void
    OptimizerLM<Functor>::
    verboseWrite(const char* message, int verbosity)
    {
      if(verbosity <= this->m_verbosity) {
        std::cout << message << std::flush;
      }
    }


    template <class Functor> template <class Type>
    inline void
    OptimizerLM<Functor>::
    verboseWrite(const char* intro, const Type& subject, int verbosity)
    {
      if(verbosity <= this->m_verbosity) {
        std::cout << intro << subject << std::endl;
      }
    }

  } // namespace optimization

} // namespace dlr

#endif // #ifndef _DLR_OPTIMIZERLM_H_

---