Differential Evolution

Why Differential Evolution Matters

Suppose you must tune two real-valued settings for a system whose quality can be measured but whose derivatives are unavailable. Trying a dense grid becomes expensive as the number of settings grows, while a single local search may explore only one region. Differential evolution, or DE, maintains several candidate vectors and uses differences between them to construct new proposals at a scale suggested by the current population.

The distinctive idea is simple: take a base vector, add a scaled difference of two other vectors, mix the result with the current target, and keep the trial only when its objective is no worse. Here we minimize. The algorithm can explore a bounded continuous domain without a gradient, but a finite run normally returns the best point it evaluated rather than an optimality certificate.

We will implement DE/rand/1/bin with a supplied population, three distinct donor indices excluding the target, coordinate clipping, forced binomial crossover, and deferred generation replacement. The worked function is a shifted sphere with a known answer, allowing independent verification. Python, C++, and Java will produce identical sample output and the same recorded numerical decisions.

Core Concepts and Terminology

Let each candidate be a real vector $\mathbf x\in\mathbb R^d$ inside a box with bounds $\ell_j\le x_j\le h_j$. The population $P_g=(\mathbf x_{0,g},\ldots,\mathbf x_{N-1,g})$ is an ordered list of $N$ such vectors, and $f(\mathbf x)$ is a deterministic objective to minimize. The vectors may coincide, but their indices remain distinct.

For target $i$, choose ordered indices $a,b,c$ that are mutually distinct and different from $i$. Thus this variant requires $N\ge4$. Its raw mutant or donor vector is

$$\mathbf v_{i,g}=\mathbf x_{a,g}+F(\mathbf x_{b,g}-\mathbf x_{c,g}),\qquad F\gt0.$$

“rand” means the base index $a$ is sampled rather than fixed at the current best. “1” means one difference vector is used; it does not mean one-dimensional optimization. “bin” means coordinate-wise, binomial-style crossover. The scale $F$ is not a mutation probability: the unbounded differential displacement has length $F\lVert\mathbf x_b-\mathbf x_c\rVert_2$. Clipping and crossover can subsequently change its direction and length relative to the target.

The worked objective is

$$f(x,y)=(x-1)^2+(y+1)^2,\qquad -5\le x,y\le5.$$

Its unique real-valued minimum is $f(1,-1)=0$. The dashed path in the diagram connects successive best discoveries, which can belong to different population slots; DE does not have a single search trajectory or a fixed step size. The two-generation run stops at objective $2$, so it has improved substantially without reaching the known optimum.

How Differential Evolution Works, Step by Step

The input fixes the initial population, lower and upper bounds, objective center, seed, generation limit $G$, differential scale $F$, crossover threshold $m$, and evaluation budget $B$. The programs use the deterministic shifted-sphere objective

$$f(\mathbf x)=\sum_{j=0}^{d-1}(x_j-c_j)^2.$$

1. Validate and initialize. Copy all input vectors. Require enough budget to evaluate all $N$ initial members; otherwise return EVAL_LIMIT with no evaluated point and no random draws. Evaluate every initial member, even duplicates, and archive the earliest minimum. Record generation $0$.
2. Check limits before a generation. If $G$ generations have completed, return GEN_LIMIT. Otherwise require budget for $N$ trial evaluations; if it is insufficient, return EVAL_LIMIT before drawing randomness. No partial generation is produced.
3. Choose three donors without replacement. Exclude the target index, then draw a rank among the remaining indices. Remove that chosen index and repeat twice. The code maps each rank past at most three sorted exclusions, avoiding repeated whole-population shuffling or rejection of duplicate indices. Under uniform rank draws, this samples ordered donor triples uniformly.
4. Draw the crossover choices. Draw a forced coordinate $j_{\mathrm{rand}}$ in $0,\ldots,d-1$, then one integer $r_j$ in $0,\ldots,999$ for every coordinate. These draws are consumed even for the forced coordinate and at crossover rates $0$ and $1$.
5. Form and repair the mutant. Compute each difference, multiply by $F$, and add the base coordinate in that order. Clip every raw mutant coordinate into its interval:

$$\widehat v_j=\min\bigl(h_j,\max(\ell_j,v_j)\bigr).$$

The code implements clipping with comparisons so an in-range signed zero is retained consistently across languages. The clip counter includes every out-of-range raw coordinate, even if crossover later discards that coordinate.

1. Build the trial. Use the repaired mutant coordinate if its crossover draw is below $m$ or if it is the forced coordinate; otherwise retain the target coordinate:

$$u_j=\begin{cases}\widehat v_j,&j=j_{\mathrm{rand}}\text{ or }r_j\lt m,\\x_{i,g,j},&\text{otherwise}.\end{cases}$$

1. Evaluate and select. Evaluate the trial once, including when it equals its target. Accept when $f(\mathbf u_i)\le f(\mathbf x_{i,g})$. Accepting a tie permits neutral movement. A rejected trial leaves the target and its stored objective unchanged. Update the archive only on a strictly smaller objective.
2. Commit the generation. Repeat for every target using the same frozen old population, storing winners in a separate next population. After all $N$ trials, replace the population and its score array together, record a snapshot, and repeat.

At least one coordinate is donor-sourced, but the trial can still equal the target numerically. This happens with duplicate vectors, fixed coordinates, zero differences, clipping, or rounding. With ideal independent crossover draws, the expected number of donor-sourced coordinates is $1+(d-1)CR$; it is not necessarily the number of changed values. At $CR=0$ only the forced coordinate uses the donor, and at $CR=1$ every coordinate does. For $d=1$, crossover rate does not affect the chosen value.

