Why scipy.optimize.differential_evolution does not return last population?

Question

Differential evolution is population based algorithm. However, scipy.optimize.differential_evolution returns result as OptimizeResult which gives only the best elite of the population ignoring the rest of its generation.

How can someone preserve information about the last population?

Answer 1

One way I found to return the last population is to download the source code and work with the class DifferentialEvolutionSolver directly:

solver = DifferentialEvolutionSolver(fun, bounds, arg**)
solver.init_population_random()
solver.__next__() # for each generation

last_pop = solver.population
las_pop_costs = solver.population_energies

Answer 2

It's possible to keep the trial population at each iteration using the ability of the workers keyword to accept a map-like callable that is sent the entire trial population and is expected to return an array with the function values evaluated for the entire trial population:

from scipy.optimize import rosen, differential_evolution
bounds=[(0, 10), (0, 10)]

# pop will retain the trial population at each iteration of the minimisation
pop = []
energies = []

def maplike_fun(func, x):
    # x.shape == (S, N), where S is the size of the population and N
    # is the number of parameters. The rosen function is vectorised,
    # so calling rosen with x.T will return an array of shape (S,)
    # you could also do:
    # v = np.array(list(map(func, x)))
    pop.append(np.copy(x))
    e = func(x.T)
    energies.append(e)
    return e

res = differential_evolution(rosen, bounds, workers=maplike_fun, polish=False, updating='deferred')
# the initial evaluation of the population is not counted in the number
# of iterations, so pop[0] is the initial population.
assert res.nit == len(pop) - 1

To get the actual population at each iteration you then need to iterate through the list of trial populations and successively update pop[0]:

pop_up = [pop[0]]
energies_up = [energies[0]]
for i in range(1, len(pop)):
    new_energies = np.copy(energies_up[-1])
    new_pop = np.copy(pop_up[-1])

    pos = energies[i] < new_energies
    new_energies[pos] = energies[i][pos]
    new_pop[pos] = pop[i][pos]
    
    pop_up.append(new_pop)
    energies_up.append(new_energies)

The actual population evolution is then described by pop_up.

From scipy 1.9 there will also be a vectorized keyword, which will send the entire trial population to the objective function at each iteration.