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

@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.