A base class for handling optimisation of scalar functions. More...

#include <GslOptimizer.h>

Inheritance diagram for QUESO::GslOptimizer:

Collaboration diagram for QUESO::GslOptimizer:

Public Types
enum	SolverType { FLETCHER_REEVES_CG, POLAK_RIBIERE_CG, BFGS, BFGS2, STEEPEST_DESCENT, NELDER_MEAD, NELDER_MEAD2, NELDER_MEAD2_RAND }

Public Member Functions
	GslOptimizer (const BaseScalarFunction< GslVector, GslMatrix > &objectiveFunction)
	Constructs an object that will maximize a scalar function. More...

virtual	~GslOptimizer ()
	Destructor. More...

virtual void	minimize (OptimizerMonitor *monitor=NULL)
	Minimize the objective function, starting at `m_initialPoint`. More...

const BaseScalarFunction < GslVector, GslMatrix > &	objectiveFunction () const
	Returns the objective function. More...

void	setInitialPoint (const GslVector &intialPoint)
	Set the point at which the optimization starts. More...

const GslVector &	minimizer () const
	Return the point that minimizes the objective function. More...

void	set_solver_type (SolverType solver)

void	set_solver_type (std::string &solver)

SolverType	string_to_enum (std::string &solver)

void	set_step_size (const GslVector &step_size)
	Sets step size used in gradient-free solvers. More...

void	set_step_size (double step_size)
	Sets step size used in gradient-based solvers. More...

void	set_line_tol (double tol)
	Set GSL line minimization tolerance. More...

Public Member Functions inherited from QUESO::BaseOptimizer
	BaseOptimizer ()
	Default constructor. More...

virtual	~BaseOptimizer ()
	Destructor. More...

unsigned int	getMaxIterations () const
	Returns the maximum number of iterations the optimizer will do. More...

double	getTolerance () const
	Returns the tolerance used to test for an extremum in the optimizer. More...

double	getFiniteDifferenceStepSize () const
	Returns the step size used in the finite difference formula. More...

void	setMaxIterations (unsigned int maxIterations)
	Sets the maximum number of iterations to be used by the optimizer. More...

void	setTolerance (double tolerance)
	Sets the tolerance the optimizer will use to test for an extremum. More...

void	setFiniteDifferenceStepSize (double h)
	Sets the step to use in the finite difference derivative. More...

Private Member Functions
bool	solver_needs_gradient (SolverType solver)
	Helper function. More...

void	minimize_with_gradient (unsigned int dim, OptimizerMonitor *monitor)

void	minimize_no_gradient (unsigned int dim, OptimizerMonitor *monitor)

Private Attributes
const BaseScalarFunction < GslVector, GslMatrix > &	m_objectiveFunction

GslVector *	m_initialPoint

GslVector *	m_minimizer

SolverType	m_solver_type

GslVector	m_fstep_size
	For use in gradient-free algorithms. More...

double	m_fdfstep_size
	For use in gradient-based algorithms. More...

double	m_line_tol
	Line minimization tolerance in gradient-based algorithms. More...

Additional Inherited Members
Protected Attributes inherited from QUESO::BaseOptimizer
unsigned int	m_maxIterations

double	m_tolerance

double	m_finiteDifferenceStepSize

Detailed Description

A base class for handling optimisation of scalar functions.

WRITE DOCS HERE

Definition at line 50 of file GslOptimizer.h.

Member Enumeration Documentation

enum QUESO::GslOptimizer::SolverType

Enumerator
FLETCHER_REEVES_CG
POLAK_RIBIERE_CG
BFGS
BFGS2
STEEPEST_DESCENT
NELDER_MEAD
NELDER_MEAD2
NELDER_MEAD2_RAND

Definition at line 71 of file GslOptimizer.h.

                   { FLETCHER_REEVES_CG,
                     POLAK_RIBIERE_CG,
                     BFGS,
                     BFGS2,
                     STEEPEST_DESCENT,
                     NELDER_MEAD,
                     NELDER_MEAD2,
                     NELDER_MEAD2_RAND };