flowchart TD
    A["Validate bounds and copy population"] --> B{"Budget covers initialization?"}
    B -->|No| Z["Return no evaluated point"]
    B -->|Yes| C["Evaluate population and archive best"]
    C --> D{"Generation limit reached?"}
    D -->|Yes| Y["Return completed population and archive"]
    D -->|No| E{"Budget covers N trials?"}
    E -->|No| Y
    E -->|Yes| F["Begin separate next population"]
    F --> H["Choose three distinct old donors excluding target"]
    H --> I["Mutate, clip, and cross with forced coordinate"]
    I --> J{"Trial objective no worse than target?"}
    J -->|Yes| K["Store trial and update archive if better"]
    J -->|No| L["Store unchanged target and objective"]
    K --> M{"More targets?"}
    L --> M
    M -->|Yes| H
    M -->|No| N["Commit generation and record snapshot"]
    N --> D

Worked Example with a Full Trace

Use the box $[-5,5]^2$, center $(1,-1)$, population size $N=4$, seed $3$, scale $F=0.5$, crossover threshold $m=500$, generation limit $G=2$, and budget $B=100$. Initialize in this exact order:

The first tied minimum, slot $1$, initializes the archive. The portable generator uses

$$z_{k+1}=48271z_k\bmod2147483647,\qquad z_0=3.$$

Raw states lie in $1,\ldots,2147483646$. To draw an integer below bound $b$, set $R=2147483646$, let $v=z-1$, and accept only if $v\lt R-(R\bmod b)$; return $v\bmod b$ after acceptance and retry otherwise. Under a uniform raw-input model this avoids modulo bias. The seeded generator itself is deterministic and does not supply independent mathematical randomness.

For each target the mapped draw order is three donor ranks, the forced coordinate, and two crossover draws. Even the last donor rank, whose bound is $1$ when $N=4$, consumes a draw. For example, the first raw state is $144813$, and $(144813-1)\bmod3=2$, selecting the third eligible index for target $0$: donor $3$. The donor-rank bounds are $3,2,1$; the forced-coordinate bound is $2$; crossover bounds are $1000$.

The table records every proposal and selection. In a/b/c, $a$ is the base and $b-c$ defines the difference. j is the forced zero-based coordinate. Targets, donor vectors, and incumbent scores always come from the preceding completed generation.

For the first row, $\mathbf x_2-\mathbf x_1=(4,6)$, so the raw mutant is $(5,-2)+0.5(4,6)=(7,1)$. Clip it to $(5,1)$. Coordinate $0$ is forced even though draw $828$ is not below $500$; coordinate $1$ uses the donor because draw $48$ is below $500$. The trial objective is $20$, so it replaces target objective $41$ while the archive remains $13$.

In generation $1$, target $2$ has raw mutant $(-5.5,-1.5)$, clipped to $(-5,-1.5)$. Its first crossover draw $612$ retains the target's first coordinate $3$, so the clipped $-5$ is not used. The final trial $(3,-1.5)$ has objective $4.25$ and becomes the archive. Generation $1$ still uses old slot $0=(-4,3)$ when forming its later trials; the accepted $(5,1)$ is not available as a donor until generation $2$.

Every mapped random choice and each child-end generator state is recorded here:

There are no rejection retries in this sample. Each trial consumes $3+1+2=6$ raw draws, so each generation consumes $24$ and both generations consume $48$. The accepted populations are:

The final result is GEN_LIMIT, best $(0,-2)$ with objective $2$, $12$ evaluations, $8$ trials, $4$ acceptances, $4$ strict target improvements, $4$ accepted point changes, $5$ clipped raw coordinates, and final generator state $1236041904$. The archive is consistent with every evaluated value, yet remains above the analytic minimum $0$.

Implementation in Python, C++, and Java

The three panels are complete programs implementing the same policy. Use Python $3.10+$, C++17, or Java $17+$. For reproducible C++ floating-point decisions, compile with g++ -std=c++17 -O2 -ffp-contract=off main.cpp -o de and avoid fast-math transformations. Save the Java panel as Main.java.

Inputs are an ordered $N\times d$ initial population; bound and center vectors of length $d$; a seed in $1,\ldots,2147483646$; $0\le G\le200$; finite $0\lt F\le2$; integer $0\le m\le1000$; and $0\le B\le10^6$. Require $1\le d\le16$, $4\le N\le32$, finite bounds and center coordinates of magnitude at most $10^6$, ordered bounds, and finite initial points inside their box. Equal lower and upper bounds are allowed. The center may lie outside the box. Python accepts bounded integer or float coordinates but rejects booleans; integer control parameters must be actual integers. Java arrays must be non-null, and Python inputs must be ordinary finite sequences.

Validation occurs before searching. The numeric limits keep raw differential coordinates and the sphere's squared differences finite. Operations are performed in binary64 in a fixed order: subtract, multiply by $F$, add the base; for the objective, subtract the center, square, and accumulate coordinates left to right. The test corpus compared full binary64 traces across languages; this does not make general floating-point arithmetic exact. Tiny differences or products can round to zero, and close real-valued objectives can compare equal.

Result contains status, has_value, best-ever vector and score, final population and scores, completed generations, evaluations, trials, acceptances, strict target improvements, accepted point changes, clipped coordinates, generator statistics, generation snapshots, and per-trial events. An accepted tie need not count as an improvement or as a changed point. A clipped coordinate is counted even when its donor component is discarded by crossover. Snapshots and result arrays are copied, preserving their correspondence with scores and preventing later mutation from rewriting history.

When $B\lt N$, the result has has_value=false, no population, no archive vector, no history, and zero evaluations and draws; the zero score is a dummy field. A request for zero generations still evaluates the initial population when its budget permits. After initialization, GEN_LIMIT takes precedence when both checks would stop. Budget $7$ can evaluate the four initial points but cannot fund another four trials, so it stops at generation $0$ without consuming any randomness. All trials, including unchanged points, are evaluated; unchanged rejected targets retain cached scores.

