Visualization with xgraphics

LOPTI is mathematical optimization C++ library. There are three major parts. 1st and main part is solver-driver which either wraps external solver or uses LOPTI native solvers and presents to library user consistent and simple optimization API.

2nd part is collection of solvers. Currently 4 solvers are included (or supported): NEWUOA, CONDOR, Nelder-Mead (simplex) and Hook-Jeevs. All these solvers are derivative free and have open source (GPL or BSD like) licences. Some of them (Condor) are unchanged and treated as black-box. Some of them (NEWUOA - fortran source) went through heavy modification to interface with LOPTI. And other (Hook-Jeevs) were completely rewritten into LOPTI native solver.

Though I am primarily interested in derivative-free, unconstrained solvers for expensive objective, LOPTI inherently not limited to these methods and in future there might be other types of solvers and constraints interface. Also I am developing my own solver, preliminary results are very promising, stay tuned.

3rd part is collection of object function and object function modifiers.

Quick Start

#include <lopti/hook-jeevs.h>
using lopti::hook_jeevs_minimizer;                      // minimizer class
using lopti::make_objective;

// Array-type
typedef   lvv::array<double,2>   V;                     // == double    X[2]

// Objective Function
double   rosenbrock_fn (V& X)    {  return  100 * (X[1]-X[0]*X[0])*(X[1]-X[0]*X[0]) + (1-X[0])*(1-X[0]);  };

int main() {

        V       X = {{ -1.2,  1   }};                             // X0 - starting point.
        V       S = {{  0.2,  0.2 }};                             // step0 - initial step

        hook_jeevs_minimizer<V>   m;

        // minimizer config
        m.objective     ( make_objective<V> (&rosenbrock_fn,"rosenbrock") );
        m.x0            (X);
        m.step0         (S);
        m.max_iter      (500);

        // Result
        V       X_min   =  m.argmin();
        double  y_min   =  m.ymin();
        cout  <<  "Minimum: "  <<  y_min  <<  "    at point: "  <<  X_min  <<  endl;
 }

Source of above example is at doc/example.cc. Header-library lvvlib needs to be installed. The doc/Makefile will look for for it in default location /usr/local/include/lvv.

Install

It is possible to use just LOPTI headers. LOPTI library (liblopti.so file) used only by NEWUOA. Running make install installs by default into /usr/local/include and /usr/local/lib For installation external solvers see corresponding section bellow.

You will need

  • Compiler. LOPTI tested only on x86-64 with:

    • GNU GCC with gfortran included (if you want NEWUOA solver), tested with v4.7

    • INTEL ICC, tested with vervion 11

  • lvvlib - you will need array.h for lvv::array class.

  • BOOST C++ libraries, tested with 1.3.37

  • GIT - your distribution probably have package for GIT.

  • XGraphic (cmap.polytechnique.fr/~jouve/xd3d) - simple visualisation utility. This is optional.

  • External solvers (optional), see solvers docs below.

to install:

git clone git://github.com/lvv/lopti.git
cd lopti.git
make install
Warning
A GCC compilation flags

Do no use:

  • -fast-math — newuoa and hook-jeevs can get FP errors

  • -fargument-noalias-anything — newuoa will segfault on exit

try -fno-strict-aliasing -O2 if you encounter stability problem.

Minimizer Class

Data Types

All solvers (and objective) are template with 1st parameter being parameter vector type (named V in source). As of now, the only usable vector type is lvv::array<> which is part of lvvlib. Documentation for array: http://volnitsky.com/project/lvvlib/array.html. It is planned in future to support any random access container having members:

  • value_type

  • size()

  • begin()

  • end()

Common Interface

  • Default constructor (no arguments)

  • X0(V) - starting point

  • S0(V) - starting step size

  • V argmin() - solution optimum point

  • T ymin() - solution minimum value

  • int iter() - number of iterations

Termination Criteria

  • max_iter(int) - defaults to 10000

  • rho_end(T) - defaults to 1/infinity, specific to TR solvers

  • characteristic_size(T) - used only by Nelder-Mead

Solvers

Solver Type Use Derivative Origin Dependency

NEWUOA

Trust Region

No

2002[3] latest from Powell, successor to UOBYQA

none (converted fortran source included)

GSL Nelder-Mead (Simplex)

Direct

No

1965[1]

GNU Scientific Library (package available for most linux distro)

Condor

Trust Region

No

Written by Dr. Ir. Frank Vanden Berghen, based on Powell’s UOBYQA[4]

You need to get source or binary from author

Hook-Jeevs

Direct

No

1962[2], C implementation by Sergey Kiselev

none (converted to native LOPTI implementation)

GNU Scientific Library

GSL (gnu.org/software/gsl) is high quality scientific computing C library. But using it is not simple, at least for simple things. Look at www.gnu.org/software/gsl/manual/html_node/Multimin-Examples.html for Nelder-Mead minimizer. It is 62 line of code without object function, you will need to use cumbersome gsl_vector.