Constructor & Destructor Documentation

QUESO::GslOptimizer::GslOptimizer ( const BaseScalarFunction< GslVector, GslMatrix > & objectiveFunction )

Constructs an object that will maximize a scalar function.

The function objectiveFunction is the function that will be maximized.

Definition at line 159 of file GslOptimizer.C.

References QUESO::GslVector::cwSet(), m_fstep_size, and m_minimizer.

   : BaseOptimizer(),
     m_objectiveFunction(objectiveFunction),
     m_initialPoint(new GslVector(objectiveFunction.domainSet().
           vectorSpace().zeroVector())),
     m_minimizer(new GslVector(this->m_objectiveFunction.domainSet().
         vectorSpace().zeroVector())),
     m_solver_type(BFGS2),
     m_fstep_size(this->m_objectiveFunction.domainSet().vectorSpace().zeroVector()),
     m_fdfstep_size(1.0),
     m_line_tol(0.1)
 {
   // We initialize the minimizer to GSL_NAN just in case the optimization fails
   m_minimizer->cwSet(GSL_NAN);
 
   // Set to documented default value.
   m_fstep_size.cwSet(0.1);
 }

QUESO::GslOptimizer::~GslOptimizer ( )

virtual

Destructor.

Definition at line 179 of file GslOptimizer.C.

References m_initialPoint.

 {
   delete this->m_initialPoint;
 }

Member Function Documentation

void QUESO::GslOptimizer::minimize ( OptimizerMonitor * monitor = NULL )

virtual

Minimize the objective function, starting at m_initialPoint.

m_initialPoint is handled in the derived class

Implements QUESO::BaseOptimizer.

Definition at line 184 of file GslOptimizer.C.

References dim, m_initialPoint, m_objectiveFunction, m_solver_type, minimize_no_gradient(), minimize_with_gradient(), queso_error, and solver_needs_gradient().

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::solveWithBayesMetropolisHastings().

                                                      {
 
   // First check that initial guess is reasonable
   if (!this->m_objectiveFunction.domainSet().contains(*(this->m_initialPoint)))
     {
       if( m_objectiveFunction.domainSet().env().fullRank() == 0 )
         {
           std::cerr << "Minimization was given initial point outside of domain"
                     << std::endl;
         }
       queso_error();
     }
 
   unsigned int dim = this->m_objectiveFunction.domainSet().vectorSpace().
     zeroVector().sizeLocal();
 
   if( this->solver_needs_gradient(m_solver_type) )
     {
       this->minimize_with_gradient( dim, monitor );
     }
   else
     {
       this->minimize_no_gradient( dim, monitor );
     }
 
   return;
 }

void QUESO::GslOptimizer::minimize_no_gradient	(	unsigned int	dim,
		OptimizerMonitor *	monitor
	)

private

Definition at line 365 of file GslOptimizer.C.

References QUESO::OptimizerMonitor::append(), BFGS, BFGS2, QUESO::c_evaluate(), dim, FLETCHER_REEVES_CG, QUESO::BaseOptimizer::getMaxIterations(), QUESO::BaseOptimizer::getTolerance(), m_fstep_size, m_initialPoint, m_objectiveFunction, m_solver_type, NELDER_MEAD, NELDER_MEAD2, NELDER_MEAD2_RAND, POLAK_RIBIERE_CG, queso_error, and STEEPEST_DESCENT.