The supplied objective is deliberately simple and has a known reference answer. For a general deterministic evaluator, replace sphere, check every objective result is finite before selection, and review its cost, constraints, and numerical behavior. Real DE implementations also need explicit policies for failed, noisy, or time-varying evaluations. The sample does not claim to implement an arbitrary objective callback, nonlinear constraints, polishing, or a library's default settings.

The last sample has a fixed one-dimensional box $[2,2]$ and center $0$. Every trial equals $2$, so all four ties are accepted with no changed points. That case correctly returns the box-constrained optimum $4$; it is separate from the suboptimal collapsed-population counterexample below.

from dataclasses import dataclass, field
import math


MOD = 2147483647
RANGE = MOD - 1


class RandomSource:
    def __init__(self, seed):
        self.state, self.draws, self.rejections = seed, 0, 0

    def below(self, bound):
        limit = RANGE - RANGE % bound
        while True:
            self.state = 48271 * self.state % MOD
            self.draws += 1
            value = self.state - 1
            if value < limit:
                return value % bound
            self.rejections += 1


def sphere(point, center):
    total = 0.0
    for x, c in zip(point, center):
        delta = x - c
        term = delta * delta
        total = total + term
    return total


@dataclass
class Event:
    generation: int
    target: int
    donor_ranks: list
    donors: list
    forced: int
    crossover_draws: list
    mutant: list
    repaired: list
    trial: list
    trial_fitness: float
    accepted: bool
    rng_state: int
    rng_draws: int


@dataclass
class Snapshot:
    generation: int
    population: list
    fitness: list
    best: list
    best_fitness: float
    evaluations: int
    rng_draws: int


@dataclass
class Result:
    status: str = "EVAL_LIMIT"
    has_value: bool = False
    best: list = field(default_factory=list)
    best_fitness: float = 0.0
    population: list = field(default_factory=list)
    fitness: list = field(default_factory=list)
    generations: int = 0
    evaluations: int = 0
    trials: int = 0
    accepted: int = 0
    improved: int = 0
    changed: int = 0
    clipped_coordinates: int = 0
    rng_state: int = 0
    rng_draws: int = 0
    rng_rejections: int = 0
    history: list = field(default_factory=list)
    events: list = field(default_factory=list)


def differential_evolution(initial, lower, upper, center, seed, generations,
                           scale, crossover_per_mille, evaluation_budget):
    integers = [seed, generations, crossover_per_mille, evaluation_budget]
    if any(type(x) is not int for x in integers):
        raise ValueError("integer parameters required")
    if not (1 <= seed < MOD and 0 <= generations <= 200 and
            0 <= crossover_per_mille <= 1000 and 0 <= evaluation_budget <= 1000000):
        raise ValueError("parameter out of range")

    def finite_number(x):
        return type(x) is int or (type(x) is float and math.isfinite(x))

    if not finite_number(scale) or not 0.0 < scale <= 2.0:
        raise ValueError("invalid scale")
    dimension, size = len(lower), len(initial)
    if not 1 <= dimension <= 16 or not 4 <= size <= 32:
        raise ValueError("invalid dimension or population size")
    if len(upper) != dimension or len(center) != dimension:
        raise ValueError("vector size mismatch")
    if any(not finite_number(x) or abs(x) > 1000000 for row in
           [lower, upper, center] for x in row):
        raise ValueError("invalid bound or center")
    if any(lower[j] > upper[j] for j in range(dimension)):
        raise ValueError("reversed bounds")
    for row in initial:
        if len(row) != dimension or any(not finite_number(row[j]) or not
                lower[j] <= row[j] <= upper[j] for j in range(dimension)):
            raise ValueError("invalid initial point")
    lo, hi, origin = list(map(float, lower)), list(map(float, upper)), list(map(float, center))
    population = [list(map(float, row)) for row in initial]
    scale = float(scale)
    rng = RandomSource(seed)
    result = Result(rng_state=seed)
    if evaluation_budget < size:
        return result
    fitness = [sphere(row, origin) for row in population]
    elite = min(range(size), key=lambda i: fitness[i])
    result.has_value, result.evaluations = True, size
    result.best, result.best_fitness = population[elite][:], fitness[elite]

    def record():
        result.history.append(Snapshot(result.generations,
            [row[:] for row in population], fitness[:], result.best[:],
            result.best_fitness, result.evaluations, rng.draws))

    record()
    while True:
        if result.generations == generations:
            result.status = "GEN_LIMIT"
            break
        if evaluation_budget - result.evaluations < size:
            result.status = "EVAL_LIMIT"
            break
        next_population, next_fitness = [], []
        for i in range(size):
            excluded, donor_ranks, donors = [i], [], []
            for k in range(3):
                rank = rng.below(size - 1 - k)
                donor_ranks.append(rank)
                index = rank
                for e in sorted(excluded):
                    if index >= e:
                        index += 1
                donors.append(index)
                excluded.append(index)
            a, b, c = donors
            forced = rng.below(dimension)
            draws = [rng.below(1000) for _ in range(dimension)]
            mutant, repaired, trial = [], [], []
            for j in range(dimension):
                difference = population[b][j] - population[c][j]
                scaled = scale * difference
                value = population[a][j] + scaled
                bounded = lo[j] if value < lo[j] else hi[j] if value > hi[j] else value
                result.clipped_coordinates += int(value < lo[j] or value > hi[j])
                mutant.append(value)
                repaired.append(bounded)
                trial.append(bounded if j == forced or draws[j] < crossover_per_mille
                             else population[i][j])
            value = sphere(trial, origin)
            result.evaluations += 1
            result.trials += 1
            accepted = value <= fitness[i]
            if accepted:
                result.accepted += 1
                result.improved += int(value < fitness[i])
                result.changed += int(trial != population[i])
                next_population.append(trial[:])
                next_fitness.append(value)
            else:
                next_population.append(population[i][:])
                next_fitness.append(fitness[i])
            if value < result.best_fitness:
                result.best, result.best_fitness = trial[:], value
            result.events.append(Event(result.generations + 1, i, donor_ranks,
                donors, forced, draws, mutant, repaired, trial, value,
                accepted, rng.state, rng.draws))
        population, fitness = next_population, next_fitness
        result.generations += 1
        record()
    result.population = [row[:] for row in population]
    result.fitness = fitness[:]
    result.rng_state, result.rng_draws = rng.state, rng.draws
    result.rng_rejections = rng.rejections
    return result


