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.

Caution Though LOPTI is quite usable and everything works for me, this is still work in progress, documentation is incomplete, API have not stabilized, it requires GCC-4.4 or ICC-11 and some C++0x features are used.

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

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:

Common Interface

Termination Criteria

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 discontinues 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 segfaulted. Lowerer 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(boost::function<T (*)(V& X)> callable_object)

It can convert any callable object (www.boost.org/doc/libs/1_37_0/doc/html/function/tutorial.html) 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 boost::function as its argument, it is converted to boost::function internally.

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

See also boost function_tutorial - www.boost.org/doc/libs/1_37_0/doc/html/function/tutorial.html.

Objective should have one argument - optimization vector. If your original function have more arguments use boost::bind() or std::bind1st(), bind2end()

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, 4 ) .

There is no hard set style rules. I am trying to keep tabstops=8 and expandtab (tab character replaced with 8 spaces).

TODO

Microformat

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

References