Referenced by minimize().

   {
     // Set initial point
     gsl_vector* x = gsl_vector_alloc(dim);
     for (unsigned int i = 0; i < dim; i++) {
       gsl_vector_set(x, i, (*m_initialPoint)[i]);
     }
 
     // Tell GSL which solver we're using
     const gsl_multimin_fminimizer_type* type = NULL;
 
     switch(m_solver_type)
       {
       case(NELDER_MEAD):
         type = gsl_multimin_fminimizer_nmsimplex;
         break;
       case(NELDER_MEAD2):
         type = gsl_multimin_fminimizer_nmsimplex2;
         break;
       case(NELDER_MEAD2_RAND):
         type = gsl_multimin_fminimizer_nmsimplex2rand;
         break;
       case(FLETCHER_REEVES_CG):
       case(POLAK_RIBIERE_CG):
       case(BFGS):
       case(BFGS2):
       case(STEEPEST_DESCENT):
       default:
         // Wat?!
         queso_error();
       }
 
     // Init solver
     gsl_multimin_fminimizer* solver =
       gsl_multimin_fminimizer_alloc(type, dim);
 
     // Point GSL at the right functions
     gsl_multimin_function minusLogPosterior;
     minusLogPosterior.n = dim;
     minusLogPosterior.f = &c_evaluate;
     minusLogPosterior.params = (void *)(this);
 
     // Needed for these gradient free algorithms.
     gsl_vector* step_size = gsl_vector_alloc(dim);
 
     for(unsigned int i = 0; i < dim; i++) {
       gsl_vector_set(step_size, i, m_fstep_size[i]);
     }
 
     gsl_multimin_fminimizer_set(solver, &minusLogPosterior, x, step_size);
 
     int status;
     size_t iter = 0;
     double size = 0.0;
 
     do
       {
         iter++;
         status = gsl_multimin_fminimizer_iterate(solver);
 
         if (status) {
           if( m_objectiveFunction.domainSet().env().fullRank() == 0 )
             {
               std::cerr << "Error while GSL does optimisation. "
                         << "See below for GSL error type." << std::endl;
               std::cerr << "Gsl error: " << gsl_strerror(status) << std::endl;
             }
           break;
         }
 
         size = gsl_multimin_fminimizer_size(solver);
 
         status = gsl_multimin_test_size (size, this->getTolerance());
 
         if(monitor)
         {
           gsl_vector* x = gsl_multimin_fminimizer_x(solver);
           std::vector<double> x_min(dim);
           for( unsigned int i = 0; i < dim; i++)
             x_min[i] = gsl_vector_get(x,i);
 
           double f = gsl_multimin_fminimizer_minimum(solver);
 
           monitor->append( x_min, f, size );
         }
 
       }
 
     while ((status == GSL_CONTINUE) && (iter < this->getMaxIterations()));
 
     for (unsigned int i = 0; i < dim; i++) {
       (*m_minimizer)[i] = gsl_vector_get(solver->x, i);
     }
 
     // We're being good human beings and cleaning up the memory we allocated
     gsl_vector_free(step_size);
     gsl_multimin_fminimizer_free(solver);
     gsl_vector_free(x);
 
     return;
   }

void QUESO::GslOptimizer::minimize_with_gradient	(	unsigned int	dim,
		OptimizerMonitor *	monitor
	)

private

Definition at line 268 of file GslOptimizer.C.

References QUESO::OptimizerMonitor::append(), BFGS, BFGS2, QUESO::c_evaluate(), QUESO::c_evaluate_derivative(), QUESO::c_evaluate_with_derivative(), dim, FLETCHER_REEVES_CG, QUESO::BaseOptimizer::getMaxIterations(), m_fdfstep_size, m_initialPoint, m_line_tol, m_objectiveFunction, m_solver_type, NELDER_MEAD, NELDER_MEAD2, NELDER_MEAD2_RAND, POLAK_RIBIERE_CG, queso_error, and STEEPEST_DESCENT.