def vector(point):
    return "(" + ",".join(f"{x:.6f}" for x in point) + ")"


def show(label, result):
    print(label)
    for row in result.history:
        print(f"g={row.generation} best={vector(row.best)} "
              f"f={row.best_fitness:.6f} eval={row.evaluations} draws={row.rng_draws}")
    best = f"{vector(result.best)}:{result.best_fitness:.6f}" if result.has_value else "NONE"
    print(f"status={result.status} best={best} generations={result.generations} "
          f"eval={result.evaluations} trials={result.trials} accepted={result.accepted} "
          f"improved={result.improved} changed={result.changed} "
          f"clips={result.clipped_coordinates} draws={result.rng_draws} "
          f"rejections={result.rng_rejections} rng={result.rng_state}")


def main():
    initial = [[-4, 3], [-1, -4], [3, 2], [5, -2]]
    for label, budget in [("sphere", 100), ("budget", 7), ("zero-budget", 0)]:
        show(label, differential_evolution(initial, [-5, -5], [5, 5], [1, -1],
            3, 2, 0.5, 500, budget))
    show("fixed", differential_evolution([[2], [2], [2], [2]], [2], [2], [0],
        7, 1, 0.5, 0, 10))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
using namespace std;
using I64 = int64_t;
constexpr I64 MOD = 2147483647, RANGE = MOD - 1;

struct RandomSource {
    I64 state;
    int draws = 0, rejections = 0;
    explicit RandomSource(I64 seed) : state(seed) {}
    int below(int bound) {
        I64 limit = RANGE - RANGE % bound;
        while (true) {
            state = 48271 * state % MOD; ++draws;
            I64 value = state - 1;
            if (value < limit) return int(value % bound);
            ++rejections;
        }
    }
};
double sphere(const vector<double>& point, const vector<double>& center) {
    double total = 0.0;
    for (size_t j = 0; j < point.size(); ++j) {
        double delta = point[j] - center[j];
        double term = delta * delta;
        total = total + term;
    }
    return total;
}
struct Event {
    int generation, target;
    vector<int> donor_ranks, donors;
    int forced;
    vector<int> crossover_draws;
    vector<double> mutant, repaired, trial;
    double trial_fitness;
    bool accepted;
    I64 rng_state;
    int rng_draws;
};
struct Snapshot {
    int generation;
    vector<vector<double>> population;
    vector<double> fitness, best;
    double best_fitness;
    int evaluations, rng_draws;
};
struct Result {
    string status = "EVAL_LIMIT";
    bool has_value = false;
    vector<double> best;
    double best_fitness = 0.0;
    vector<vector<double>> population;
    vector<double> fitness;
    int generations = 0, evaluations = 0, trials = 0, accepted = 0;
    int improved = 0, changed = 0, clipped_coordinates = 0;
    I64 rng_state = 0;
    int rng_draws = 0, rng_rejections = 0;
    vector<Snapshot> history;
    vector<Event> events;
};

