LineSearch.cc

Go to the documentation of this file.

/*
 * $Id: LineSearch.cc 219 2014-11-09 08:59:54Z jdl3 $
 * Copyright (C) 2011, 2013 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifdef HAVE_CONFIG_H
#  include<config.h>
#endif

#include"LineSearch.hpp"

using namespace ipo::detail;

LineSearch::Parameters::Parameters( double const alpha, double const beta,
                    double const epsilon, size_t const maxNoImproveSteps ){
  setAlpha( alpha );
  setBeta( beta );
  setEpsilon( epsilon );
  setMaxNoImproveSteps( maxNoImproveSteps );
  setVectorDigits( 0 );
  setOutputStream( nullptr );
}

LineSearch::Parameters&
LineSearch::getParameters(){ return parameters; }

LineSearch::LineSearch( ipo_function::Function& function, double const alpha,
            double const beta,
            double const epsilon, size_t const maxNoImproveSteps )
  : parameters { alpha, beta, epsilon, maxNoImproveSteps }, function( function ){}

void
LineSearch::Parameters::setAlpha( double alpha ){
  if( alpha <= 0 ) alpha = std::numeric_limits<double>::epsilon();
  if( alpha >= 0.5 ) alpha = 0.5 - std::numeric_limits<double>::epsilon();
  this->alpha = alpha;
}

void
LineSearch::Parameters::setBeta( double beta ){
  if( beta <= 0 ) beta = std::numeric_limits<double>::epsilon();
  if( beta >= 1 ) beta = 1 - std::numeric_limits<double>::epsilon();
  this->beta = beta;
}

void
LineSearch::Parameters::setEpsilon( double epsilon ){
  if( epsilon <= 0 ) epsilon = std::numeric_limits<double>::epsilon();
  this->epsilon = epsilon;
}

double
LineSearch::Parameters::getAlpha() const {
  return alpha;
}

double
LineSearch::Parameters::getBeta() const {
  return beta;
}

double
LineSearch::Parameters::getEpsilon() const {
  return epsilon;
}

double
LineSearch::operator()( gsl::vector& vector, gsl::vector const direction,
            double const functionValue, gsl::vector const gradientValue ){
  // Initialise
  double const& epsilon { parameters.getEpsilon() };
  size_t const SIZE { vector.size() };
  auto const outputStream = parameters.getOutputStream();
  // storage
  gsl::vector testPoint { gsl::vector( SIZE ) };
  // Local function
  auto fDeltat = [&testPoint,this]( gsl::vector const& x, gsl::vector const& delta,
             double const t )->double {
    // compute testPoint
    testPoint.memcpy( delta );
    testPoint.scale( t );
    testPoint.add( x );
    return this->function( testPoint );
  };

  // Get gradient and delta
  auto gradient = [this,&epsilon,&vector,&gradientValue]()->gsl::vector {
    if( gradientValue.get() == nullptr ){
      // Have to compute a gradient value
      ipo_function::detail::ForwardDifferenceGradientEstimate
      gradientEstimate( this->function, epsilon );
      gradientEstimate.setVector( vector );
      return gradientEstimate.gradient();
    } else
      return gradientValue;
  }();

  double d { 0 };
  gsl::blas::ddot( gradient, direction, &d );
  double const delta { parameters.getAlpha() * d };

  // get initial function value
  double const f0 { functionValue > -std::numeric_limits<double>::infinity() ?
      functionValue : function( vector ) };
                                 
  // initialise t to 1
  double t { 1 };

  // Set up to check improvement
  double fBest { f0 };
  double fLast { f0 };
  double tBest { 1 };
  size_t nonImproving { 0 };
  bool feasibleStep { false };

  // loop: check for infinite case
  bool positiveStepSize { false };
  double f { fDeltat( vector, direction, t ) };
  // Possible output
  size_t iteration { 0 };
  size_t const& digits { parameters.getVectorDigits() };
  if( outputStream != nullptr ){
    *outputStream << "Line search:" << std::endl;
    *outputStream << "step "
          << ( digits > 0 ? "x " : "") << "f(x)" << std::endl;
    *outputStream << "  " << std::setw( 4  ) << iteration << " ";
    if( digits > 0 ){
      *outputStream << "( ";
      for( size_t i = 0; i < SIZE; ++i ){
    *outputStream << std::setw( 4 + digits ) << std::setprecision( digits ) << std::fixed
              << testPoint.get( i ) << " ";
      }
      *outputStream << ")";
    }
    size_t const d { std::max( static_cast<size_t>( 6 ), digits ) };
    *outputStream << std::setw( 4 + d ) << std::setprecision( d ) << std::fixed
          << f * functionScale << std::endl;
  }
  while( f > f0 + t * delta or f != f ){ // or condition forces check for NaN
    ++iteration;
    positiveStepSize = true;
    // compute testPoint
    if( f < f0 ) feasibleStep = true;
    // Decide whether improvement has occurred
    if( f < fBest ){
      nonImproving = 0;
      tBest = t;
      fBest = f;
    } else if( f < fLast ){
      nonImproving = 0;
    } else {
      if( feasibleStep ){
    if( ++nonImproving > parameters.getMaxNoImproveSteps() )
      break; // and return
      }
    }
    constexpr double std_eps { std::numeric_limits<double>::epsilon() };
    if( t < std_eps ) break;
    // Record last observed value
    fLast = f;
    // New value of t
    t *= parameters.getBeta();
    // New value of f
    f = fDeltat( vector, direction, t );
    // Possibly put some output
    if( parameters.getOutputStream() != nullptr ){
      *outputStream << "  " << std::setw( 4  ) << iteration << " ";
      if( digits > 0 ){
    *outputStream << "( ";
    for( size_t i { 0 }; i < SIZE; ++i ){
      *outputStream << std::setw( 4 + digits ) << std::setprecision( digits ) << std::fixed
            << testPoint.get( i ) << " ";
    }
    *outputStream << ")";
      }
      size_t const d { std::max( static_cast<size_t>( 6 ), digits ) };
      *outputStream << std::setw( 4 + d ) << std::setprecision( d ) << std::fixed
            << f * functionScale << std::endl;
    }
  }
  if( fBest < f ){
    t = tBest;
    f = fDeltat( vector, direction, t ); // Changes testPoint.
  }

  if( f < f0 and (positiveStepSize or true) ){
    // return testpoint
    vector = testPoint;
    
    // return t
    return t;
  } else {
    return 0;
  }

}

void
LineSearch::Parameters::setMaxNoImproveSteps( size_t maxNoImproveSteps ){
  this->maxNoImproveSteps = maxNoImproveSteps;
}

size_t
LineSearch::Parameters::getMaxNoImproveSteps() const {
  return maxNoImproveSteps;
}