GSL Nelder-Mead (simplex)

#include <lopti/gsl-nelder-mead-wrap.h>
class gsl_nelder_mead_minimizer<V>

LDFLAGS=-lgsl

Old, somewhat slow, can be used with non-smooth functions. Termination criteria is set with member function characteristic_size(double x) Solver iteration is not equal to eval-count, it is about 1.5*eval-count. Rarely but it can happens, algorithm can get stuck, you can unstuck it by restarting form stuck point. Only double vector, index should start from 0.

Trust Region solvers

The Trust-Region (TR) solvers are derived from newton method. TR solvers work by approximating OF with quadratic model function, which is constructed in such way that it passes through evaluation points near current minimum (basis). Near means: at distance less then rho, which is dynamically changed by solver for each iteration. To build an exact model, NP=(N+1)*(N+2) evaluation points are needed, so TR solvers at first NP iterations just evaluate NP points near X0. You can notice this slow start on convergence charts. For high N (>10) this can be especially noticeable. Instead of exact model it is possible with NEWUOA to use approximate model with less then (N+1)*(N+2) points.

Uses member rho_begin(T), instead of S0(V), defaults to 0.2. Because rho is scalar (applied equally to all directions) it is recommended to normalize object junction with [rescale]. Termination criteria is rho_end(T), defaults to 1/infinity.

NEWUOA

#include <lopti/newuoa.h>
class newuoa_minimizer<V,[NP=2*N+1]>

LDFLAGS=-llopti -lgfortran

One of the best solvers in this collection. Optional second template parameters NP - is Number-of-Points used for building model, defaults to 2*N+1 (approximate model). max_iter should be bigger then NP.

Caution If second template paramter specifies exact model, with N > 14, NEWUOA crashes (index out of range).

When this solver is used LDFLAGS=-llopti needed. Any lvv::array<> vector.

CONDOR