Result differential_evolution(const vector<vector<double>>& initial,
        const vector<double>& lower, const vector<double>& upper,
        const vector<double>& center, I64 seed, int generations, double scale,
        int crossover_per_mille, int evaluation_budget) {
    if (seed < 1 || seed >= MOD || generations < 0 || generations > 200 ||
        crossover_per_mille < 0 || crossover_per_mille > 1000 ||
        evaluation_budget < 0 || evaluation_budget > 1000000)
        throw invalid_argument("parameter out of range");
    if (!isfinite(scale) || scale <= 0.0 || scale > 2.0)
        throw invalid_argument("invalid scale");
    int dimension = int(lower.size()), size = int(initial.size());
    if (dimension < 1 || dimension > 16 || size < 4 || size > 32)
        throw invalid_argument("invalid dimension or population size");
    if (upper.size() != lower.size() || center.size() != lower.size())
        throw invalid_argument("vector size mismatch");
    for (int j = 0; j < dimension; ++j) {
        if (!isfinite(lower[j]) || !isfinite(upper[j]) || !isfinite(center[j]) ||
            abs(lower[j]) > 1000000 || abs(upper[j]) > 1000000 || abs(center[j]) > 1000000)
            throw invalid_argument("invalid bound or center");
        if (lower[j] > upper[j]) throw invalid_argument("reversed bounds");
    }
    for (const auto& row : initial) {
        if (row.size() != lower.size()) throw invalid_argument("vector size mismatch");
        for (int j = 0; j < dimension; ++j)
            if (!isfinite(row[j]) || row[j] < lower[j] || row[j] > upper[j])
                throw invalid_argument("invalid initial point");
    }
    auto lo = lower, hi = upper, origin = center;
    auto population = initial;
    RandomSource rng(seed); Result result; result.rng_state = seed;
    if (evaluation_budget < size) return result;
    vector<double> fitness;
    for (const auto& row : population) fitness.push_back(sphere(row, origin));
    int elite = int(min_element(fitness.begin(), fitness.end()) - fitness.begin());
    result.has_value = true; result.evaluations = size;
    result.best = population[elite]; result.best_fitness = fitness[elite];
    auto record = [&]() {
        result.history.push_back(Snapshot{result.generations, population, fitness,
            result.best, result.best_fitness, result.evaluations, rng.draws});
    };
    record();
    while (true) {
        if (result.generations == generations) { result.status = "GEN_LIMIT"; break; }
        if (evaluation_budget - result.evaluations < size) { result.status = "EVAL_LIMIT"; break; }
        vector<vector<double>> next_population;
        vector<double> next_fitness;
        for (int i = 0; i < size; ++i) {
            vector<int> excluded{i}, donor_ranks, donors;
            for (int k = 0; k < 3; ++k) {
                int rank = rng.below(size - 1 - k); donor_ranks.push_back(rank);
                int index = rank; auto ordered = excluded; sort(ordered.begin(), ordered.end());
                for (int e : ordered) if (index >= e) ++index;
                donors.push_back(index); excluded.push_back(index);
            }
            int a = donors[0], b = donors[1], c = donors[2];
            int forced = rng.below(dimension);
            vector<int> draws;
            for (int j = 0; j < dimension; ++j) draws.push_back(rng.below(1000));
            vector<double> mutant, repaired, trial;
            for (int j = 0; j < dimension; ++j) {
                double difference = population[b][j] - population[c][j];
                double scaled = scale * difference;
                double value = population[a][j] + scaled;
                double bounded = value < lo[j] ? lo[j] : value > hi[j] ? hi[j] : value;
                result.clipped_coordinates += value < lo[j] || value > hi[j];
                mutant.push_back(value); repaired.push_back(bounded);
                trial.push_back(j == forced || draws[j] < crossover_per_mille ? bounded : population[i][j]);
            }
            double value = sphere(trial, origin); ++result.evaluations; ++result.trials;
            bool accepted = value <= fitness[i];
            if (accepted) {
                ++result.accepted; result.improved += value < fitness[i];
                result.changed += trial != population[i];
                next_population.push_back(trial); next_fitness.push_back(value);
            } else {
                next_population.push_back(population[i]); next_fitness.push_back(fitness[i]);
            }
            if (value < result.best_fitness) { result.best = trial; result.best_fitness = value; }
            result.events.push_back(Event{result.generations + 1, i, donor_ranks, donors,
                forced, draws, mutant, repaired, trial, value, accepted, rng.state, rng.draws});
        }
        population = move(next_population); fitness = move(next_fitness);
        ++result.generations; record();
    }
    result.population = population; result.fitness = fitness;
    result.rng_state = rng.state; result.rng_draws = rng.draws;
    result.rng_rejections = rng.rejections;
    return result;
}
string vector_text(const vector<double>& point) {
    ostringstream out; out << fixed << setprecision(6) << '(';
    for (size_t i = 0; i < point.size(); ++i) { if (i) out << ','; out << point[i]; }
    out << ')'; return out.str();
}
void show(const string& label, const Result& result) {
    cout << fixed << setprecision(6) << label << '\n';
    for (const auto& row : result.history)
        cout << "g=" << row.generation << " best=" << vector_text(row.best)
             << " f=" << row.best_fitness << " eval=" << row.evaluations << " draws=" << row.rng_draws << '\n';
    cout << "status=" << result.status << " best=";
    if (result.has_value) cout << vector_text(result.best) << ':' << result.best_fitness;
    else cout << "NONE";
    cout << " generations=" << result.generations << " eval=" << result.evaluations
         << " trials=" << result.trials << " accepted=" << result.accepted << " improved=" << result.improved
         << " changed=" << result.changed << " clips=" << result.clipped_coordinates
         << " draws=" << result.rng_draws << " rejections=" << result.rng_rejections
         << " rng=" << result.rng_state << '\n';
}
int main() {
    vector<vector<double>> initial;
    initial.push_back(vector<double>{-4, 3}); initial.push_back(vector<double>{-1, -4});
    initial.push_back(vector<double>{3, 2}); initial.push_back(vector<double>{5, -2});
    vector<double> lo{-5, -5}, hi{5, 5}, center{1, -1};
    show("sphere", differential_evolution(initial, lo, hi, center, 3, 2, .5, 500, 100));
    show("budget", differential_evolution(initial, lo, hi, center, 3, 2, .5, 500, 7));
    show("zero-budget", differential_evolution(initial, lo, hi, center, 3, 2, .5, 500, 0));
    vector<vector<double>> fixed_population(4, vector<double>{2});
    show("fixed", differential_evolution(fixed_population, vector<double>{2}, vector<double>{2},
         vector<double>{0}, 7, 1, .5, 0, 10));
    return 0;
}

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;