Referenced by minimize().

   {
     // Set initial point
     gsl_vector * x = gsl_vector_alloc(dim);
     for (unsigned int i = 0; i < dim; i++) {
       gsl_vector_set(x, i, (*m_initialPoint)[i]);
     }
 
     // Tell GSL which solver we're using
     const gsl_multimin_fdfminimizer_type* type = NULL;
 
     switch(m_solver_type)
       {
       case(FLETCHER_REEVES_CG):
         type = gsl_multimin_fdfminimizer_conjugate_fr;
         break;
       case(POLAK_RIBIERE_CG):
         type = gsl_multimin_fdfminimizer_conjugate_pr;
         break;
       case(BFGS):
         type = gsl_multimin_fdfminimizer_vector_bfgs;
         break;
       case(BFGS2):
         type = gsl_multimin_fdfminimizer_vector_bfgs2;
         break;
       case(STEEPEST_DESCENT):
         type = gsl_multimin_fdfminimizer_steepest_descent;
         break;
       case(NELDER_MEAD):
       case(NELDER_MEAD2):
       case(NELDER_MEAD2_RAND):
       default:
         // Wat?!
         queso_error();
       }
 
     // Init solver
     gsl_multimin_fdfminimizer * solver =
       gsl_multimin_fdfminimizer_alloc(type, dim);
 
     // Point GSL to the right functions
     gsl_multimin_function_fdf minusLogPosterior;
     minusLogPosterior.n = dim;
     minusLogPosterior.f = &c_evaluate;
     minusLogPosterior.df = &c_evaluate_derivative;
     minusLogPosterior.fdf = &c_evaluate_with_derivative;
     minusLogPosterior.params = (void *)(this);
 
     gsl_multimin_fdfminimizer_set(solver, &minusLogPosterior, x, m_fdfstep_size, m_line_tol);
 
     int status;
     size_t iter = 0;
 
     do {
       iter++;
       status = gsl_multimin_fdfminimizer_iterate(solver);
 
       if (status) {
         if( m_objectiveFunction.domainSet().env().fullRank() == 0 )
           {
             std::cerr << "Error while GSL does optimisation. "
                       << "See below for GSL error type." << std::endl;
             std::cerr << "Gsl error: " << gsl_strerror(status) << std::endl;
           }
         break;
       }
 
       status = gsl_multimin_test_gradient(solver->gradient, this->getTolerance());
 
       if(monitor)
         {
           gsl_vector* x = gsl_multimin_fdfminimizer_x(solver);
           std::vector<double> x_min(dim);
           for( unsigned int i = 0; i < dim; i++)
             x_min[i] = gsl_vector_get(x,i);
 
           double f = gsl_multimin_fdfminimizer_minimum(solver);
 
           gsl_vector* grad = gsl_multimin_fdfminimizer_gradient(solver);
           double grad_norm = gsl_blas_dnrm2(grad);
 
           monitor->append( x_min, f, grad_norm );
         }
 
     } while ((status == GSL_CONTINUE) && (iter < this->getMaxIterations()));
 
     for (unsigned int i = 0; i < dim; i++) {
       (*m_minimizer)[i] = gsl_vector_get(solver->x, i);
     }
 
     // We're being good human beings and cleaning up the memory we allocated
     gsl_multimin_fdfminimizer_free(solver);
     gsl_vector_free(x);
 
     return;
   }

const GslVector & QUESO::GslOptimizer::minimizer ( ) const

Return the point that minimizes the objective function.

This state will be filled with GSL_NAN if, for some reason, the optimization failed

Definition at line 227 of file GslOptimizer.C.

References m_minimizer.

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::solveWithBayesMetropolisHastings().

 {
   return *(this->m_minimizer);
 }

const BaseScalarFunction< GslVector, GslMatrix > & QUESO::GslOptimizer::objectiveFunction ( ) const

Returns the objective function.

Definition at line 213 of file GslOptimizer.C.

References m_objectiveFunction.

Referenced by QUESO::c_evaluate(), and QUESO::c_evaluate_derivative().

 {
   return this->m_objectiveFunction;
 }

void QUESO::GslOptimizer::set_line_tol ( double tol )

Set GSL line minimization tolerance.

Applicable only to gradient-based solvers. Default is 0.1, as recommended by GSL documentation. See GSL documentation for more details.

void QUESO::GslOptimizer::set_solver_type ( SolverType solver )

Definition at line 232 of file GslOptimizer.C.

References m_solver_type.

Referenced by set_solver_type().

   {
     m_solver_type = solver;
   }

void QUESO::GslOptimizer::set_solver_type ( std::string & solver )

Definition at line 517 of file GslOptimizer.C.

References set_solver_type(), and string_to_enum().

   {
     this->set_solver_type( this->string_to_enum(solver) );
   }

void QUESO::GslOptimizer::set_step_size ( const GslVector & step_size )

