Run a multi-dimensional optimization

This page contains a very simple tutorial to optimize a multi-dimensional cost function. Since problems with larger domains, such as multi-dimensional problems, will take many iterations of optimizers to search and find the minima, we will demonsrate how to execute these searches in parallel to reduce the amount of time it takes.

The key concepts you should gain from this page:

• How to define a multi-dimensional cost function.

• How to use an optimization search to find the minimum of the multi-dimensional cost function.

• How to call the optimization in parallel.

Define a multi-dimensional cost function

We can increase the number of parameters in our cost function like using the 2-dimensional volcano function below.

import matplotlib.pyplot as plt
import mystic
import numpy
from mystic import models
from scipy import stats

def cost_func(p):
    """ Executed for each set of drawn parameters in the optimization search.
    """

    # get the x and y values from Mystic
    # p is a list so we map each index of the list to our variables (x, y)
    x, y = p

    # get value at Gaussian function x and y
    var = stats.multivariate_normal(mean=[0, 0], cov=[[0.5, 0], [0, 0.5]])
    gauss = -50.0 * var.pdf([x, y])

    # get value at volcano function x and y
    r = numpy.sqrt(x**2 + y**2)
    mu, sigma = 5.0, 1.0
    stat = 25.0 * (numpy.exp(-r / 35.0) + 1.0 /
                   (sigma * numpy.sqrt(2.0 * numpy.pi)) *
                    numpy.exp(-0.5 * ((r - mu) / sigma) ** 2)) + gauss

    # whether to flip sign of function
    # a positive lets you search for minimum
    # a negative lets you search for maximum
    stat *= 1.0

    return stat

mystic.model_plotter(cost_func, depth=True, fill=True, bounds='-9.5:9.5,-9.5:9.5')

Run an opitimizer

We can run a single optimizer with the multi-dimensional cost function like we did before. Below we set the initial points, strict ranges, and termination condition for a solver to find the minimum. This should look familar from the prior page in the tutorial.

from mystic import tools
from mystic.solvers import PowellDirectionalSolver
from mystic.ensemble import BuckshotSolver
from mystic.termination import VTR

# set random seed so we can reproduce results
tools.random_seed(0)

# create a solver
solver = PowellDirectionalSolver(dim=2)

# set the initial position to (0, 0)
solver.SetInitialPoints([0, 0])

# set the range to search for both parameters between -9.5 and 9.5
solver.SetStrictRanges((-9.5, -9.5), (9.5, 9.5))

# find the minimum
# pass the termination condition
solver.Solve(cost_func, termination=VTR())

# print the best parameters
print(f"The best solution is {solver.bestSolution}")

The best solution is [-2.24303552e-02  9.64325450e-07]

Run an ensemble of optimizers using multi-processing

In the prior examples, we only ran a single optimizer and found a nearby local minimum. However, we can use an ensemble of optimizers to map the minima of the cost function. Below, we show how to execute an optimization search using the BuckshotSolver. This is an ensemble solver available in the solver module. A short list of some ensemble sovlers are:

• BuckshotSolver: Ensemble from a uniform distribution of solvers in the parameter space.

• LatticeSolver: Ensemble from a grid of solvers in the parameter space.

Recall, we did use SetInitialPoint in prior examples. So there are other ways to set initial positions.

from mystic import tools
from mystic.solvers import PowellDirectionalSolver
from mystic.ensemble import BuckshotSolver
from mystic.termination import VTR

# set random seed so we can reproduce results
tools.random_seed(0)

# create a solver
solver = BuckshotSolver(dim=2, npts=8)

# since we have an search solver
# we specify what optimization algorithm to use within the search
solver.SetNestedSolver(PowellDirectionalSolver)

# set the range to search for both parameters between -9.5 and 9.5
solver.SetStrictRanges((-9.5, -9.5), (9.5, 9.5))

# find the minimum
# pass the termination condition
solver.Solve(cost_func, termination=VTR())

