Structural Bioinformatics Library
Template C++ / Python API for developping structural bioinformatics applications.
User Manual

Real_value_function_minimizer

Authors: A. Roth and F. Cazals and T. Dreyfus

This package provides wrappers for libraries minimizing real valued functions at a given point. The two main libraries that are wrapped are :

  • Ipopt : the Interior Point Optimizer library, a compelte toolbox for minimization with or without constraints;
  • LBFGS++ : a header-only C++ library using the LBFGS method method for minimizing.

Introduction

Goals: minima, maxima, and critical points

Given a function $f : A \rightarrow \mathbb{R} $, a global minimum of $f$ is a point $\bar{x} \in A$ such that $\forall x \in A, f(\bar{x}) \leq f(x)$.


Given a function $f : A \rightarrow \mathbb{R} $, a local minimum of $f$ is a point $\bar{x} \in A$ with a neighborhood $V \subset A$ such that $\forall x \in V, f(\bar{x}) \leq f(x)$


When $f$ is convex, any local minimum that is in the interior of $A$ is also a global minimum.However, if $f$ is nonconvex, $f$ may have multiple local minima, with (generically) a single global minimum.

Mathematically, we note that a local minimum is a particular case of a so-called critical point, that is a point where the gradient of the function vanishes. Furthermore, following the Morse lemma, the so-called index of the critical point i.e. the number of negative eigenvalues of the Hessian qualifies its status, from local minimum, to local maximum. Thus, in all generality, methods tracking critical points need to access exact values both for the gradient and the Hessian.

Following classical terminology, a method approximating the Hessian from exact gradients is termed quasi-Newton. For example, the Broyden–Fletcher–Goldfarb–Shanno method is a quasi-Newton method.

In the sequel, when describing methods, by exact partial derivatives, we refer to partial derivatives computed using symbolic differentiation or via automatic differentiation.

Optimization versus biophysics

Optimization and biophysics differ on one fundamental point.

In optimization, one is often satisfied with a (good enough) local minimum (or maximum).

In biophysics, the enumeration of local minima with sufficiently low potential energy is a critical endeavor, as from statistical physics, macroscopic properties emerge from ensemble properties of conformations located in the catchment basins of such local minima. This enumeration task is precisely that tackled by functions from the package Landscape_explorer, where the domain is a conformational space and the function is the potential energy of a target molecule.

We now explain the main algorithms used to minimize a function.

Methods

A large number of methods exist for looking for global or local extrema. Some of them are dedicated to specific functionals (e.g. linear functionals), while some are more general. In this package, our focus is on the latter, and more specifically on iterative methods exploiting 1st and 2nd derivatives .

Two methods are wrapped:

  • IPOPT : implements a primal-dual interior point method.

    • First order derivatives must be exact (as defined in Introduction);
    • Second order derivatives may be exact or may be approximated numerically (whence the quasi-Newton status).

  • L-BFGS : a quasi-Newton optimization algorithm well suited to handle cases with a large number of variables. Indeed, as opposed to BFGS, which stores a fully dimensional inverse Hessian matrix, L-BFGS uses an estimate of this matrix based on a sparse set of vectors.

    • First order derivatives must be exact;
    • Second order derivatives are approximated.

Implementation

This package provides a C++ concept RealValueFunctionMinimizer defining a functor parametrized by the function to optimize : it takes an initial point as input, and returns a pair (bool, point) such that that boolean value is true iff a local minimum has been found, and the point is the local minimum if the local minimum has been found. The requirements of the concept RealValueFunctionMinimizer are :

  • the type Point representing a D-dimensional point;
  • the type FT representing the number type of the coordinates and of the function values;
  • the functor performing the calculations :
std::pair<bool, Point*> operator()(const Point* p)const;

The package provides two models of this concept :

In both cases, the template parameters are :

  • PointType : represents the type of the point,
  • RealValueFunction : functor that computes the function value from an input point
IFT operator()(const InternalPoint* p)const;
  • GradientFunction : functor that computes the gradient q of the input function from an input point p :
void operator()(const InternalPoint* p, InternalPoint* q)const;

Note that both for the function and the gradient functors, the types for numbers and points refer to the internal types to the libraries(Ipopt or LBFGS++), not the types used outside of these libraries. This is important in particular when a particular data structure is required in a program that is not the one used in one of these libraries : in such a case, one has to count the cost of converting its own data structures toward and backward these internal data structures.

Both models share identical tuning parameters :

  • set_tolerance(double) : tunes the tolerance of the algorithm, i.e the epsilon value under which a point is considered a minimum;
  • set_maximum_number_of_iterations(unsigned) : tunes the maximal number of iterations before the algorithm quits with no optimized solution;
  • set_verbose_mode(unsigned) : tunes the verbosity level within the different libraries;
