Toy problems

The following subsections show how to use OMADS on selected 2D test functions.

Tip

The videos below show the poll trajectory, frame size and mesh coarseness during an optimization. They were recorded with an earlier release of OMADS, so the path differs from what the current version takes, but the mechanics are the same.

2D Rosenbrock

\[f(x_{1}, x_{2})=(1-x_{1})^{2}+100(x_{2}-x_{1}^{2})^{2}.\]

The Rosenbrock function has curved contours, and its global minimum \(f(1, 1) = 0\) lies in a parabolic, banana-shaped valley. Many algorithms find the valley easily but struggle to reach the minimum, because the function values in the valley barely change.

The script below solves it with each of the three OMADS solvers:

Listing 7 rosenbrock.py

"""2D Rosenbrock with the three OMADS solvers."""
import numpy as np
from OMADS import mads, poll, search


def rosenbrock(x):
  x = np.asarray(x)
  f = np.sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)
  return [f, [0]]


data = {
    "evaluator": {"blackbox": rosenbrock},
    "param": {
        "name": "rosen",
        "baseline": [-2.0, -2.0],
        "lb": [-5, -5],
        "ub": [10, 10],
        "var_names": ["x1", "x2"],
        "scaling": [15.0, 15.0],
        "post_dir": "./post",
    },
    "options": {
        "seed": 0, "budget": 1000, "tol": 1e-9,
        "check_cache": True, "store_cache": True, "rich_direction": True,
    },
    "search": {"type": "sampling", "s_method": "ACTIVE", "ns": 10},
}

if __name__ == "__main__":
  for solver in (mads, poll, search):
    out = solver.main(data)[0]
    print(f"{solver.__name__.split('.')[-1]:>6}: xmin = {out['xmin']}, "
          f"fmin = {out['fmin'][0]:.3e}, iterations = {out['niterations']}")

  mads: xmin = [1.0, 1.0], fmin = 0.000e+00, iterations = 21
  poll: xmin = [0.9390255808830261, 0.8823902606964111], fmin = 3.756e-03, iterations = 253
search: xmin = [-0.602294921875, 0.409423828125], fmin = 2.785e+00, iterations = 169

Full MADS reaches the exact minimum in 223 evaluations. Its quadratic model search makes large, well-aimed steps along the valley. The poll alone crawls along the valley with orthogonal steps and is still on its way when the 1000-evaluation budget runs out, and the sampling search alone is an exploration tool, not a local solver. Use mads unless you have a reason not to.

2D Alpine

\[f({\bf x})= \sum_{i=1}^{d} |x_{i} \sin(x_{i})+0.1x_{i}|.\]

The Alpine function is non-smooth, which can stop gradient-based algorithms from reaching the global minimum. Its design space is non-smooth, non-convex, separable and multimodal, and it has several global minima with \(f = 0\), for example at the origin and near \(x_i \approx 6.18\).

Listing 8 alpine.py

"""2D Alpine function (non-smooth, multimodal)."""
import numpy as np
from OMADS import mads


def alpine(x):
  f = sum(abs(xi * np.sin(xi) + 0.1 * xi) for xi in x)
  return [f, [0]]


data = {
    "evaluator": {"blackbox": alpine},
    "param": {"name": "alpine", "baseline": [5.0, 5.0], "lb": [0.0, 0.0],
              "ub": [10.0, 10.0], "scaling": [10.0, 10.0], "post_dir": "./post"},
    "options": {"seed": 0, "budget": 1000, "tol": 1e-9,
                "check_cache": True, "store_cache": True, "rich_direction": True},
    "search": {"type": "sampling", "s_method": "ACTIVE", "ns": 10},
}

if __name__ == "__main__":
  out, *_ = mads.main(data)
  print("xmin =", out["xmin"])
  print("fmin =", out["fmin"], "after", out["nbb_evals"], "evaluations")

xmin = [6.183017810806632, 6.183017225936055]
fmin = [4.5235426833212244e-06] after 270 evaluations

2D Ackley 3

Ackley 3 combines a large funnel with a bumpy exponential term, which can trap hill-climbing algorithms in one of its many sub-optimal solutions:

