Differential Evolution › Implementation in NumPy and SciPy
Machine Learning from ScratchAll chapters Differential Evolution › Implementation in NumPy and SciPySection view · Sections in this chapter6 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
The primitive implementation uses NumPy arrays for the population and objective values. We compare its result with scipy.optimize.differential_evolution, then evaluate it alongside uniform random search and bounded multistart Nelder-Mead under the stated budget.
The two-dimensional Rastrigin objective adds a cosine ripple to a quadratic bowl:
where:
- : the number of decision variables
- : the -th coordinate of
Rewrite the objective as . Every summand is nonnegative, and zero is attained only at the origin, which lies inside our box. This proves the known global minimum used to interpret this synthetic example. The ripples create local minima with different objective values. This smooth benchmark illustrates search behavior; it does not establish performance on arbitrary black-box or non-differentiable objectives.
PythonIn [9]Expand
def rastrigin(X):
"""Vectorized Rastrigin function; accepts a (n, d) array and returns (n,) values."""
X = np.atleast_2d(X)
return np.sum(X**2 + 10.0 * (1.0 - np.cos(2.0 * np.pi * X)), axis=1)
BOUNDS = np.array([[-5.12, 5.12], [-5.12, 5.12]])
def de_rand_1_bin(
f, bounds, popsize=20, F=0.7, CR=0.9, max_nfev=3000, seed=0, record=False
):
"""Synchronous DE/rand/1/bin; count candidate evaluations, including initialization."""
bounds = np.asarray(bounds, dtype=float)
if (
bounds.ndim != 2
or bounds.shape[1] != 2
or not np.isfinite(bounds).all()
):
raise ValueError("bounds must be a finite (dimension, 2) array")
lo, hi = bounds[:, 0], bounds[:, 1]
if lo.size == 0 or np.any(lo >= hi):
raise ValueError("each lower bound must be less than its upper bound")
if not isinstance(popsize, (int, np.integer)) or popsize < 4:
raise ValueError("DE/rand/1 requires at least four population members")
if not isinstance(max_nfev, (int, np.integer)) or max_nfev < popsize:
raise ValueError("the budget must cover the initial population")
if not 0 <= F < 2 or not 0 <= CR <= 1:
raise ValueError("require 0 <= F < 2 and 0 <= CR <= 1")
rng = np.random.default_rng(seed)
d = lo.size
pop = rng.uniform(lo, hi, size=(popsize, d))
fit = np.asarray(f(pop), dtype=float)
if fit.shape != (popsize,) or not np.isfinite(fit).all():
raise ValueError(
"the batch objective must return one finite value per member"
)
nfev = popsize
best = int(np.argmin(fit))
best_x, best_f = pop[best].copy(), float(fit[best])
trace_best = np.minimum.accumulate(fit).tolist()
snaps = [pop.copy()] if record else None
snap_evaluations = [nfev] if record else None
while nfev < max_nfev:
# Every donor and target comes from the same frozen generation.
old_pop, old_fit = pop.copy(), fit.copy()
next_pop, next_fit = old_pop.copy(), old_fit.copy()
for i in range(popsize):
if nfev >= max_nfev:
break # A partial final generation spends exactly the remaining budget.
pool = np.delete(np.arange(popsize), i)
r1, r2, r3 = rng.choice(pool, size=3, replace=False)
mutant = np.clip(
old_pop[r1] + F * (old_pop[r2] - old_pop[r3]), lo, hi
)
mask = rng.random(d) < CR
mask[rng.integers(d)] = True
trial = np.where(mask, mutant, old_pop[i])
trial_values = np.asarray(f(trial[None, :]), dtype=float)
if trial_values.shape != (1,) or not np.isfinite(trial_values[0]):
raise ValueError(
"the objective must return one finite trial value"
)
ft = float(trial_values[0])
nfev += 1 # Rejected and numerically duplicate trials also cost one evaluation.
if ft <= old_fit[i]:
next_pop[i], next_fit[i] = trial, ft
if ft < best_f:
best_f, best_x = ft, trial.copy()
trace_best.append(best_f)
pop, fit = next_pop, next_fit
if record:
snaps.append(pop.copy())
snap_evaluations.append(nfev)
return dict(
best_x=best_x,
best_f=best_f,
nfev=nfev,
trace_nfev=np.arange(1, nfev + 1),
trace_best=np.asarray(trace_best),
snaps=snaps,
snap_evaluations=snap_evaluations,
population=pop,
)The old population is copied before each sweep, so every trial uses the same generation of donors and targets. Every initialized individual and every trial contributes one candidate evaluation, including ties, repeats, and rejected trials. The final sweep can be partial: a budget of 21 with a population of 20 permits one trial, not a second full population. The trace records the incumbent after every candidate, including initialization. We retain the best candidate separately from population membership.
Now we run it, and also run the SciPy reference with a comparable budget.
ConsoleOut [11]Expand
our DE/rand/1/bin : best f = 3.346e-15 at [ 0. -0.] using 3000 evaluations scipy reference : best f = 9.950e-01 at [0. 0.995] using 3000 evaluations difference in best value: 9.95e-01
The reference uses the same mutation factor, crossover rate, uniform initialization scheme, and deferred convention, but its random draws need not match our implementation. In this two-variable setup, SciPy's popsize=10 means 20 members; 149 generations plus initialization allow at most scalar evaluations when polishing is disabled. Zero tolerances can still stop at exactly zero population-energy spread. We display the actual counts and results. A single agreement or disagreement cannot prove an implementation correct; the worked updates, budget tests, and feasibility checks are separate checks. These API conventions are documented in the version-matched SciPy source.
Key Parameters
These are the choices used in the chapter. Their effects depend on the objective and the evaluation budget:
popsize: our function takes the actual number of members, at least four. Initialization spends that many evaluations; a larger population leaves fewer generations under a fixed budget. SciPy instead uses a dimension-dependent multiplier unless an explicit initial population is supplied.F: scales a donor difference. We use 0.7 and test other values below. A fixed does not imply a fixed absolute step length. SciPy's mutation interval samples one value per generation, not per mutant.CR: controls non-forced coordinate inheritance. We use 0.9, while forcing one selected coordinate at every setting. The sweep is an example rather than a universal tuning rule.max_nfev: includes all initialization and trial evaluations. Our function rejects a budget smaller than the initial population and stops partway through a sweep when necessary.strategy: this chapter implementsrand1bin; alternatives change donor or crossover rules and require their own evaluation.
Key Methods
The reference and teaching functions have different purposes:
differential_evolution(): Runs the full DE loop, including initialization, mutation, crossover, and selection. Use it when a reference implementation or built-in bound and polishing options are needed.de_rand_1_bin(): The chapter's NumPy implementation, useful for inspecting the four ingredients and for budget-matched experiments.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!