public class Main {
    static final long MOD=2147483647L, RANGE=MOD-1;
    static class RandomSource {
        long state; int draws, rejections;
        RandomSource(long seed) { state=seed; }
        int below(int bound) {
            long limit=RANGE-RANGE%bound;
            while (true) {
                state=48271*state%MOD; ++draws; long value=state-1;
                if (value<limit) return (int)(value%bound);
                ++rejections;
            }
        }
    }
    static double sphere(double[] point, double[] center) {
        double total=0.0;
        for (int j=0;j<point.length;++j) {
            double delta=point[j]-center[j]; double term=delta*delta;
            total=total+term;
        }
        return total;
    }
    static double[][] copy(double[][] rows) {
        double[][] result=new double[rows.length][];
        for (int i=0;i<rows.length;++i) result[i]=rows[i].clone();
        return result;
    }
    static class Event {
        int generation, target;
        int[] donor_ranks, donors;
        int forced;
        int[] crossover_draws;
        double[] mutant, repaired, trial;
        double trial_fitness;
        boolean accepted;
        long rng_state; int rng_draws;
        Event(int g,int t,int[] ranks,int[] ds,int forced,int[] draws,
              double[] m,double[] r,double[] tr,double f,boolean ok,long rs,int rd) {
            generation=g;target=t;donor_ranks=ranks;donors=ds;this.forced=forced;
            crossover_draws=draws;mutant=m;repaired=r;trial=tr;trial_fitness=f;
            accepted=ok;rng_state=rs;rng_draws=rd;
        }
    }
    static class Snapshot {
        int generation; double[][] population;
        double[] fitness, best; double best_fitness;
        int evaluations, rng_draws;
        Snapshot(int g,double[][] p,double[] f,double[] b,double bf,int ev,int rd) {
            generation=g;population=copy(p);fitness=f.clone();best=b.clone();
            best_fitness=bf;evaluations=ev;rng_draws=rd;
        }
    }
    static class Result {
        String status="EVAL_LIMIT"; boolean has_value;
        double[] best=new double[0]; double best_fitness;
        double[][] population=new double[0][]; double[] fitness=new double[0];
        int generations,evaluations,trials,accepted,improved,changed,clipped_coordinates;
        long rng_state; int rng_draws,rng_rejections;
        List<Snapshot> history=new ArrayList<>(); List<Event> events=new ArrayList<>();
    }
    static Result differential_evolution(double[][] initial,double[] lower,double[] upper,
            double[] center,long seed,int generations,double scale,int crossover_per_mille,
            int evaluation_budget) {
        if (seed<1 || seed>=MOD || generations<0 || generations>200 ||
            crossover_per_mille<0 || crossover_per_mille>1000 || evaluation_budget<0 || evaluation_budget>1000000)
            throw new IllegalArgumentException("parameter out of range");
        if (!Double.isFinite(scale) || scale<=0.0 || scale>2.0)
            throw new IllegalArgumentException("invalid scale");
        int dimension=lower.length,size=initial.length;
        if (dimension<1 || dimension>16 || size<4 || size>32)
            throw new IllegalArgumentException("invalid dimension or population size");
        if (upper.length!=dimension || center.length!=dimension)
            throw new IllegalArgumentException("vector size mismatch");
        for (int j=0;j<dimension;++j) {
            if (!Double.isFinite(lower[j]) || !Double.isFinite(upper[j]) || !Double.isFinite(center[j]) ||
                Math.abs(lower[j])>1000000 || Math.abs(upper[j])>1000000 || Math.abs(center[j])>1000000)
                throw new IllegalArgumentException("invalid bound or center");
            if (lower[j]>upper[j]) throw new IllegalArgumentException("reversed bounds");
        }
        for (double[] row : initial) {
            if (row.length!=dimension) throw new IllegalArgumentException("vector size mismatch");
            for (int j=0;j<dimension;++j)
                if (!Double.isFinite(row[j]) || row[j]<lower[j] || row[j]>upper[j])
                    throw new IllegalArgumentException("invalid initial point");
        }
        double[] lo=lower.clone(),hi=upper.clone(),origin=center.clone();
        double[][] population=copy(initial);
        RandomSource rng=new RandomSource(seed);Result result=new Result();result.rng_state=seed;
        if (evaluation_budget<size) return result;
        double[] fitness=new double[size];int elite=0;
        for (int i=0;i<size;++i) {
            fitness[i]=sphere(population[i],origin);
            if (fitness[i]<fitness[elite]) elite=i;
        }
        result.has_value=true;result.evaluations=size;
        result.best=population[elite].clone();result.best_fitness=fitness[elite];
        result.history.add(new Snapshot(0,population,fitness,result.best,result.best_fitness,result.evaluations,rng.draws));
        while (true) {
            if (result.generations==generations) { result.status="GEN_LIMIT";break; }
            if (evaluation_budget-result.evaluations<size) { result.status="EVAL_LIMIT";break; }
            double[][] next_population=new double[size][];double[] next_fitness=new double[size];
            for (int i=0;i<size;++i) {
                int[] excluded=new int[4],donor_ranks=new int[3],donors=new int[3];excluded[0]=i;
                for (int k=0;k<3;++k) {
                    int rank=rng.below(size-1-k);donor_ranks[k]=rank;int index=rank;
                    int[] ordered=Arrays.copyOf(excluded,k+1);Arrays.sort(ordered);
                    for (int e : ordered) if (index>=e) ++index;
                    donors[k]=index;excluded[k+1]=index;
                }
                int a=donors[0],b=donors[1],c=donors[2],forced=rng.below(dimension);
                int[] draws=new int[dimension];
                for (int j=0;j<dimension;++j) draws[j]=rng.below(1000);
                double[] mutant=new double[dimension],repaired=new double[dimension],trial=new double[dimension];
                for (int j=0;j<dimension;++j) {
                    double difference=population[b][j]-population[c][j];double scaled=scale*difference;
                    double value=population[a][j]+scaled;
                    double bounded=value<lo[j] ? lo[j] : value>hi[j] ? hi[j] : value;
                    if (value<lo[j] || value>hi[j]) ++result.clipped_coordinates;
                    mutant[j]=value;repaired[j]=bounded;
                    trial[j]=j==forced || draws[j]<crossover_per_mille ? bounded : population[i][j];
                }
                double value=sphere(trial,origin);++result.evaluations;++result.trials;
                boolean accepted=value<=fitness[i];
                if (accepted) {
                    ++result.accepted;if (value<fitness[i]) ++result.improved;
                    boolean changed=false;
                    for (int j=0;j<dimension;++j) if (trial[j]!=population[i][j]) changed=true;
                    if (changed) ++result.changed;
                    next_population[i]=trial.clone();next_fitness[i]=value;
                } else {
                    next_population[i]=population[i].clone();next_fitness[i]=fitness[i];
                }
                if (value<result.best_fitness) { result.best=trial.clone();result.best_fitness=value; }
                result.events.add(new Event(result.generations+1,i,donor_ranks,donors,forced,draws,
                    mutant,repaired,trial,value,accepted,rng.state,rng.draws));
            }
            population=next_population;fitness=next_fitness;++result.generations;
            result.history.add(new Snapshot(result.generations,population,fitness,result.best,
                result.best_fitness,result.evaluations,rng.draws));
        }
        result.population=copy(population);result.fitness=fitness.clone();
        result.rng_state=rng.state;result.rng_draws=rng.draws;result.rng_rejections=rng.rejections;
        return result;
    }
    static String vector(double[] point) {
        StringBuilder s=new StringBuilder("(");
        for (int j=0;j<point.length;++j) {
            if (j>0) s.append(',');s.append(String.format(Locale.ROOT,"%.6f",point[j]));
        }
        return s.append(')').toString();
    }
    static void show(String label,Result result) {
        System.out.println(label);
        for (Snapshot row : result.history)
            System.out.printf(Locale.ROOT,"g=%d best=%s f=%.6f eval=%d draws=%d%n",
                row.generation,vector(row.best),row.best_fitness,row.evaluations,row.rng_draws);
        String best=result.has_value ? vector(result.best)+":"+String.format(Locale.ROOT,"%.6f",result.best_fitness) : "NONE";
        System.out.println("status="+result.status+" best="+best+" generations="+result.generations+
            " eval="+result.evaluations+" trials="+result.trials+" accepted="+result.accepted+
            " improved="+result.improved+" changed="+result.changed+" clips="+result.clipped_coordinates+
            " draws="+result.rng_draws+" rejections="+result.rng_rejections+" rng="+result.rng_state);
    }
    public static void main(String[] args) {
        double[][] initial=new double[4][];
        initial[0]=new double[]{-4,3};initial[1]=new double[]{-1,-4};
        initial[2]=new double[]{3,2};initial[3]=new double[]{5,-2};
        double[] lo={-5,-5},hi={5,5},center={1,-1};
        show("sphere",differential_evolution(initial,lo,hi,center,3,2,.5,500,100));
        show("budget",differential_evolution(initial,lo,hi,center,3,2,.5,500,7));
        show("zero-budget",differential_evolution(initial,lo,hi,center,3,2,.5,500,0));
        double[][] fixed=new double[4][];for (int i=0;i<4;++i) fixed[i]=new double[]{2};
        show("fixed",differential_evolution(fixed,new double[]{2},new double[]{2},new double[]{0},7,1,.5,0,10));
    }
}

