OptimLib: Differential Evolution

Differential Evolution (DE) is a stochastic genetic search algorithm for global optimization of potentially ill-behaved nonlinear functions.

Definition and Syntax

bool de(arma::vec& init_out_vals, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)> opt_objfn, void* opt_data);
bool de(arma::vec& init_out_vals, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)> opt_objfn, void* opt_data, algo_settings_t& settings);

Function arguments:

init_out_vals a column vector of initial values; will be replaced by the solution values.
opt_objfn the function to be minimized, taking three arguments:
- vals_inp a vector of inputs;
- grad_out an empty vector, as DE does not require the gradient to be known/exist; and
- opt_data additional parameters passed to the function.
opt_data additional parameters passed to the function.
settings parameters controlling the optimization routine; see below.

Optimization control parameters:

bool vals_bound whether the search space is bounded. If true, then

arma::vec lower_bounds this defines the lower bounds.
arma::vec upper_bounds this defines the upper bounds.

int de_n_pop population size of each generation.
int de_n_gen number of vectors to generate.
int de_check_freq number of generations between convergence checks.

int de_mutation_method which method of mutation to apply:

de_mutation_method = 1 applies the 'rand' policy.
de_mutation_method = 2 applies the 'best' policy.

double de_par_F the mutation parameter $F$ in the details section below.
double de_par_CR the crossover parameter $CR$ in the details section below.

arma::vec de_initial_lb lower bounds on the initial population; defaults to init_out_vals $- \ 0.5$.
arma::vec de_initial_ub upper bounds on the initial population; defaults to init_out_vals $+ \ 0.5$.

Details

The DE method in OptimLib comes in two varieties: the simple version, outlined first, and the more advanced method of Zamuda and Brest (2015).

Basic Differential Evolution (DE)

Let $X^{(i)}$ denote a $N_p \times d$-dimensional array of values at stage $i$ of the algorithm. The basic DE algorithm is comprised of three steps.

The Mutation Step. For unique random indices $a,b,c$:

if de_mutation_method = 1, use the 'rand' method:
if de_mutation_method = 2, use the 'best' method:

The Crossover Step. For a random integer $r_k \in \{1, \ldots, d\}$,
The Update Step.

The algorithm stops when the relative improvement in the objective function is less than err_tol between de_check_freq number of generations, or when the total number of 'generations' exceeds n_gen.

Differential Evolution with Population Reduction and Multiple Mutation Strategies (DE-PRMM)

Let $X^{(i)}$ denote a $N_p \times d$-dimensional array of values at stage $i$ of the algorithm.

Let $M$ denote the maximum number of function evaluations (of which there are $N_p$ per generation).
Each generation is split into two sub-populations of size $N_{p}^{(m)}$ and $N_{p}^{(b)} = N_p - N_{p}^{(m)}$, labelled the 'main' population and the 'best' population, respectively.
The population is reduced $p_{\text{max}}-1$ number of times.

The steps are as follows.

Population Reduction Step.
If $$i = \dfrac{M}{p_{\text{max}} N_p}$$ reduce the population by half ($N_p /2$) and double the number of generations ($2 \times N_g$).
Values are chosen based on a simple pairwise comparison: $$X_{i,\text{new}} (j,:) = \begin{cases} X^{(i)} (j,:) &\text{ if } f(X^{(i)} (j,:)) < f(X^{(i)} (j+N_p/2,:)) \\ X^{(i)} (j+N_p/2,:) & \text{ else} \end{cases}$$ Reset $i = 1$.
The Mutation Step.
Sample $F_i$ and $CR_i$ values according to: $$F_i = \begin{cases} F_l + (F_u - F_l) \times u_2 & \text{ if } u_1 < \tau_F \\ F_{i-1} & \text{ else} \end{cases}$$ $$CR_i = \begin{cases} u_4 & \text{ if } u_3 < \tau_{CR} \\ CR_{i-1} & \text{ else} \end{cases}$$ If $i \leq N_p - N_p^{(b)}$:
- then for unique random indices $a,b,c$:
Else, if $N_p - N_p^{(b)} < i \leq N_p$:
- then for unique random indices $a,b$:
The Crossover Step.
For a random integer $r_k \in \{1, \ldots, d\}$,
The Update Step.
Setting the 'best' vector.
Let $X_{i+1,\text{best}}^{(m)} := \arg \min \{ f(X^{(i+1)}(1,:)), \ldots, f(X^{(i+1)}(N_p^{(m)},:)) \}$ and $X_{i+1,\text{best}}^{(b)} := \arg \min \{ f(X^{(i+1)}(N_p^{(m)}+1,:)), \ldots, f(X^{(i+1)}(N_p,:)) \}$.
- If $i = N_p^{(m)}$ and $f(X_{i+1,\text{best}}^{(m)}) < f(X_{i,\text{best}}^{(b)})$, set
- If $i = N_p$ and $f(X_{i+1,\text{best}}^{(b)}) < f(X_{i+1,\text{best}}^{(m)})$, and

Examples

To illustrate the effectiveness of DE, we will tackle two well-known performance tests from the numerical optimization literature:

The Rastrigin function:
Ackley's function:

Both functions are 'bumpy' and contain many local minima. Plots of both functions are given below.

Rastrigin

Ackley

Both minimization problems possess a unique global minimum: the zero vector. Code implementing these examples using OptimLib is given below.

#include "optim.hpp"

struct rastrigin_fn_data {
    double A;
};

double rastrigin_fn(const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)
{
    const int n = vals_inp.n_elem;

    rastrigin_fn_data* objfn_data = reinterpret_cast<rastrigin_fn_data*>(opt_data);
    const double A = objfn_data->A;

    double obj_val = A*n + arma::accu( arma::pow(vals_inp,2) - A*arma::cos(2*arma::datum::pi*vals_inp) );
    //
    return obj_val;
}

double ackley_fn(const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)
{
    const double x = vals_inp(0);
    const double y = vals_inp(1);
    const double pi = arma::datum::pi;

    double obj_val = -20*std::exp( -0.2*std::sqrt(0.5*(x*x + y*y)) ) - std::exp( 0.5*(std::cos(2*pi*x) + std::cos(2*pi*y)) ) + std::exp(1) + 20;
    //
    return obj_val;
}

int main()
{
    //
    // Rastrigin function

    rastrigin_fn_data test_data;
    test_data.A = 10;

    arma::vec x = arma::ones(2,1) + 1.0; // initial values: (2,2)

    bool success = optim::de(x,rastrigin_fn,&test_data);

    if (success) {
        std::cout << "de: Rastrigin test completed successfully." << std::endl;
    } else {
        std::cout << "de: Rastrigin test completed unsuccessfully." << std::endl;
    }

    arma::cout << "de: solution to Rastrigin test:\n" << x << arma::endl;

    //
    // Ackley function

    arma::vec x_2 = arma::ones(2,1) + 1.0; // initial values: (2,2)

    bool success_2 = optim::de(x_2,ackley_fn,nullptr);

    if (success_2) {
        std::cout << "de: Ackley test completed successfully." << std::endl;
    } else {
        std::cout << "de: Ackley test completed unsuccessfully." << std::endl;
    }

    arma::cout << "de: solution to Ackley test:\n" << x_2 << arma::endl;

    return 0;
}