\[f(x_{1}, x_{2})= -200\, e^{-0.2 \sqrt{x^{2}_{1} + x^{2}_{2}}} + 5\, e^{\cos(3x_{1})+\sin(3x_{2})}.\]

Its design space is non-convex, non-separable and multimodal. The global minimum is \(f \approx -186.4112\) near \((0, \pm 0.0068)\).

Listing 9 ackley3.py

"""2D Ackley 3 function."""
import numpy as np
from OMADS import mads


def ackley3(x):
  f = (-200.0 * np.exp(-0.2 * np.sqrt(x[0] ** 2 + x[1] ** 2))
       + 5.0 * np.exp(np.cos(3.0 * x[0]) + np.sin(3.0 * x[1])))
  return [f, [0]]


data = {
    "evaluator": {"blackbox": ackley3},
    "param": {"name": "ackley3", "baseline": [5.0, 5.0], "lb": [-32.0, -32.0],
              "ub": [32.0, 32.0], "scaling": [64.0, 64.0], "post_dir": "./post"},
    "options": {"seed": 0, "budget": 1000, "tol": 1e-9,
                "check_cache": True, "store_cache": True, "rich_direction": True},
    "search": {"type": "sampling", "s_method": "ACTIVE", "ns": 10},
}

if __name__ == "__main__":
  out, *_ = mads.main(data)
  print("xmin =", out["xmin"])
  print("fmin =", out["fmin"], "after", out["nbb_evals"], "evaluations")

xmin = [4.852190613746643e-06, -0.006760597229003906]
fmin = [-186.41121263368007] after 253 evaluations

2D Eggholder

The Eggholder function is a classic benchmark for global optimizers. Its many deep, irregular basins make the single global minimum, \(f(512, 404.2319) \approx -959.6407\), hard to find:

\[f(x_{1}, x_{2})= -(x_{2}+47)\sin\left(\sqrt{\left|\tfrac{x_{1}}{2}+x_{2}+47\right|}\right)-x_{1}\sin\left(\sqrt{|x_{1}-(x_{2}+47)|}\right).\]

This example runs MADS twice, with the default ACTIVE search and with a space-filling Latin hypercube search of 30 points per step:

Listing 10 eggholder.py

"""2D Eggholder function: compare an ACTIVE and a Latin hypercube search."""
import copy
import numpy as np
from OMADS import mads


def eggholder(x):
  a, b = x
  f = (-(b + 47.0) * np.sin(np.sqrt(abs(a / 2.0 + (b + 47.0))))
       - a * np.sin(np.sqrt(abs(a - (b + 47.0)))))
  return [f, [0]]


data = {
    "evaluator": {"blackbox": eggholder},
    "param": {"name": "eggholder", "baseline": [50.0, 280.0],
              "lb": [-512.0, -512.0], "ub": [512.0, 512.0],
              "scaling": [1024.0, 1024.0], "post_dir": "./post"},
    "options": {"seed": 0, "budget": 2000, "tol": 1e-9,
                "check_cache": True, "store_cache": True, "rich_direction": True},
    "search": {"type": "sampling", "s_method": "ACTIVE", "ns": 10},
}

if __name__ == "__main__":
  out, *_ = mads.main(data)
  print(f"ACTIVE search: fmin = {out['fmin'][0]:.4f} at {out['xmin']}")

  lhs = copy.deepcopy(data)
  lhs["search"] = {"type": "sampling", "s_method": "LH", "ns": 30}
  out, *_ = mads.main(lhs)
  print(f"LH search:     fmin = {out['fmin'][0]:.4f} at {out['xmin']}")

ACTIVE search: fmin = -493.9605 at [150.2320980578661, 301.3137746155262]
LH search:     fmin = -959.6407 at [512.0, 404.2317980080843]

The ACTIVE search concentrates on the neighbourhood of the incumbent, which is efficient on most problems but converges here to a local basin. The Latin hypercube search keeps sampling the whole box, lands in the right basin and the poll finishes the job on the bound \(x_1 = 512\).

Choosing the search for multimodal functions

On strongly multimodal functions, try "s_method": "LH" with ns between \(5n\) and \(20n\), try several seeds, or start from several points by giving baseline as a list of points.