diff --git a/es/lib.py b/es/lib.py index 00cea09..2ca2f44 100644 --- a/es/lib.py +++ b/es/lib.py @@ -37,6 +37,20 @@ def default_learning_rate_B_exponential(dimensions): return 0.005 * (9 + 3. * np.log(dimensions)) / (5. * dimensions * np.sqrt(dimensions)) +def default_population_size_elitist(dimensions): + """ + """ + return 10 + + +def default_learning_rates_sigma_elitist(dimensions): + return 1/5 * dimensions**(-3/2.), dimensions**(-3/2.) + + +def default_learning_rate_A_elitist(dimensions): + return 1/4. * dimensions**(-3/2.) + + def utility(fitness): """ Utility function for fitness shaping. diff --git a/es/mo_natural_es.py b/es/mo_natural_es.py new file mode 100644 index 0000000..84d817a --- /dev/null +++ b/es/mo_natural_es.py @@ -0,0 +1,136 @@ +import numpy as np +import pygmo + +from . import lib + + +def is_positive_definite(B): + """Returns true when input is positive-definite, via Cholesky""" + try: + np.linalg.cholesky(B) + return True + except np.linalg.LinAlgError: + return False + + +def optimize(func, cov, + learning_rates_sigma=None, learning_rate_A=None, + population_size=None, max_iter=2000, + fitness_shaping=True, record_history=False, + rng=None): + """ + Multi-objective natural evolution strategies. + See Glasmachers, T., Schaul, T., & Schmidhuber, J. (2010). A + natural evolution strategy for multi-objective optimization. + Parallel Problem Solving from …, 8–11. + """ + if not isinstance(cov, np.ndarray): + raise TypeError('sigma needs to be of type np.ndarray') + + # Dimensionality of parameter space + N = cov.shape[0] + + if learning_rates_sigma is None: + (learning_rate_sigma_minus, + learning_rate_sigma_plus) = lib.default_learning_rates_sigma_elitist(N) + if learning_rate_A is None: + learning_rate_A = lib.default_learning_rate_A_elitist(N) + if population_size is None: + population_size = lib.default_population_size_elitist(N) + if not is_positive_definite(cov): + raise ValueError('covariance matrix needs to be positive semidefinite') + + if rng is None: + rng = np.random.RandomState() + elif isinstance(rng, int): + rng = np.random.RandomState(seed=rng) + + generation = 0 + history_mu = [] + history_cov = [] + history_pop = [] + history_fitness = [] + history_ranks = [] + + # Find Cholesky decomposition of covariance matrix + A = np.linalg.cholesky(cov).T + assert(np.sum(np.abs(np.dot(A.T, A) - cov)) < 1e-12), 'Chochelsky decomposition failed' + + # Decompose A into scalar step and normalized covariance factor A + sigma = (abs(np.linalg.det(A)))**(1. / N) + A /= sigma + + if record_history: + cov = sigma**2 * np.dot(A.T, A) + # history_cov.append(cov.copy()) + # history_pop.append(np.empty((population_size, *np.shape(N)))) + + mu = np.zeros((population_size, N)) + A = np.array([A for i in range(population_size)]) + sigma = np.array([[sigma] for i in range(population_size)]) + # Draw first population + s = rng.normal(0, 1, size=(population_size, N)) + mu = mu + sigma * np.array([np.dot(A[i], s[i]) for i in range(population_size)]) + while True: + s = rng.normal(0, 1, size=(population_size, N)) + mu_prime = mu + sigma * np.array([np.dot(A[i], s[i]) for i in range(population_size)]) + sigma_prime = sigma + A_prime = A + z = np.array([mu, mu_prime]) + + z_flattened = z.reshape(-1, np.shape(mu)[1]) + fitness = np.array([func(zi) for zi in z_flattened]).reshape(2, population_size, -1) + fitness_flattened = fitness.reshape(-1, fitness.shape[-1]) + # Sort fitness values by non_dominant_sorting + # We multiply by -1. because the pygmo routine is written for minimzation + pareto_sets = pygmo.fast_non_dominated_sorting(-1. * fitness_flattened)[0] + hv = pygmo.hypervolume(-1. * fitness_flattened) + ref_point = hv.refpoint(offset=0.1) + + sorted_fitness = np.zeros((0, fitness_flattened.shape[1])) + sorted_parameters = np.zeros((0, *(np.shape(z_flattened)[1:]))) + for pset in pareto_sets: + pfront = fitness_flattened[pset] + par = z_flattened[pset] + hv_i = pygmo.hypervolume(-1. * pfront) + delta_S = hv_i.contributions(ref_point) + sorted_fitness = np.vstack((sorted_fitness, + pfront[np.argsort(delta_S)[::-1]])) + sorted_parameters = np.vstack((sorted_parameters, + par[np.argsort(delta_S)[::-1]])) + # Compute ranks for individuals + ranks = np.array([np.where(sorted_parameters == + z_flattened[i])[0][0] for i in + range(z_flattened.shape[0])]).reshape(2, -1) + # Update step sizes and covariance + for i in range(population_size): + if ranks[1][i] < ranks[0][i]: # Offspring is more successful than parent + sigma[i] *= np.exp(learning_rate_sigma_plus) + sigma_prime[i] *= np.exp(learning_rate_sigma_plus) + A_prime[i] *= (np.exp(learning_rate_A * + (np.outer(s[i], s[i]) - np.eye(z.shape[-1])))) + else: + sigma[i] /= np.exp(learning_rate_sigma_minus) + sigma_prime[i] /= np.exp(learning_rate_sigma_minus) + + # Copy N best individuals into the next generations + mu = np.vstack((mu, mu_prime))[np.argsort(ranks.flatten())[:population_size]] + sigma = np.vstack((sigma, sigma_prime))[np.argsort(ranks.flatten())[:population_size]] + A = np.vstack((A, A_prime))[np.argsort(ranks.flatten())[:population_size]] + + if record_history: + history_mu.append(mu.copy()) + cov = [sigma[i] ** 2 * np.dot(A[i].T, A[i]) for i in range(population_size)] + history_cov.append(cov.copy()) + history_pop.append(z.copy()) + history_fitness.append(fitness.copy()) + history_ranks.append(ranks) + generation += 1 + + # exit if max iterations reached + if generation > max_iter or np.min(sigma) ** 2 < 1e-20: + break + + return {'mu': mu, 'sigma': sigma, 'history_mu': history_mu, + 'history_cov': history_cov, 'history_pop': history_pop, + 'history_fitness': history_fitness, 'history_ranks': history_ranks} diff --git a/test/test_mo_nes.py b/test/test_mo_nes.py new file mode 100644 index 0000000..ac71da2 --- /dev/null +++ b/test/test_mo_nes.py @@ -0,0 +1,47 @@ +import numpy as np +import pylab as pl +from mpl_toolkits.mplot3d import Axes3D +from matplotlib.patches import Ellipse +import sys +from functools import partial +from matplotlib.colors import LogNorm +import seaborn as sns +from matplotlib import gridspec +cp = sns.color_palette() +sys.path.append('../') +from es import separable_natural_es as snes +from es import mo_natural_es as mo_nes +from es import exponential_natural_es as xnes +np.random.seed(123) + +""" +This script tests whether the multi-objective NES finds +multiple solutions to a multi-dimensional fitness function. +""" + +def test_2d_fitness(): + def fitness_add(x): + x = np.array(x).reshape(-1, 2) + return np.sin(0.75*np.sum(x, axis=1))**2 + + + def fitness_subtract(x): + x = np.array(x).reshape(-1, 2) + return -10.*(np.diff(x, axis=1)[:, 0])**2 + + + def fitness(x): + return np.array([fitness_add(x), + fitness_subtract(x)]).T + + + mu = np.array([1., 1.]) + cov = np.array([[0.2, 0.], [0., 0.2]])*10. + + num_iter = 100 + population_size = 10 + + res = mo_nes.optimize(fitness, cov, max_iter=num_iter, + record_history=True, population_size=population_size, + rng=12354) + assert(np.allclose(res['mu'], np.pi/3., atol=0.1))