In the case of SBL::CSB::T_Real_value_function_minimizer_Ipopt, the supplemental template parameter MinimizerConstraints is the optimization algorithm from the Ipopt library. It extends the type Ipopt::TNLP as explained in the Ipopt documentation. The default class for this parameter is SBL::CSB::T_Minimizer_Ipopt_parameters_without_constraint , which implements the Ipopt solver with a quasi-Newton algorithm and without any constraint. In order to use Ipopt with a Newton algorithm, or with some constraints, it suffices to provide a model of MinimizerConstraints using the Hessian matrix and / or implementing some constraints, as explained in the documentation of Ipopt.


Examples

Using Ipopt

The following example instantiates the Ipopt algorithm without any constraint for the norm function, and minimizes an input D-dimensional point. It repeats the operation one thousand times for computing an average time of computation.

#include <SBL/CSB/Real_value_function_minimizer_Ipopt.hpp>
#include <SBL/Models/Conformation_traits_vector.hpp>
#include <map>
#include <iostream>
#include <fstream>
#include <chrono>
typedef std::vector<double> Point;
class Energy_evaluator
{
public:
double operator()(double* p)const
{
double e = 0;
for(std::size_t i = 0; i < 2; i++)
e += p[i]*p[i];
return e;
}
};
class Gradient_evaluator
{
public:
void operator()(double* p, double* G)const
{
for(std::size_t i = 0; i < 2; i++)
G[i] = 2*p[i];
}
};
<Point, Energy_evaluator, Gradient_evaluator> Quencher;
int main(int argv, char** argc){
if(argv < 2)
return -1;
Point* p = new Point();
for(std::size_t i = 1; i < 3; i++)
p->push_back(atof(argc[i]));
Quencher quencher;
std::pair<bool, Point*> res = quencher(p);
if(res.first)
std::cout << *res.second << std::endl;
unsigned nb_repeat = 1000;
auto start = std::chrono::steady_clock::now();
for(std::size_t j = 0; j < nb_repeat; j++){
std::pair<bool, Point*> res = quencher(p);
}
auto finish = std::chrono::steady_clock::now();
double elapsed_seconds = std::chrono::duration_cast<std::chrono::duration<double> >(finish - start).count();
std::cout << "average time: " << elapsed_seconds/nb_repeat << "s" << std::endl;
delete p;
return 0;
}

Using LBFGS++

The following example instantiates the LBFGS algorithm, and minimizes an input D-dimensional point. It repeats the operation one thousand times for computing an average time of computation.

#include <SBL/CSB/Real_value_function_minimizer_LBFGS.hpp>
#include <SBL/Models/Conformation_traits_vector.hpp>
#include <map>
#include <iostream>
#include <fstream>
#include <chrono>
typedef std::vector<double> Point;
class Energy_evaluator
{
public:
double operator()(const Point& p)const
{
double e = 0;
for(std::size_t i = 0; i < 2; i++)
e += p[i]*p[i];
return e;
}
};
class Gradient_evaluator
{
public:
void operator()(const Point& p, Point& G)const
{
G.reserve(p.size());
G.resize(p.size());
for(std::size_t i = 0; i < 2; i++)
G[i] = 2*p[i];
}
};
<Point, Energy_evaluator, Gradient_evaluator> Quencher;
int main(int argv, char** argc){
if(argv < 2)
return -1;
Point* p = new Point();
for(std::size_t i = 1; i < 3; i++)
p->push_back(atof(argc[i]));
Quencher quencher;
std::pair<bool, Point*> res = quencher(p);
if(res.first)
std::cout << *res.second << std::endl;
unsigned nb_repeat = 1000;
auto start = std::chrono::steady_clock::now();
for(std::size_t j = 0; j < nb_repeat; j++){
std::pair<bool, Point*> res = quencher(p);
}
auto finish = std::chrono::steady_clock::now();
double elapsed_seconds = std::chrono::duration_cast<std::chrono::duration<double> >(finish - start).count();
std::cout << "average time: " << elapsed_seconds/nb_repeat << "s" << std::endl;
delete p;
return 0;
}

Applications

This package offers base applications for quenching molecular conformations, i.e finding local minima in the vicinity of the given conformations :

  • sbl-bln-quencher.exe : quenches conformations of BLN69
  • sbl-amber-quencher.exe : quenches conformations of proteins using Amber force field
  • sbl-charmm-quencher.exe : quenches conformations of proteins using CHARMM force field