Visualize n-dimensional bayesian optimization results

Question

I am working on a 6-dimensional bayesian optimization problem using (skopt's gp_minimize). After the optimizer ran for j iterations I would like to somehow visualize the "progress/result" of the optimization. As I am new to Bayesian optimization I would like to ask for input on how and what to visualize. What are good parameters to visualize to show the improvement and maybe even the parameter dependency of the optimized parameters?

Answer 1

One way to do it is to look at the "trajectory" of the best point as a function of the number of functions evaluations. When the objective function is 2D, you can actually look at contours of the objective functions and plot the best point trajectory at the same time, just like I did here:

http://infinity77.net/go_2021/thebenchmarks.html#solver-solvers-example

For higher dimensionalities, you can try and track a one-dimensional evolution of the objective function and its parameters. Example below with N=3:

import numpy
import matplotlib.pyplot as plt

from scipy.optimize import dual_annealing


def rosen(x, function_history, params_history):
    """
    Rosenbrock objective function.

    This class defines the Rosenbrock global optimization problem. This is a
    minimization problem defined as follows:

    .. math::

       f_{\\text{Rosenbrock}}(x) = \\sum_{i=1}^{n-1} [100(x_i^2 - x_{i+1})^2 + (x_i - 1)^2]


    Here, :math:`n` represents the number of dimensions and
    :math:`x_i \\in [-5, 10]` for :math:`i = 1, ..., n`.

    *Global optimum*: :math:`f(x) = 0` for :math:`x_i = 1` for :math:`i = 1, ..., n`

    Jamil, M. & Yang, X.-S. A Literature Survey of Benchmark Functions
    For Global Optimization Problems Int. Journal of Mathematical Modelling
    and Numerical Optimisation, 2013, 4, 150-194.
    """

    f = numpy.sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)

    params_history.append(x)
    function_history.append(f)

    return f


def solve():

    N = 3
    numpy.random.seed(1)
    bounds = list(zip([-5.] * N, [10.0] * N))
    minimizer_kwargs = {'method': 'L-BFGS-B', 'options': {'maxiter': 100, 'disp': False, 'ftol': 1e-6},
                        'bounds': bounds}

    x0 = numpy.random.uniform(-5.0, 10.0, N)
    function_history, params_history = [], []
    res = dual_annealing(rosen, bounds, x0=x0, seed=1, local_search_options=minimizer_kwargs,
                         args=(function_history, params_history))

    fig = plt.figure()
    ax = fig.add_subplot(2, 2, 1)

    start_f = function_history[0]
    nfun = [1]
    f_progression = [start_f]
    params_progression = [params_history[0]]

    for k, f in enumerate(function_history):
        if f < start_f:
            nfun.append(k+1)
            f_progression.append(f)
            params_progression.append(params_history[k])
            start_f = f

    ax.plot(nfun, f_progression, lw=2)
    ax.set_ylabel('$F(X_1, X_2, X_3)$', fontsize=16)
    ax.set_yscale('log')
    params_progression = numpy.array(params_progression)

    for i in range(3):
        ax = fig.add_subplot(2, 2, i+2)
        ax.plot(nfun, params_progression[:, i], lw=2)
        ax.set_ylabel('$X_%d$'%(i+1), fontsize=16)

    for ax in fig.axes:
        ax.grid()
        ax.set_xlabel('Functions Evaluations', fontsize=16)

    plt.show()


if __name__ == '__main__':
    solve()

You get a picture like the one below:

Not perfect, but one way to do it :-) .

Good luck.