# print the best parameters
print(f"The best solution is {solver.bestSolution}")

The best solution is [-3.16087505e-05 -2.24303335e-02]

That took awhile, we can run things in parallel. Let’s cut the execution time by using multi-processing. We can specify a pool like below. Running the command below, you should see it took a shorter time to complete.

from pathos.pools import ProcessPool as Pool

# set random seed so we can reproduce results
tools.random_seed(0)

# create a solver
solver = BuckshotSolver(dim=2, npts=8)

# since we have an search solver
# we specify what optimization algorithm to use within the search
solver.SetNestedSolver(PowellDirectionalSolver)

# set multi-processing
solver.SetMapper(Pool().map)

# set the range to search for both parameters between -9.5 and 9.5
solver.SetStrictRanges((-9.5, -9.5), (9.5, 9.5))

# find the minimum
# pass the termination condition
solver.Solve(cost_func, termination=VTR())

# print the best parameters
print(f"The best solution is {solver.bestSolution}")

The best solution is [-3.16087505e-05 -2.24303335e-02]

Sample and store cached results

We can store results to an archive of sampled points. Below is an example of creating an archive, using a sampler, and then plotting the results.

A list of available samplers are in the samplers modules. These samplers are Mystic’s interface to its optimizer-driven sampling interface. An example using the LatticeSampler is shown below. Its execution is similar to the LatticeSolver, however, we can connect the archive to it for storing the sampled points.

from mystic import cache
from mystic.math import interpolate
from mystic.samplers import LatticeSampler

# create an archive
model = cache.cached(archive="volcano")(cost_func)

# sample initial points from a Uniform distribution then search for minima
sampler = LatticeSampler([(-9.5, 9.5), (-9.5, 9.5)], model, npts=8)
sampler.sample_until(terminated=all)

# invert the model
imodel = model.__inverse__

# sample initial points from a Uniform distribution then search for maxima
isampler = LatticeSampler([(-9.5, 9.5), (-9.5, 9.5)], imodel, npts=8)
isampler.sample_until(terminated=all)

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import mystic.cache as mc

# create a figure
figure = plt.figure()
ax = plt.axes(projection="3d")
ax.autoscale(tight=True)

# read the archive and show the sampled points
c = mc.archive.read('volcano')
x, y = np.array(list(c.keys())), np.array(list(c.values()))
ax.plot(x.T[0], x.T[1], y, 'ko', linewidth=2, markersize=4)
plt.show()

import pox
pox.rmtree("volcano")

We can increase the number of samples and get a better interpolation of the cost function.

for i, n in enumerate([4, 8, 16]):

    # create an archive
    model = cache.cached(archive=f"volcano_{n}")(cost_func)

    # sample initial points from a Uniform distribution then search for minima
    sampler = LatticeSampler([(-9.5, 9.5), (-9.5, 9.5)], model, npts=n)
    sampler.sample_until(terminated=all)

    # invert the model
    imodel = model.__inverse__

    # sample initial points from a Uniform distribution then search for maxima
    isampler = LatticeSampler([(-9.5, 9.5), (-9.5, 9.5)], imodel, npts=n)
    isampler.sample_until(terminated=all)

    # create a figure
    figure = plt.figure(i)
    ax = plt.axes(projection="3d")
    ax.autoscale(tight=True)

    # read the archive and show the sampled points
    c = mc.archive.read(f"volcano_{n}")
    x, y = np.array(list(c.keys())), np.array(list(c.values()))
    ax.plot(x.T[0], x.T[1], y, 'ko', linewidth=2, markersize=4)

    pox.rmtree(f"volcano_{n}")

plt.show()

You should now have an understanding of how to setup a optimization search and execute it in parallel using a pool. We have also shown how to attach an archive to a sampler and plot the results.

On the next page, we will discuss how we can use a surrogate model to find the extrema of expensive cost functions using a surrogate model.