Sets step size used in gradient-free solvers.

By default, the step size used will be a vector of 0.1. Use this method to reset the step_size to the desired values.

Definition at line 467 of file GslOptimizer.C.

References m_fstep_size.

   {
     m_fstep_size = step_size;
   }

void QUESO::GslOptimizer::set_step_size ( double step_size )

Sets step size used in gradient-based solvers.

GSL doesn't document this parameter well, but it seems to be related to the line search, so we default to 1.0 for full step.

Definition at line 472 of file GslOptimizer.C.

References m_fdfstep_size.

   {
     m_fdfstep_size = step_size;
   }

void QUESO::GslOptimizer::setInitialPoint ( const GslVector & intialPoint )

Set the point at which the optimization starts.

Definition at line 219 of file GslOptimizer.C.

References m_initialPoint, and QUESO::GslVector::sizeLocal().

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::solveWithBayesMetropolisHastings().

 {
   for (unsigned int i = 0; i < initialPoint.sizeLocal(); i++) {
     (*(this->m_initialPoint))[i] = initialPoint[i];
   }
 }

bool QUESO::GslOptimizer::solver_needs_gradient ( SolverType solver )

private

Helper function.

Definition at line 237 of file GslOptimizer.C.

References BFGS, BFGS2, FLETCHER_REEVES_CG, NELDER_MEAD, NELDER_MEAD2, NELDER_MEAD2_RAND, POLAK_RIBIERE_CG, queso_error, and STEEPEST_DESCENT.

Referenced by minimize().

   {
     bool gradient_needed = false;
 
     switch(solver)
       {
       case(FLETCHER_REEVES_CG):
       case(POLAK_RIBIERE_CG):
       case(BFGS):
       case(BFGS2):
       case(STEEPEST_DESCENT):
         {
           gradient_needed = true;
           break;
         }
       case(NELDER_MEAD):
       case(NELDER_MEAD2):
       case(NELDER_MEAD2_RAND):
         {
           break;
         }
       default:
         {
           // Wat?!
           queso_error();
         }
       } // switch(solver)
 
     return gradient_needed;
   }

GslOptimizer::SolverType QUESO::GslOptimizer::string_to_enum ( std::string & solver )

Definition at line 477 of file GslOptimizer.C.

References BFGS, BFGS2, FLETCHER_REEVES_CG, m_objectiveFunction, NELDER_MEAD, NELDER_MEAD2, NELDER_MEAD2_RAND, POLAK_RIBIERE_CG, queso_error, and STEEPEST_DESCENT.

Referenced by set_solver_type().

   {
     SolverType solver_type;
 
     if( solver == std::string("fletcher_reeves_cg") )
       solver_type = FLETCHER_REEVES_CG;
     else if( solver == std::string("polak_ribiere_cg") )
       solver_type = POLAK_RIBIERE_CG;
     else if( solver == std::string("bfgs") )
       solver_type = BFGS;
     else if( solver == std::string("bfgs2") )
       solver_type = BFGS2;
     else if( solver == std::string("steepest_decent") )
       solver_type = STEEPEST_DESCENT;
     else if( solver == std::string("nelder_mead") )
       solver_type = NELDER_MEAD;
     else if( solver == std::string("nelder_mead2") )
       solver_type = NELDER_MEAD2;
     else if( solver == std::string("nelder_mead2_rand") )
       solver_type = NELDER_MEAD2_RAND;
     else
       {
         if( m_objectiveFunction.domainSet().env().fullRank() == 0 )
           {
             std::cerr << "Error: Invalid GslOptimizer solver name: " << solver << std::endl
                       << "       Valids choices are: fletcher_reeves_cg" << std::endl
                       << "                           polak_ribiere_cg" << std::endl
                       << "                           bfgs" << std::endl
                       << "                           bfgs2" << std::endl
                       << "                           steepest_decent" << std::endl
                       << "                           nelder_mead" << std::endl
                       << "                           nelder_mead2" << std::endl
                       << "                           nelder_mead2_rand" << std::endl;
           }
         queso_error();
       }
 
     return solver_type;
   }