Benchmarking different studies against each other

Drivers:

BayesianOptimization, ActiveLearning, CMAES, DifferentialEvolution, ScipyMinimizer

Download script: benchmark.m

In some cases, it is not clear, which driver or which driver configuration is able to minimize an objective function with the smallest number of iterations. If similar minimizations are performed on a regular basis, it can be worthwhile to find the best choice by benchmarking several studies with different drivers or configurations against each other.

As an example, the usual 2D Rastrigin function on a circular domain shall be minimized:

\[ \begin{align}\begin{aligned}&\text{min.}\,& f(x_1,x_2) = 2\cdot10 + \sum_{i=1,2} \left(x_i^2 - 10\cos(2\pi x_i)\right)\\&\text{s.t.}\,& \sqrt{x_1^2 + x_2^2} \leq 1.5.\end{aligned}\end{align} \]

We compare the heuristic minimization methods CMAES and DifferentialEvolution, the multi-start local minimization method ScipyMinimizer, and Bayesian optimization. For the latter, we use the standard BayesianOptimization driver and the ActiveLearning driver that is configured to use Neural Networks instead of a Gaussian process to learn the objective function.

We assume that we can evaluate the function 2 times in parallel and that each evaluation takes 2 seconds. Hence, the driver has to compute new samples in less than a second in order to avoid significant overhead. This is an edge case for the Bayesian optimization approach, that has a considerable overhead typically larger than a second. I this case the other methods with much faster sample computation times can be more appropriate. Nevertheless, in this case of a multi-modal objective Bayesian optimization reaches the global minim still well before the other optimization methods. Gaussian process regression is computational efficient for this small number of iterations. In comparison the Neural Network training time of the ActiveLearning driver leads to a larger overhead.

jcm_optimizer_path = '<JCM_OPTIMIZER_PATH>';
addpath(fullfile(jcm_optimizer_path, 'interface', 'matlab'));
 
server = jcmoptimizer.Server();
client = jcmoptimizer.Client(server.port); 

% Definition of the search domain
design_space = {...
   struct('name', 'x1', 'type', 'continuous', 'domain', [-1.5,1.5]), ...
   struct('name', 'x2', 'type', 'continuous', 'domain', [-1.5,1.5]), ...
};

% Definition of fixed environment parameter
environment = {...     
   struct('name', 'radius', 'type', 'fixed', 'domain', 1.5) ...
};

% Definition of a constraint on the search domain
constraints = {...
   struct('name', 'circle', 'expression', 'sqrt(x1^2 + x2^2) <= radius')...
};

%Creation of studies to benchmark against each other
studies = struct();
drivers = {"BayesianOptimization", "ActiveLearning", "CMAES", ...
           "DifferentialEvolution", "ScipyMinimizer"};
for i = 1:length(drivers)
   studies.(drivers{i}) = client.create_study( ...
        'design_space', design_space, ...
        'environment', environment, ...
        'constraints', constraints, ...
        'driver', drivers{i}, ...
        'name', drivers{i}, ...
        'study_id', sprintf("benchmark_%s", drivers{i}), ...
        'open_browser' , false ...
   );
end

%Configuration of all studies
config = {'num_parallel', 1, 'max_iter', 250};
min_val = 1e-3; %Stop study when this value was observed
for i = 1:length(drivers)
    driver = drivers{i};
    study = studies.(driver);
    if strcmp(driver, "ActiveLearning")
        study.configure( ...
           'surrogates', {struct('type', "NN")}, ...
           'objectives', {struct('type', "Minimizer", 'min_val', min_val)}, ...
           config{:});
    elseif strcmp(driver, "ScipyMinimizer")
        %For a more global search, 3 initial gradient-free Nelder-Mead are started
        study.configure( ...
           'method', "Nelder-Mead", 'num_initial', 3, ...
           'min_val', min_val, config{:});
    else
        study.configure('min_val', min_val, config{:});
    end
end
% Evaluation of the black-box function for specified design parameters
function observation = evaluate(study, sample)

    pause(2); % make objective expensive
    observation = study.new_observation();
    x1 = sample.x1;
    x2 = sample.x2;
    observation.add(10*2 ...
                    + (x1.^2-10*cos(2*pi*x1)) ...
                    + (x2.^2-10*cos(2*pi*x2)) ...
                   );
end  

% Creation of a benchmark with 6 repetitions and add all 4 studies
benchmark = client.create_benchmark('num_average', 6);
for i = 1:length(drivers)
   benchmark.add_study(studies.(drivers{i}));
end

% Run the benchmark - this will take a while
benchmark.set_evaluator(@evaluate);
benchmark.run();

% Plot cummin convergence w.r.t. number of evaluations and time
fig = figure('Position', [0, 0, 1000, 500]);
subplot(1, 2, 1) 
data = benchmark.get_data('x_type', "num_evaluations");

for i = 1:length(data.names)
  std_error = data.sdev(i,:) ./ sqrt(6);
  p = plot(data.X(i,:), data.Y(i,:), 'LineWidth', 2.0);
  hold on
  patch([data.X(i,:), fliplr(data.X(i,:))], ...
        [data.Y(i,:) - std_error, fliplr(data.Y(i,:) + std_error)], ...
        p(1).Color, 'FaceAlpha', 0.2, 'EdgeAlpha', 0.2);
end
grid()
xlabel("Number of Evaluations", 'FontSize', 12)
ylabel("Average Cummulative Minimum", 'FontSize', 12)
ylim([-0.4, 10])

subplot(1, 2, 2)
data = benchmark.get_data('x_type', "time");
plots=[];
for i = 1:length(data.names)
  std_error = data.sdev(i,:) ./ sqrt(6);
  p = plot(data.X(i,:), data.Y(i,:), 'LineWidth', 2.0, 'DisplayName', data.names{i});
  plots(i) = p;
  hold on
  patch([data.X(i,:), fliplr(data.X(i,:))], ...
        [data.Y(i,:) - std_error, fliplr(data.Y(i,:) + std_error)], ...
        p(1).Color, 'FaceAlpha', 0.2, 'EdgeAlpha', 0.2);
end
legend(plots, data.names{:});
grid()
xlabel("Time (sec)", 'FontSize', 12)
ylim([-0.4, 10])
saveas(fig, "benchmark.png")

Left: The left graph shows the cummulative minimum as a function of the number of evaluations of the objective function. The results are averaged over six independent runs of each driver. The shaded area shows the uncertainty of the mean value. Clearly, the BayesianOptimization driver performs the best and finds the global minimum zero after about 70 iterations while the other methods did not converge to zero even after 250 iterations. Right: The convergence of the cummulative minimum as a function of the total optimization time looks slightly different. As expected, the ActiveLearning and BayesianOptimization drivers show an overhead for computing new samples. If one would have an optimization budget on less than 50 seconds, the multi-start local ScipyMinimizer performs better. It can be also seen that the training of Neural Networks leads to a larger overhead compared to the Gaussian process regression used by the BayesianOptimization driver.