Expected standard output, identical in all three languages:

• sphere
• g=0 best=(-1.000000,-4.000000) f=13.000000 eval=4 draws=0
• g=1 best=(3.000000,-1.500000) f=4.250000 eval=8 draws=24
• g=2 best=(0.000000,-2.000000) f=2.000000 eval=12 draws=48
• status=GEN_LIMIT best=(0.000000,-2.000000):2.000000 generations=2 eval=12 trials=8 accepted=4 improved=4 changed=4 clips=5 draws=48 rejections=0 rng=1236041904
• budget
• g=0 best=(-1.000000,-4.000000) f=13.000000 eval=4 draws=0
• status=EVAL_LIMIT best=(-1.000000,-4.000000):13.000000 generations=0 eval=4 trials=0 accepted=0 improved=0 changed=0 clips=0 draws=0 rejections=0 rng=3
• zero-budget
• status=EVAL_LIMIT best=NONE generations=0 eval=0 trials=0 accepted=0 improved=0 changed=0 clips=0 draws=0 rejections=0 rng=3
• fixed
• g=0 best=(2.000000) f=4.000000 eval=4 draws=0
• g=1 best=(2.000000) f=4.000000 eval=8 draws=20
• status=GEN_LIMIT best=(2.000000):4.000000 generations=1 eval=8 trials=4 accepted=4 improved=0 changed=0 clips=0 draws=20 rejections=0 rng=1645572783

Complexity, Correctness, and Invariants

Let $g$ be completed generations. Initialization costs $N$ objective evaluations and every full generation costs another $N$, so for $B\ge N$,

$$g=\min\left(G,\left\lfloor\frac{B-N}{N}\right\rfloor\right),\qquad E=N(g+1),\qquad T=gN,$$

where $T$ is the number of trials. If $B\lt N$, no objective is evaluated. Because this implementation has no early convergence exit, every objective landscape uses the same evaluation count for fixed $N,G,B$. Best, average, and worst solution quality can differ widely even when the work count is identical.

Choosing donors maps ranks past at most three exclusions, so it costs constant work per trial rather than scanning or shuffling all $N$ members. Mutation, clipping, crossover, score evaluation, and trace copying cost $\Theta(d)$ per target for the sphere. If $Q$ raw generator draws include $J$ rejection retries, then

$$Q=T(4+d)+J.$$

With this specific generator and mapped bounds no larger than $1000$, a rejected raw state is followed by an accepted one: a rejected state is $2147483647-k$ for $1\le k\le999$, whose successor $2147483647-48271k$ is outside all these rejection intervals. Thus each mapped draw takes at most two raw draws. Including input validation and copies, best-, average-, and worst-case runtime for this bounded implementation is $\Theta((g+1)Nd)$. Replacing the generator changes that particular retry guarantee; independent uniform raw draws give expected constant mapping cost instead.

