N-dimensional test functions for optimization, in Python

logsumexp.py

""" logsumexp ~ Su Boyd Candes, Diff eq modeling Nesterov accelerated gradient, 2014, p. 7

    Logsumexp( seed: int / RandomState / env SEED  / 0 )
"""

from __future__ import division
import os
import numpy as np

__version__ = "2015-01-31 jan  denis-bz-py at t-online.de"


# logsumexp = Logsumexp() below: logsumexp(x), logsumexp.gradient(x)

#...............................................................................
class Logsumexp( object ):
    __doc__ = globals()["__doc__"]

    def __init__( s, seed=None, m=1, n=1, rho=20 ):
        s.rho = rho
        s.__name__ = "logsumexp"
        s.reset( seed )
        s._initrandomAb( n, m )

    #...........................................................................
    def __call__( s, x ):
        ax_b = s._ax_b( x )
        exp = np.exp( ax_b / s.rho )
        return s.rho * np.log( np.sum( exp ))

    def gradient( s, x ):
        ax_b = s._ax_b( x )
        exp = np.exp( ax_b / s.rho )
        return s.A.T.dot( exp ) / np.sum( exp )

    def _ax_b( s, x ):
        x = np.asarray_chkfinite(x)
        if len(x) != s.A.shape[1]:
            s._initrandomAb( len(x) )
        return s.A.dot(x) - s.b

    def reset( s, seed ):
        if seed is None:
            seed = int( os.getenv( "SEED", 0 ))
        if not isinstance(seed, np.random.RandomState):
            seed = np.random.RandomState( seed=seed )
        s.randomstate = seed
        s._initrandomAb( 1, 1 )

    def _initrandomAb( s, n, m=0 ):
        """ init random A b on first call each size """
            # bug: calls 3d 4d 3d 4d ... todo: dict n -> A b
        if m == 0:
            m = 4 * n  # 200 x 50
        s.A = s.randomstate.normal( size=(m, n) )
        s.b = s.randomstate.normal( scale=np.sqrt(2), size=m )


logsumexp = Logsumexp()  # logsumexp(x), logsumexp.gradient(x)


#...............................................................................
if __name__ == "__main__":
    import sys

    def avmax( x ):
        absx = np.fabs(x)
        return "|av| %.3g  max %.3g" % (absx.mean(), absx.max())

