Benchmark against scipy

As biteopt is implemented in C++ its inner loops can be expected to run a lot faster than scipy’s solvers which are coded in Python. Here we will benchmark biteopt against scipy’s most powerful stochastic optimizers on a famous test function for optimization which can be computed very quickly: the Rastrigin function.

import numpy as np
from scipybiteopt import biteopt
from scipy.optimize import dual_annealing, differential_evolution
import timeit

The rules are: No local optimizations allowed (therefore no_local_search=True for dual_annealing and polish = False for differential_evolution) as we want to compare the performance of the stochastic optimizers. We keep the default number of function evaluations ad termination criteria. We use a moderately large number of dimensions: d = 10.

def rastrigin(x):
    return np.sum(x*x - 10*np.cos(2*np.pi*x)) + 10*x.shape[0]

lower_bounds = [-8.12] * 10
upper_bounds = [5.12] * 10
bounds=list(zip(lower_bounds, upper_bounds))

res_da = dual_annealing(rastrigin, bounds, no_local_search=True)
print("Scipy Dual Annealing: minimal value={}, number of function evaluations={}".format(res_da.fun, res_da.nfev))
%timeit -n 1 res_da = dual_annealing(rastrigin, bounds, no_local_search=True)

print()
res_de = differential_evolution(rastrigin, bounds, polish = False)
print("Scipy Differential Evolution: minimal value={}, number of function evaluations={}".format(res_de.fun, res_de.nfev))
%timeit -n 1 res_de = differential_evolution(rastrigin, bounds, polish = False)

print()
res_bo = biteopt(rastrigin, bounds)
print("Biteopt: minimal value={}, number of function evaluations={}".format(res_bo.fun, res_bo.nfev))
%timeit -n 1 res_bo = biteopt(rastrigin, bounds)

Scipy Dual Annealing: minimal value=1.2482702970828541e-05, number of function evaluations=20001
953 ms ± 11.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Scipy Differential Evolution: minimal value=0.0, number of function evaluations=124800
4.05 s ± 825 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Biteopt: minimal value=4.5713477447861806e-10, number of function evaluations=10000
86 ms ± 2.72 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

We see that biteopt is app. 10 times faster than dual_annealing and 50 times faster than differential-evolution.

Let’s see how this changes when the objective function is also compiled and no longer uses the Python interpreter. For that we will use numba’s convenient just in time compiler.

from numba import njit

#let numba JIT compile
rastrigin_compiled = njit(rastrigin)
out = rastrigin_compiled(np.zeros((10, )))

res_da = dual_annealing(rastrigin_compiled, bounds, no_local_search=True)
print("Scipy Dual Annealing: minimal value={}, number of function evaluations={}".format(res_da.fun, res_da.nfev))
%timeit -n 1 res_da = dual_annealing(rastrigin_compiled, bounds, no_local_search=True)

print()
res_de = differential_evolution(rastrigin_compiled, bounds, polish = False)
print("Scipy Differential Evolution: minimal value={}, number of function evaluations={}".format(res_de.fun, res_de.nfev))
%timeit -n 1 res_de = differential_evolution(rastrigin_compiled, bounds, polish = False)

print()
res_bo = biteopt(rastrigin_compiled, bounds)
print("Biteopt: minimal value={}, number of function evaluations={}".format(res_bo.fun, res_bo.nfev))
%timeit -n 1 res_bo = biteopt(rastrigin_compiled, bounds)

Scipy Dual Annealing: minimal value=8.391295509113661e-06, number of function evaluations=20001
734 ms ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Scipy Differential Evolution: minimal value=0.0, number of function evaluations=119700
4.2 s ± 972 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Biteopt: minimal value=0.0, number of function evaluations=10000
11.7 ms ± 346 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Now biteopt is app. 70 times faster than dual_annealing and 400 times faster than differential-evolution!

Note that this comes with a tradeoff though: biteopt does not perform any sanity checks. Using scipybiteopt can be like using raw C/C++: Errors or exceptions occuring during evaluation of the objective are not caught and might crash the Python interpreter. Help is always welcome to make this wrapper more robust.