For an external evaluator costing $C_f$ per vector, expected search work is $O(NC_f+gN(d+C_f))$, plus input construction and validation. Expensive simulations often dominate variation. No generic polynomial bound follows for the generations needed to approach a global optimum. Runtime to a fixed evaluation budget is a different question from runtime to a requested quality.

Two live populations and their scores use $O(Nd)$ memory. The supplied programs keep full history: snapshots and trial events retain $\Theta((g+1)Nd)$ auxiliary words, including raw mutants, repaired mutants, trials, and crossover draws. Retaining only the current population and archive reduces storage to $O(Nd)$ plus evaluator state. Printing a full high-dimensional trace also has a real cost.

Correctness here means maintaining feasible candidates, valid objective associations, budget accounting, and monotonic progress in the computed deterministic objective. The proof is an induction. Initial points are validated and evaluated. Every donor index exists, differs from the target, and differs from the other donors. Repair places every donor coordinate inside its interval, and crossover chooses between two feasible values, so each trial is feasible. Selection pairs each next vector with the score actually computed for that vector, retaining an old pair on rejection. Finally,

$$f(\mathbf x_{i,g+1})=\min\bigl(f(\mathbf x_{i,g}),f(\mathbf u_{i,g})\bigr)\le f(\mathbf x_{i,g})$$

for every fixed slot $i$. Consequently population minimum, mean, and maximum objective values cannot increase in exact arithmetic over those computed values. Any newly evaluated global record is necessarily accepted, because it is better than every current target; therefore the archive's score equals the population minimum at each completed generation. The archived vector need not still occupy a slot if equal-fitness replacements later change the witnesses. Strict archive updates retain the first evaluated witness of each new record.

Requiring a full generation's evaluation budget before any random draws ensures the result ends with a consistent completed population. None of these facts proves global convergence. For instance, if all population vectors equal a feasible suboptimal point, all differences are zero and this variant remains there forever, regardless of the positive scale or crossover rate.

Variants, Pitfalls, and Optimizations

Donor indices must be distinct and exclude the target. Accidentally selecting $b=c$ makes the differential term zero; allowing repeated indices changes the intended DE/rand/1 policy. Distinct indices can still hold identical vectors, so index correctness alone does not ensure diversity. A donor pool that is too small must be rejected rather than handled by an endless resampling loop.

Force one donor source. Without the forced coordinate, ideal independent crossover takes every coordinate from the target with probability $(1-CR)^d$. At $CR=0$, this happens every time. The forced coordinate prevents that all-target mask, but not a numerically identical trial: a collapsed population or a fixed coordinate can make donor and target values equal.

Use one update policy consistently. Deferred updates read all parents from the old generation and write winners separately. Immediate updates allow later trials to use earlier winners; that is a legitimate different variant but produces different trajectories and affects parallel reproducibility. A best-ever archive may update during this deferred implementation because its values are never used to construct donors. In a best/1 strategy, when the best base becomes available is itself part of the algorithm.

Treat boundary handling as an algorithmic choice. Clipping, reflection, random resampling, and transformed coordinates create different proposal distributions. Clipping can concentrate vectors on box faces and erase variation. It does not enforce a nonlinear feasibility constraint such as a budget or a collision-free route. Specify when repair happens; here it occurs on the mutant before crossover. Do not evaluate an unrepaired point and store its score alongside a repaired point.

Watch population collapse and numerical stagnation. If every vector is $(4,4)$ in the worked box, all donor differences vanish and every trial stays there with objective $34$, although $(1,-1)$ is feasible with objective $0$. Zero spread or unchanged fitness is therefore not an optimality certificate. Very small $F$, cancellation, or floating-point rounding can also suppress steps. Restarts, random immigrants, or adaptation can restore exploration, but they change the method and consume additional evaluations.

Choose the representation and scales deliberately. Bound widths may differ by orders of magnitude. Normalizing coordinates to $[0,1]$ can simplify initialization, tolerances, and variation; return to physical units before evaluating. Axis-wise crossover is not invariant to an arbitrary rotation of the coordinate system, even though differential mutation has a simple vector form. Integer and categorical settings need their own decoding, repair, or an alternative method.

Separate strategy from defaults. best/1 uses the best vector as its base; rand/2 uses two differences and requires enough distinct indices; exponential crossover copies a contiguous run of coordinates rather than making independent coordinate choices. Dithering varies $F$, while adaptive schemes update search parameters from observed outcomes. The SciPy reference documents several strategies and immediate/deferred updating. Its defaults, boundary policy, stopping criteria, random source, and optional polishing differ from this teaching implementation, so matching a strategy name alone does not imply matching traces.

Assess evaluations and uncertainty. Use several independent seeds, compare at equal objective budgets, validate final candidates independently, and include random search and a suitable local method as baselines. If fitness is noisy, accepting one apparently better score and caching it indefinitely can be misleading; use a stated reevaluation or statistical policy. Parallel evaluation fits deferred generations well, but assign randomness and results by target index so completion order does not silently change the run. A smooth convex sphere is a verification case, not evidence that DE is the most efficient solver for that objective.

Practice Problems and Real-World Use

Begin with the shifted sphere and its analytic reference answer. For a box-constrained sphere with arbitrary center, the exact real-valued minimizer clips each center coordinate into its interval; this follows because the squared-coordinate terms minimize independently. Compare that reference point with the algorithm's returned point and objective, allowing for finite-precision evaluation. Then try nonseparable or multimodal functions where the known answer remains available for checking.

In a practical study, document bounds and units, initialization, population size, donor strategy, differential scale, crossover rule, repair policy, update timing, stopping limits, random streams, evaluator failures, and any post-search polishing. Report objective evaluations and variation across seeds as well as the best result. A continuous optimizer can propose useful calibration or engineering parameters, but the evaluator and feasibility checks must reflect the real system.