#include <lopti/condor.h>
class condor_minimizer<V,[NP=2*N+1>

LDFLAGS=-lcondor CXXFLAGS=-I/path/to/dir/with/condorlib-dir

Fast for low dimensions. Constraints capability not used. Condor internal rescaler not used (use LOPTI rescaler). Can use only exact model, so it use might be problematic for big N. LOPTI will look for <condor/Vector.h> from condorlib subdirectory in Condor source directory.

Caution Sometimes Condor can get stuck at high precision, it is not clear if it is rounding error or algorithm limitation. Also it sometimes aborts at FP illegal ops in the middle of optimization run.

Only double vector, index should start from 0.

Hook-Jeevs

Native LOPTI solver, converted from Oleg Kiselev, 1998 implementation (http://netlib.org/c++/linalg.tgz) . Old and usually slow solver. But sometimes can outperform all other at not trivial high dimension tasks.

Caution Might segfault. Try lowering optimization from -O3 to -O2 to fix.

Code still needs some style cleanup, and thread safety.

Object Function

Object function are defined in class derived from class objective0. Object function also can be wrapped in multiple wrappers, such as rescaler, logger, tracer, adapter, etc. Planned wrappers are also noise-adder, constrain-adder-by-barrier-function,

make_objective() - create objective from any callable object

Is optional convenience class which can convert plain function to LOPTI objective class.

make_objective(std::function<T (*)(V& X)> callable_object)

It can convert any callable object to objective class. Callable object can be plain function, class member function, or functor. A callable object passed (and possibly casted) to make_objective constructor. Because make_objective expects a std::function as its argument, it is converted to std::function internally.

double   plain_function(lvv::array<double,2> X);
...
mzr.objective ( make_objective<V> (&plain_function);

trace() - progress tracer

trace(objective o,)

Prints on screen eval-number, seconds of last eval, X and f-value

xg_log() - xgraphic logger

xg_log(objective o, minimizer<V> m)

Logs eval-number, f-value and X on every objective evaluation. Program xgraphic is open source visualization utility which can directly input log file. Created log name is composed from objective name, number of dimentions and minimizer name. Logs are created under ./log/

rescale() - normalizer

rescale<V>(objective o, V rescale_vector)

Most solvers work much better if function first derivative (and second?) have about the same order in all directions. The rescale converts object function into such normalized function. When solver request to evaluate objective - f(X), rescale gives to rescale-wrapped objective modified X. Modification is simple multiplication by rescaling factor: X*rescale_vector. When you use rescale, don’t forget also rescale X0 by deviding on rescale vector and answer argmin() by multiplying by rescale vector.

Rescale
typedef  array<double,2>   V;
V   R  = {{ 1, 0.0100 }};               // rescale vector
V   X0 = {{ -1.2, 1 }};
V   X_opt;                              // answer
condor_minimizer<V0>    mzr;
        mzr     .objective              (rescale<V0>( trace<V0>(rosenbrock<V0>()), R));
        mzr     .x0             (X/=R);
        mzr     .rho_begin      (1);
        mzr     .rho_end        (0.01);
X_opt = mzr     .argmin();
X_opt *= R;                             // un-rescale
cout << "optimum at: " << X_opt << endl;

How to submit patch

You can just email patch to leonid@volnitsky.com. All feedback is much appreciated.

If you are on Github it is even easier. See github patch-submit HOWTOs: ( 1, 2, 3 ) .

There is no hard set style rules.

TODO

  • fix all type-puning

  • add box and polyhedron constraints

  • remove using from headers

  • unify and simplify termination criteria.

  • add mkl TR solver solver

  • add zib nlqe1 solver

  • add zib newton_test functions

  • add gsl bfgs solver

  • add assert / throw on use uninitialized data in minimizer class

Microformat

LOPTI
C++ library
Mathematical Optimization Labrary
http://volnitsky.com/project/lopti
http://volnitsky.com/project/lopti/xgr001.gif

References

  • [1] J.A. Nelder and R. Mead, A simplex method for function minimization, Computer Journal vol. 7 (1965), 308-315.

  • [2] R. Hooke and T. A. Jeeves, Direct Search Solution of Numerical and Statistical Problems, Journal of the ACM, Vol. 8, April 1961, pp. 212-229

  • [3] M. J. D. Powell. The NEWUOA software for unconstrained optimization with derivatives. DAMTP Report 2004/NA05, University of Cambridge, 2004.

  • [4] Frank Vanden Berghen, Hugues Bersini, CONDOR, a new parallel, constrained extension of Powell’s UOBYQA algorithm: Experimental results and comparison with the DFO algorithm, Journal of Computational and Applied Mathematics, Elsevier, Volume 181, Issue 1, September 2005, Pages 157-175

  • [5] Intel Math Kernel Library Reference Manual. Optimization Solver Routines: www.intel.com/software/products/mkl/docs/WebHelp/osr/osr_Intro.html

  • [6] Stefen Boyd, Stanford Engineering Open Courseware, Convex Optimization - see.stanford.edu/see/lecturelist.aspx?coll=2db7ced4-39d1-4fdb-90e8-364129597c87

    • Backtraking Line Search, Armijo-Goldstein condition (ee364a lecture 15, 18min)

    • Gradient Descent method (ee364a lecture 15, 24 min)

    • Steepest Descent method (ee364a lecture 15, 32.1 min)

  • [7] Anita H.M, Numerical Methods for Scientists and Engineers MacGraw-Hill, 1991, ISBN 0074600133

    • Newton, Qausi-Newton, BFG, BFGS (p366)

    • Direction Set Methods (p370)

    • Conjugate Gradient method (p376)

  • [8] Bjoerck A., Dahlquist G.Numerical mathematics and scientific computation (web draft, 1999) Vols.2,3.

    • Steepest Descent (p403, §11.2.2)

    • Newton (p401, §11.2.3)

    • Quasi-Newton Methods (p406, §11.2.3)

  • [9] S.D.Conte, Carl de Boor Elementary Numerical Analysis - An Algorithmic Approach MacGraw-Hill 1980 3rd-Ed

    • Steepest Descent (p211 5.1)

  • [10] W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes - The Art of Scientific Computing 3nd ed, Cambridge University Press 2007(1988), ISBN-13 978-0-511-33555-6, ISBN-13 978-0-521-88068-8

    • Simplex Method (p502 §10.5)

    • Direction Set (Powell’s) method (p509 §10.7)

    • Quasi-Newton (BFGS) (p521 §10.9)

  • [11] Quarteroni A., Sacco R., Saleri F. Numerical mathematics, 2nd ed., Springer 2007

    • The Hooke and Jeeves Method (p300 §7.2.1)

    • Descent Methods (good overview, Newton, quasi-Newton, Gradient, Conjugate Gradient method) (p306 §7.2.2)

    • Newton (p313 §7.2.6)

    • quasi-Newton (p313 §7.2.7)

  • [12] Schwartz R. Biological modeling and simulation. A survey of practical models, algorithms, and numerical methods MIT-2008, ISBN 0262195844

    • The Levenberg–Marquardt Method (p90 §5.6.2)

  • [13] Zarowsky C.J. An introduction to numerical analysis for electrical and computer engineers, Wiley 2004

    • Backtraking Line Search (p353)

    • Newton’s algorithm with the backtracking (p353 §8.4)

  • [14] Wolfram Alpha: Examples: Optimization http://www.wolframalpha.com/examples/Optimization.html ; http://www.wolframalpha.com/input/?i=local+extrema+sin+x^2

  • [15] 'Mixing C++ an Fortran http://solarianprogrammer.com/2012/05/11/mixed-language-programming-cpp-11-fortran-2008/