#...............................................................................
    nn = [50]
    nseed = 3

        # to change these params in sh or ipython, run this.py  a=1  b=None  c=[3] ...
    for arg in sys.argv[1:]:
        exec( arg )

    np.set_printoptions( threshold=100, edgeitems=3, linewidth=120,
        formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g

#...............................................................................
    for n in nn:
      for seed in range(nseed):
        print "\nlogsumexp( n %d  seed %s ) --" % (n, seed)
        logsumexp = Logsumexp( seed=seed )
        x = np.zeros(n)
        f = logsumexp( x )
        g = logsumexp.gradient( x )
        print "f(0): %.3g" % f
        print "grad: %s  %s" % (avmax(g), g)

Nd-testfuncs-python.md

N-dimensional test functions for optimization, in Python.

These are the n-dim Matlab functions by A. Hedar (2005), translated to Python-numpy.
http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO.htm
ackley dp griew levy mich perm powell power rast rosen schw sphere sum2 trid zakh .m + ellipse nesterov powellsincos

http://imgur.com/Iu9H0jk is a plot of 3 no-derivative algorithms, NLOpt LN_BOBYQA LN_PRAXIS LN_SBPLX, on these functions, in 8d, from 10 random start points. This plot shows, yet again, that

• comparing methods -- which is "best" -- is really difficult
• random effects, here the 10 random start points x0, can easily mask method differences.
  (A rule of thumb: Wendel's theorem suggests that one should take at least 2*dim random x0.)

Ideas wanted

• better ways to plot methods / treatments across many test cases
• how to find test functions similar to a given real one -- for optimization, a gallery of function landscapes ?
• how to plot convergence paths in 10d.

Comments welcome.

cheers
-- denis 8 Sept 2014

Ndtestfuncs.py

#!/usr/bin/env python
""" some n-dimensional test functions for optimization in Python.
    Example:

    import numpy as np
    from ... import ndtestfuncs  # this file, ndtestfuncs.py

    funcnames = "ackley ... zakharov"
    dim = 8
    x0 = e.g. np.zeros(dim)
    for func in ndtestfuncs.getfuncs( funcnames ):
        fmin, xmin = myoptimizer( func, x0 ... )
        # calls func( x ) with x a 1d numpy array or array-like of any size >= 2

These are the n-dim Matlab functions by A. Hedar (2005), translated to Python-numpy.
http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO.htm
    ackley.m dp.m griew.m levy.m mich.m perm.m powell.m power.m
    rast.m rosen.m schw.m sphere.m sum2.m trid.m zakh.m
    + ellipse nesterov powellsincos randomquad logsumexp

--------------------------------------------------------------------------------
    All functions appearing in this work are fictitious;
    any resemblance to real-world functions, living or dead, is purely coincidental.
--------------------------------------------------------------------------------

Notes
-----

Each `func( x )` works for `x` of any size >= 2.
Each starts off
    x = np.asarray_chkfinite(x)  # ValueError if any NaN or Inf

The values of most functions increase as O(dim) or O(dim^2),
so the convergence curves for dim 5 10 20 ... are not comparable,
and `ftol_abs` depends on func() scaling.
Better would be to scale function values to min 1, max 100 in all dimensions.
Similarly, `xtol_abs` depends on `x` scaling;
`x` should be scaled to -1 .. 1 in all dimensions.
(`ftol_rel` and `xtol_rel` can misbehave near 0; see `isclose`, abs or rel.)

Results from any optimizer depend of course on `ftol_abs xtol_abs maxeval ...`
plus hidden or derived parameters, e.g. BOBYQA rho.
Methods like Nelder-Mead that track sets of points, starting with `x0 + initstep I`,
are sensitive to `initstep`. And what are `ftol` and `xtol` for sets ?

Some functions have many local minima or saddle points (more in higher dimensions ?),
making the final fmin very sensitive to the starting x0.
Also, some have a local minimum at [0 0 ...] so starting there is boring.
*Always* look at a few points near a purported xmin.

Fun constrained problems: min f(x) over the surface of f's bounding box.


See also
--------

http://en.wikipedia.org/wiki/Test_functions_for_optimization  -- 2d plots
http://www.scipy.org/NumPy_for_Matlab_Users

nlopt/test/... runs and plots BOBYQA PRAXIS SBPLX ... on these ndtestfuncs
    from several random startpoints, in several dimensions
Nd-testfuncs-python.md
F = Funcmon(func): wrap func() to monitor and plot F.fmem F.xmem F.cost
"""
    # zillions of papers and methods for derivative-free optimization alone


#...............................................................................
from __future__ import division
import numpy as np
from numpy import abs, cos, exp, mean, pi, prod, sin, sqrt, sum

try:
    from opt.testfuncs.powellsincos import Powellsincos
except ImportError:
    Powellsincos = None
try:
    from opt.testfuncs.randomquad import randomquad
    from opt.testfuncs.logsumexp import logsumexp
except ImportError:
    randomquad = logsumexp = None

__version__ = "2015-03-06 mar  denis-bz-py t-online de"  # + randomquad logsumexp


#...............................................................................
def ackley( x, a=20, b=0.2, c=2*pi ):
    x = np.asarray_chkfinite(x)  # ValueError if any NaN or Inf
    n = len(x)
    s1 = sum( x**2 )
    s2 = sum( cos( c * x ))
    return -a*exp( -b*sqrt( s1 / n )) - exp( s2 / n ) + a + exp(1)

#...............................................................................
def dixonprice( x ):  # dp.m
    x = np.asarray_chkfinite(x)
    n = len(x)
    j = np.arange( 2, n+1 )
    x2 = 2 * x**2
    return sum( j * (x2[1:] - x[:-1]) **2 ) + (x[0] - 1) **2

#...............................................................................
def griewank( x, fr=4000 ):
    x = np.asarray_chkfinite(x)
    n = len(x)
    j = np.arange( 1., n+1 )
    s = sum( x**2 )
    p = prod( cos( x / sqrt(j) ))
    return s/fr - p + 1

#...............................................................................
def levy( x ):
    x = np.asarray_chkfinite(x)
    n = len(x)
    z = 1 + (x - 1) / 4
    return (sin( pi * z[0] )**2
        + sum( (z[:-1] - 1)**2 * (1 + 10 * sin( pi * z[:-1] + 1 )**2 ))
        +       (z[-1] - 1)**2 * (1 + sin( 2 * pi * z[-1] )**2 ))

#...............................................................................
michalewicz_m = .5  # orig 10: ^20 => underflow

def michalewicz( x ):  # mich.m
    x = np.asarray_chkfinite(x)
    n = len(x)
    j = np.arange( 1., n+1 )
    return - sum( sin(x) * sin( j * x**2 / pi ) ** (2 * michalewicz_m) )

#...............................................................................
def perm( x, b=.5 ):
    x = np.asarray_chkfinite(x)
    n = len(x)
    j = np.arange( 1., n+1 )
    xbyj = np.fabs(x) / j
    return mean([ mean( (j**k + b) * (xbyj ** k - 1) ) **2
            for k in j/n ])
    # original overflows at n=100 --
    # return sum([ sum( (j**k + b) * ((x / j) ** k - 1) ) **2
    #       for k in j ])

#...............................................................................
def powell( x ):
    x = np.asarray_chkfinite(x)
    n = len(x)
    n4 = ((n + 3) // 4) * 4
    if n < n4:
        x = np.append( x, np.zeros( n4 - n ))
    x = x.reshape(( 4, -1 ))  # 4 rows: x[4i-3] [4i-2] [4i-1] [4i]
    f = np.empty_like( x )
    f[0] = x[0] + 10 * x[1]
    f[1] = sqrt(5) * (x[2] - x[3])
    f[2] = (x[1] - 2 * x[2]) **2
    f[3] = sqrt(10) * (x[0] - x[3]) **2
    return sum( f**2 )

#...............................................................................
def powersum( x, b=[8,18,44,114] ):  # power.m
    x = np.asarray_chkfinite(x)
    n = len(x)
    s = 0
    for k in range( 1, n+1 ):
        bk = b[ min( k - 1, len(b) - 1 )]  # ?
        s += (sum( x**k ) - bk) **2  # dim 10 huge, 100 overflows
    return s

#...............................................................................
def rastrigin( x ):  # rast.m
    x = np.asarray_chkfinite(x)
    n = len(x)
    return 10*n + sum( x**2 - 10 * cos( 2 * pi * x ))

#...............................................................................
def rosenbrock( x ):  # rosen.m
    """ http://en.wikipedia.org/wiki/Rosenbrock_function """
        # a sum of squares, so LevMar (scipy.optimize.leastsq) is pretty good
    x = np.asarray_chkfinite(x)
    x0 = x[:-1]
    x1 = x[1:]
    return (sum( (1 - x0) **2 )
        + 100 * sum( (x1 - x0**2) **2 ))

#...............................................................................
def schwefel( x ):  # schw.m
    x = np.asarray_chkfinite(x)
    n = len(x)
    return 418.9829*n - sum( x * sin( sqrt( abs( x ))))

#...............................................................................
def sphere( x ):
    x = np.asarray_chkfinite(x)
    return sum( x**2 )

#...............................................................................
def sum2( x ):
    x = np.asarray_chkfinite(x)
    n = len(x)
    j = np.arange( 1., n+1 )
    return sum( j * x**2 )

#...............................................................................
def trid( x ):
    x = np.asarray_chkfinite(x)
    return sum( (x - 1) **2 ) - sum( x[:-1] * x[1:] )

#...............................................................................
def zakharov( x ):  # zakh.m
    x = np.asarray_chkfinite(x)
    n = len(x)
    j = np.arange( 1., n+1 )
    s2 = sum( j * x ) / 2
    return sum( x**2 ) + s2**2 + s2**4

#...............................................................................
    # not in Hedar --

def ellipse( x ):
    x = np.asarray_chkfinite(x)
    return mean( (1 - x) **2 )  + 100 * mean( np.diff(x) **2 )

#...............................................................................
def nesterov( x ):
    """ Nesterov's nonsmooth Chebyshev-Rosenbrock function, Overton 2011 variant 2 """
    x = np.asarray_chkfinite(x)
    x0 = x[:-1]
    x1 = x[1:]
    return abs( 1 - x[0] ) / 4 \
        + sum( abs( x1 - 2*abs(x0) + 1 ))

#...............................................................................
def saddle( x ):
    x = np.asarray_chkfinite(x) - 1
    return np.mean( np.diff( x **2 )) \
        + .5 * np.mean( x **4 )


#-------------------------------------------------------------------------------
allfuncs = [
    ackley,
    dixonprice,
    ellipse,
    griewank,
    levy,
    michalewicz,  # min < 0
    nesterov,
    perm,
    powell,
    # powellsincos,  # many local mins
    powersum,
    rastrigin,
    rosenbrock,
    schwefel,  # many local mins
    sphere,
    saddle,
    sum2,
    trid,  # min < 0
    zakharov,
    ]

if Powellsincos is not None:  # try import
    _powellsincos = {}  # dim -> func
    def powellsincos( x ):
        x = np.asarray_chkfinite(x)
        n = len(x)
        if n not in _powellsincos:
            _powellsincos[n] = Powellsincos( dim=n )
        return _powellsincos[n]( x )

    allfuncs.append( powellsincos )

if randomquad is not None:  # try import
    allfuncs.append( randomquad )
    allfuncs.append( logsumexp )

allfuncs.sort( key = lambda f: f.__name__ )


#...............................................................................
allfuncnames = " ".join([ f.__name__ for f in allfuncs ])
name_to_func = { f.__name__ : f  for f in allfuncs }

    # bounds from Hedar, used for starting random_in_box too --
    # getbounds evals ["-dim", "dim"]
ackley._bounds       = [-15, 30]
dixonprice._bounds   = [-10, 10]
griewank._bounds     = [-600, 600]
levy._bounds         = [-10, 10]
michalewicz._bounds  = [0, pi]
perm._bounds         = ["-dim", "dim"]  # min at [1 2 .. n]
powell._bounds       = [-4, 5]  # min at tile [3 -1 0 1]
powersum._bounds     = [0, "dim"]  # 4d min at [1 2 3 4]
rastrigin._bounds    = [-5.12, 5.12]
rosenbrock._bounds   = [-2.4, 2.4]  # wikipedia
schwefel._bounds     = [-500, 500]
sphere._bounds       = [-5.12, 5.12]
sum2._bounds         = [-10, 10]
trid._bounds         = ["-dim**2", "dim**2"]  # fmin -50 6d, -200 10d
zakharov._bounds     = [-5, 10]

ellipse._bounds      =  [-2, 2]
logsumexp._bounds    = [-20, 20]  # ?
nesterov._bounds     = [-2, 2]
powellsincos._bounds = [ "-20*pi*dim", "20*pi*dim"]
randomquad._bounds   = [-10000, 10000]
saddle._bounds       = [-3, 3]


#...............................................................................
def getfuncs( names, dim=0 ):
    """ for f in getfuncs( "a b ..." ):
            y = f( x )
    """
    if names == "*":
        return allfuncs
    funclist = []
    for nm in names.split():
        if nm not in name_to_func:
            raise ValueError( "getfuncs( \"%s\" ) not found" % names )
        funclist.append( name_to_func[nm] )
    return funclist

def getbounds( funcname, dim ):
    """ "ackley" or ackley -> [-15, 30] """
    funcname = getattr( funcname, "__name__", funcname )
    func = getfuncs( funcname )[0]
    b = func._bounds[:]
    if isinstance( b[0], basestring ):  b[0] = eval( b[0] )
    if isinstance( b[1], basestring ):  b[1] = eval( b[1] )
    return b

#...............................................................................
_minus = "dixonprice perm powersum schwefel sphere sum2 trid zakharov "  # nlopt ~ same ?

def allfuncs_minus( minus=_minus ):
    return [f for f in allfuncs
            if f.__name__ not in minus.split() ]

def funcnames_minus( minus=_minus ):
    return " ".join([ f.__name__
            for f in allfuncs_minus( minus=minus )])


#-------------------------------------------------------------------------------
if __name__ == "__main__":  # standalone test --
    import sys

    dims = [2, 4, 10]  # , 100]
    nstep = 11  # 11: 0 .1 .2 .. 1
    seed = 0

        # to change these params in sh or ipython, run this.py  a=1  b=None  c=[3] ...
    for arg in sys.argv[1:]:
        exec( arg )

    np.set_printoptions( threshold=20, edgeitems=5, linewidth=120, suppress=True,
        formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    np.random.seed(seed)

    #...........................................................................
    for dim in dims:
        print "\n# ndtestfuncs dim %d  along the diagonal low .. high corner --" % dim
        # cmp matlab, anyone ?

        for func in allfuncs:
            lo, hi = getbounds( func, dim )
            steps = np.linspace( lo, hi, nstep )

            Y = np.array([ func( t * np.ones(dim) ) for t in steps ])
            jmin = Y.argmin()
            Ymin = Y[jmin]
            print "%-12s %dd  min %-8.3g  at %.3g \t lo hi %4.3g %4.3g \t Y-min %s" % (
                    func.__name__, dim, Ymin, steps[jmin], lo, hi, Y - Ymin )

    # plot-ndtestfuncs.py

Thank you so much ! You solved one of my big problems brother. I was using many of these mathematical optimization functions for my PhD thesis. Searched a lot of them but you made my day. Lots of love, grace and peace to you for writing all these for free and making them publically available for free

Great work man, saved me a big time

@Bessagg, you're welcome. Fwiw, there's an old plot comparing 3 DFOs (derivative-free) here on scicomp.stack (vote +- if you like);
there are certainly more papers than testcases on the web :) even in fractal subfields of optimization.