Title: IMPROVING EVOLUTIONARY STRATEGIES WITH GEN-ERATIVE NEURAL NETWORKS

Abstract: Evolutionary Strategies (ES) are a popular family of black-box zeroth-order optimization algorithms which rely on search distributions to efficiently optimize a large variety of objective functions. This paper investigates the potential benefits of using highly flexible search distributions in ES algorithms, in contrast to standard ones (typically Gaussians). We model such distributions with Generative Neural Networks (GNNs) and introduce a new ES algorithm that leverages their expressiveness to accelerate the stochastic search. Because it acts as a plug-in, our approach allows to augment virtually any standard ES algorithm with flexible search distributions. We demonstrate the empirical advantages of this method on a diversity of objective functions.

1. INTRODUCTION:

We are interested in the global minimization of a black-box objective function, only accessible through a zeroth-order oracle. In many instances of this problem the objective is expensive to evaluate, which excludes brute force methods as a reasonable mean of optimization. Also, as the objective is potentially non-convex and multi-modal, its global optimization cannot be done greedily but requires a careful balance between exploitation and exploration of the optimization landscape (the surface defined by the objective).

The family of algorithms used to tackle such a problem is usually dictated by the cost of one evaluation of the objective function (or equivalently, by the maximum number of function evaluations that are reasonable to make) and by a precision requirement. For instance, Bayesian Optimization (Jones et al., 1998; Shahriari et al., 2016) targets problems of very high evaluation cost, where the global minimum must be approximately discovered after a few hundreds of function evaluations. When aiming for a higher precision and hence having a larger budget (e.g. thousands of function evaluations), a popular algorithm class is the one of Evolutionary Strategies (ES) (Rechenberg, 1978; Schwefel, 1977) , a family of heuristic search procedures.

ES algorithms rely on a search distribution, which role is to propose queries of potentially small value of the objective function. This search distribution is almost always chosen to be a multivariate Gaussian. It is namely the case of the Covariance Matrix Adaptation Evolution Strategies (CMA-ES) (Hansen & Ostermeier, 2001) , a state-of-the-art ES algorithm made popular in the machine learning community by its good results on hyper-parameter tuning (Friedrichs & Igel, 2005; Loshchilov & Hutter, 2016) . It is also the case for Natural Evolution Strategies (NES) (Wierstra et al., 2008) algorithms, which were recently used for direct policy search in Reinforcement Learning (RL) and shown to compete with state-of-the-art MDP-based RL techniques (Salimans et al., 2017) . Occasionally, other distributions have been used; e.g. fat-tails distributions like the Cauchy were shown to outperform the Gaussian for highly multi-modal objectives (Schaul et al., 2011) .

We argue in this paper that in ES algorithms, the choice of a standard parametric search distribution (Gaussian, Cauchy, ..) constitutes a potentially harmful implicit constraint for the stochastic search of a global minimum. To overcome the limitations of classical parametric search distributions, we propose using flexible distributions generated by bijective Generative Neural Networks (GNNs), with computable and differentiable log-probabilities. We discuss why common existing optimization methods in ES algorithms cannot be directly used to train such models and design a tailored algorithm that efficiently train GNNs for an ES objective. We show how this new algorithm can readily incorporate existing ES algorithms that operates on simple search distributions, Algorithm 1: Generic ES procedure input: zeroth-order oracle on f , distribution π 0 , population size λ repeat (Sampling) Sample x 1 , . . . , x λ i.i.d ∼ π t (Evaluation) Evaluate f (x 1 ), . . . , f (x n ). (Update) Update π t to produce x of potentially smaller objective values. until convergence; like the Gaussian. On a variety of objective functions, we show that this extension can significantly accelerate ES algorithms.

We formally introduce the problem and provide background on Evolutionary Strategies in Section 2. We discuss the role of GNNs in generating flexible search distributions in Section 3. We explain why usual algorithms fail to train GNNs for an ES objective and introduce a new algorithm in Section 4. Finally we report experimental results in Section 5.

2. PRELIMINARIES:

In what follows, the real-valued objective function f is defined over a compact X and π will generically denote a probability density function over X . We consider the global optimization of f :

EQUATION

2.1 EVOLUTIONARY STRATEGIES

The generic procedure followed by ES algorithms is presented in Algorithm 1. To make the update step tractable, the search distribution is tied to a family of distributions and parametrized by a realvalued parameter vector θ (e.g. the mean and covariance matrix of a Gaussian), and is referred to as π θ . This update step constitutes the main difference between ES algorithms.

Natural Evolution Strategies One principled way to perform that update is to minimize the expected objective value over samples x drawn from π θ . Indeed, when the search distribution is parametric and tied to a parameter θ, this objective can be differentiated with respect to θ thanks to the log-trick:

EQUATION

This quantity can be approximated from samples -it is known as the score-function or REINFORCE (Williams, 1992) estimator, and provides a direction of update for θ. Unfortunately, naively following a stochastic version of the gradient (2) -a procedure called Plain Gradient Evolutionary Strategies (PGES) -is known to be highly ineffective. PGES main limitation resides in its instability when the search distribution is concentrating, making it unable to precisely locate any local minimum. To improve over the PGES algorithm the authors of Wierstra et al. (2008) proposed to descend J(θ) along its natural gradient (Amari, 1998) . More precisely, they introduce a trust-region optimization scheme to limit the instability of PGES, and minimize a linear approximation of J(θ) under a Kullback-Leibler (KL) divergence constraint:

EQUATION

To avoid solving analytically the trust region problem (3), Wierstra et al. (2008) shows that its solution can be approximated by:

EQUATION

4) is the Fischer Information Matrix (FIM) of π θ . The parameter θ is therefore not updated along the negative gradient of J but rather along F -1 θ ∇ θ J(θ), a quantity known as the natural gradient. The FIM F θ is known analytically when π θ is a multivariate Gaussian and the resulting algorithm, Exponential Natural Evolutionary Strategies (xNES) (Glasmachers et al., 2010) has been shown to reach state-of-the-art performances on a large ES benchmark. CMA-ES Naturally, there exist other strategies to update the search distribution π θ . For instance, CMA-ES relies on a variety of heuristic mechanisms like covariance matrix adaptation and evolution paths, but is only defined when π θ is a multivariate Gaussian. Explaining such mechanisms would be out of the scope of this paper, but the interested reader is referred to the work of Hansen (2016) for a detailed tutorial on CMA-ES.

ES iterations:



2.2. LIMITATIONS OF CLASSICAL SEARCH DISTRIBUTIONS:

ES implicitly balance the need for exploration and exploitation of the optimization landscape. The exploitation phase consists in updating the search distribution, and exploration happens when samples are drawn from the search distribution's tails. The key role of the search distribution is therefore to produce a support adapted to the landscape's structure, so that new points are likely to improve over previous samples.

We argue here that the choice of a given parametric distribution (the multivariate Gaussian distribution being overwhelmingly represented in state-of-the-art ES algorithms) constitutes a potentially harmful implicit constraint for the stochastic search of a global minimum. For instance, a Gaussian distribution is not adapted to navigate a curved valley because of its inability to continuously curve its density. This lack of flexibility will lead it to drastically reduce its entropy, until the curved valley looks locally straight. At this point, the ES algorithm resembles a hill-climber and barely takes advantage of the exploration abilities of the search distribution. An illustration of this phenomenon is presented in Figure 2 on the Rosenbrock function. Another limitation of classical search distribution is their inability to follow multiple hypothesis, that is to explore at the same time different local minima. Even if mixture models can show such flexibility, hyper-parameters like the number of mixtures have optimal values that are impossible to guess a priori.

We want to introduce flexible search distributions to overcome these limitations. Such distributions should, despite their expressiveness, be easily trainable. We should also be concerned when designing them with their role in the exploration/exploitation trade off: a search distribution with too much capacity could over-fit some seemingly good samples, leading to premature convergence. To sum-up, we want to design search-distributions that are:

• more flexible than classical distributions • yet easily trainable • while keeping control over the exploration / exploitation trade-off

In the following section, we carefully investigate the class of Generative Neural Networks (GNNs) to find a parametric class of distributions satisfying such properties.

3. FLEXIBLE SEARCH DISTRIBUTIONS WITH GNNS:

Generative Neural Networks (MacKay, 1995) have been studied in the context of density estimation and shown to be able to model complex and highly multimodal distributions (Srivastava et al., 2017) .

We propose here to leverage their expressiveness for ES, and train them in a principled way thanks to the ES objective:

J(π) = E π [f (x)]

As discussed in Section 2, optimizing J(π) with gradient-based methods is possible through the score-function estimator, which requires to be able to compute and efficiently differentiate the logprobabilities of π.

3.1. GNN BACKGROUND:

The core idea behind a GNN is to map a latent variable z ∈ Z drawn from a known distribution ν ω to an output variable x = g η (z) where g η is the forward-pass of a neural network. The parameter η represents the weights of this neural network while ω describe the degrees of freedom of the latent space distribution ν ω . We denote θ = (ω, η) and π θ (x) the density of the output variable x.

For general neural network architectures, it is impossible to compute π θ (x) for samples x drawn from the GNN. This is namely why their are often trained with adversarial methods (Goodfellow et al., 2014) for sample generation purposes, bypassing the need of computing densities, but at the expense of a good density estimation (mode-dropping). An alternative to adversarial methods was proposed with variational auto-encoders (Kingma & Welling, 2013) however at the cost of learning two neural networks (an encoder and a decoder). A less computationally expensive method consists in restricting the possible architectures to build bijective GNNs, also known as Normalizing Flows (NF) (Rezende & Mohamed, 2015; Papamakarios et al., 2017) , which allows the exact computation of the distribution's density. Indeed, if g η is a bijection from Z to X with inverse h η g -1 η , the change of variable formula provides a way to compute π θ (x):

EQUATION

To have a tractable density one therefore needs to ensure that the determinant of the Jacobian |∂h η (x)/∂x| is easily computable. Several models satisfying these two properties (i.e bijectivity and computable Jacobian) have been proposed for density estimation (Rippel & Adams, 2013; Dinh et al., 2014; 2016) , and proved their expressiveness despite their relatively simple structure.

NFs therefore answer two of our needs when building our new search distribution: flexibility and easiness to train. In this work, we will focus on one NF model: the Non-Linear Independent Component Estimation (Dinh et al., 2014 ) (NICE) model, for its numerical stability and volume preserving properties.

3.2. NICE MODEL:

The authors of NICE proposed to build complex yet invertible transformations through the use of additive coupling layers. An additive coupling layer leaves half of its input unchanged, and adds a non-linear transformation of the first half to the second half. More formally, by noting

v = [v 1 , v 2 ]

the output of a coupling layer and u = [u 1 , u 2 ] its input, one has:

EQUATION

where t is an arbitrarily complex transformation -modelled by a Multi-Layer Perceptron (MLP) with learnable weights and biases. This transformation has unit Jacobian determinant and is easily invertible:

EQUATION

and only requires a feed-forward pass on the MLP t. The choice of the decomposition u = [u 1 , u 2 ] can be arbitrary, and is performed by applying a binary filter to the input. By stacking additive coupling layers, one can create complex distributions, and the inversion of the resulting mapping is independent of the complexity of the neural networks t. The density of the resulting distribution is readily computable thanks to the inverse transform theorem (5).

3.3. VOLUME PRESERVING PROPERTIES:

The transformation induced by NICE is volume preserving (it has a unitary Jacobian determinant). This is quite desirable in a ES context, as the role of concentrating the distribution on a minimum can be left to the latent space distribution ν ω . The role of the additive coupling layers is therefore only to introduce non-linearities in the inverse transform h η so that the distribution is better adapted to the optimization landscape. The fact that this fit is volume-preserving (every subset of the latent space has an image in the data space with the same probability mass) encourages the search distribution to align its tails with regions of small value of the optimization landscape, which is likely to improve the quality of future exploration steps. The NICE model therefore fits perfectly our needs for a flexible search distribution that is easy to train, and that provides enough control on the exploration / exploitation trade-off. Other bijective GNN models like the Real-NVP (Dinh et al., 2016) introduce non-volume preserving transformations, which cannot provide such a control. In practice, we observed that using such transformations for ES led to early concentration and premature convergence.

4. AN EFFICIENT TRAINING ALGORITHM:

We are now equipped with enough tools to use GNNs for ES: an adapted model (NICE) for our search distribution π θ , and an objective to train it with:

EQUATION

Here, θ describes jointly the free parameters of the latent distribution ν ω and η, the weights and biases of the MLPs forming the additive coupling layers.

We start this section by explaining why existing training strategies based on the objective (8) are not sufficient to truly leverage the flexibility of GNNs for ES, before introducing a new algorithm tailored for this task.

4.1. LIMITATIONS OF EXISTING TRAINING STRATEGIES:

We found that the PGES algorithm (naive stochastic gradient descent of ( 8) with the score-function estimator) applied to the NICE distribution suffers from the same limitations as when applied to the Gaussian; it is inable to precisely locate any local minimum. As for the Gaussian, training the NICE distribution for ES requires employing more sophisticated algorithms -such as NES.

However, using the natural gradient for the GNNs distributions is not trivial. First the Fischer Information Matrix F θ is not known analytically and must be estimated via Monte-Carlo sampling, thereby introducing approximation errors. Also, we found that the approximations justifying to follow the descent direction provided by the natural gradient are not adapted to the NICE distribution. Indeed, the assumption behind the NES update (4) is that the loss J(θ) can be (locally) well approximated by the quadratic objective:

J(θ + δθ) = J(θ) + δθ T ∇ θ J(θ) + γ 2 δθ T F θ δθ (9)

where γ is a given non-negative Lagrange multiplier. For NICE, given the highly non-linear nature of π θ this approximation is bound to fail even close to the current parameter θ and will lead to spurious updates. A classical technique (Martens, 2010) to avoid such updates is to artificially increase the curvature of the quadratic term, and is known as damping. Practically, this implies using F θ + βI instead of F θ as the local curvature metric, with β a non-negative damping parameter.

We found that to ensure continuous decrease of J(θ), and because of its highly non-linear nature when using the GNNs, the damping parameter β has to be set to such high values that the modifications of the search distribution are too small to quickly make progress and by no means reaches state-of-the-art performances. We observed that even if the training of the additive coupling layers is performed correctly (i.e the distribution has the correct shape), high damping of the latent space parameters prevents the distribution from quickly concentrating when a minimum is found.

It is unclear how the damping parameter should be adapted to avoid spurious update, while still allowing the distribution to make large step in the latent space and ensure fast concentration when needed. In the following, we present an alternated minimization scheme to bypass the issues raised by natural gradient training for GNN distributions in a ES context.

4.2. ALTERNATING MINIMIZATION:

So far, we used the parameter θ to describe both ω and η (respectively, the free parameters of the latent space distribution ν ω and the degrees of freedom of the non-linear mapping g η ), and the optimization over all these parameters was performed jointly. Separating the roles of ω and η, the initial objective (2) can be rewritten as follows:

EQUATION

Therefore, the initial objective can be rewritten as the expected value of samples drawn from the latent distribution, under the objective f • g η -that is, the representation of the objective function f in the latent space. If ν ω is a standard distribution (i.e efficiently trainable with the natural gradient) and f • g η is a well structured function (i.e one for which ν ω is an efficient search distribution), then the single optimization of ω by classical methods (such as the natural gradient) should avoid the limitations discussed in 2.2. This new representation motivates the design of a new training algorithm that optimizes the parameters ω and η separately.

Alternating Minimization In the following, we will replace the notation π θ with π ω,η to refer to the NICE distribution with parameter θ = (ω, η). We want to optimize ω and η in an alternate fashion, which means performing the following updates at every step of the ES procedure:

EQUATION

EQUATION

This means that at iteration t, samples are drawn from π ωt,ηt and serve to first optimize the latent space distribution parameters ω, and then the additive coupling layers parameters η. For the following iteration, the population is sampled under π ωt+1,ηt+1 .

The update (11a) of the latent space parameters is naturally derived from the new representation (10) of the initial objective. Indeed, ω can be updated via natural gradient ascent of J(ω, η t ) -that is with keeping η = η t fixed. Practically, this therefore reduces to applying a NES algorithm to the latent distribution ν ω on the modified objective function f • g ηt .

Once the latent space parameters updated, the coupling layers parameters should be optimized with respect to:

EQUATION

At this stage, the only available samples are drawn under π ωt,ηt . To estimate, based on these samples, expectations under π ωt+1,ηt one must use importance propensity scores:

EQUATION

The straightforward minimization of this off-line objective is known to lead to degeneracies (Swaminathan & Joachims, 2015, Section 4), and must therefore be regularized. For our application, it is also desirable to make sure that the update η does not undo the progress made in the latent space -in other words, we want to regularize the change in f • g η . To that extent, we adopt a technique proposed in Schulman et al. (2017) and minimize a modification on the initial objective with clipped propensity weights:

EQUATION

clip ε (x) clips the value of x between 1and 1 + . The parameter ε is an hyper-parameter that controls the change in distribution, and the program ( 14) can be efficiently solved via a gradient descent type algorithm, such as Adam (Kingma & Ba, 2014).

To sum up, we propose optimizing the latent distribution and the coupling layers separately. The latent space is optimized by natural gradient descent, and the coupling layers via an off-policy objective with clipped propensity weights. We call this algorithm GNN-ES for Generative Neural Networks Evolutionary Strategies.

Latent space optimization It turns out the GNN-ES can be readily modified to incorporate virtually any existing ES algorithms that operates on the simple distribution ν ω . For instance, if ν ω is set to be a multivariate Gaussian with learnable mean and covariance matrix, the latent space optimization (11a) can be performed by either xNES or CMA-ES. This holds for any standard distribution ν ω and any ES algorithm operating on that distribution. This remark allows us to place GNN-ES in a more general framework and to understand it as a way to improve existing ES algorithm, by providing a principled way to learn complex, non-linear transformations on top of rather standard search distributions (like the Gaussian). In what follows, we will use the GNN prefix in front of existing ES algorithm to describe its augmented version with our algorithm, working as a plug-in.

Pseudo-code for this general algorithm can be found in Appendix B.

4.3. ADDITIONAL TOOLS:

Using historic data ES algorithms typically use small populations of samples to estimate expectations. Such small sample sizes don't allow for enough data exposure for the GNN to build a meaningful transformation g η . To circumvent this problem, we augment the off-line program ( 14) with samples for past generations thanks to the fused importance sampling estimator (Peshkin & Shelton, 2002) . This technique is classical in similar settings like MDP-based reinforcement learning and counterfactual reasoning (Nedelec et al., 2017; Agarwal et al., 2017 ) and proves to be essential for our problem. Formally, for a given horizon T that controls how far we look in the past, this amounts to storing the samples x drawn from π θ t-T +1 , . . . , π θt (as well as their respective scores) in a buffer H T . The objective (13) can then be rewritten as:

EQUATION

This technique allows to increase the data exposure of the GNN by using past samples (and therefore does not require additional function evaluations) and to reduce the variance of the off-line estimator of the original expectation (12) (Nedelec et al., 2017) . To control the change in distribution, the fused propensity weights can then be clipped in a similar fashion than in the program (14).

Mode preserving properties:

To achieve improved exploration, the search distribution should align its tails with the level sets of the objective function. This is not guaranteed when performing the update step ( 14) since the GNN's update could simply move the mean of the search distribution without shaping the tails. One way to encourage the GNN's capacity to be allocated to the tails is to impose a mode-preserving property. If µ denotes the location of a mode of the latent distribution, then the mode of the distribution π θ generated by the NICE model is located in g η (µ) (see Appendix A for the proof). It is therefore easy to build a map f η based on the initial g η that is mode-preserving:

EQUATION

where µ t denotes the mode of the latent distribution ν ω at iteration t. Defined as such, f η preserves the mode of the previous search distribution (since f ηt+1 (µ) = f ηt (µ)), is trivially still a bijection and remains volume preserving. Using the push-forward map f η instead of g η , we explicitly push the flexibility brought by the GNN to impact only the tails of the search distribution. As detailed in an ablation study presented in Appendix F, this additional tool turns out to be essential in order to use GNNs for ES.

5. EXPERIMENTAL RESULTS:

In all that follows, we build the NICE model with three coupling layers. Each coupling layer's nonlinear mapping t is built with a one hidden layer MLP, with 128 neurons and leaky ReLU (Maas et al., 2013) activation functions. This architecture is kept constant in all our experiments.

5.1. VISUALIZATION:

We present here two-dimensional visualizations of the behavior of a GNN distribution trained with GNN-xNES -the latent distribution is therefore Gaussian. Figure 3a displays the density level lines of the resulting search distribution on the Rosenbrock function. Figure 3b displays the density level lines of the latent distribution, as well as the learned representation of the objective in the latent space. The search distribution is able to have curved density isolines, enabling better exploration. In the latent space, the global minimum can be reached without navigating a curved valley. Figures 4a and 4b provide similar visualizations on the Rastrigin function, a highly multimodal but symmetric objective. The GNN lowers the barriers between local minima, making it easier to escape a local minimum to the global minimum.

5.2. SYNTHETIC OBJECTIVES:

Experimental set-up We present experiments on both unimodal and multimodal objectives for xNES and GNN-xNES. We use the official implementation of xNES1 with default hyper-parameters (such as the population size λ), both as a baseline and as an inner optimization method for GNN-xNES. All experiments are run on the COmparing Continous Optimizers (COCO) (Hansen et al., 2016) platform, a popular framework for comparing black-box optimization algorithms. It namely allows to benchmark different algorithms on translated and rotated versions of the same objectives, in order to evaluate multiple configurations with different global minimum positions. We compare xNES and GNN-xNES on functions from the 2018 Black-Box Optimization Benchmark (BBOB) (Hansen et al., 2010) suite. When comparing these two algorithms, we impose that their initial search distributions are close in order to ensure fair comparison. We insist on the fact that the xNES algorithm has the exact same configuration whether it is used by itself or as an inner-optimization algorithm for GNN-xNES. Further experimental details, including additional hyper-parameters value for GNN-xNES are provided in Appendix C.

Unimodal landscapes:

We run the different algorithms on two unimodal landscapes where we expect GNN search distributions to bring a significant improvement over the Gaussian -as discussed in 2.2. These objectives functions are the Rotated Rosenbrock function (a curved valley with high conditioning) and the Bent Cigar (an asymmetric and curved Cigar function). Extensive details on these objective functions can be found in the BBOB documentation (Hansen et al., 2010) . Results on additional unimodal functions can be found in Appendix E.

Performance is measured through Empirical Cumulative Distribution Functions (ECDFs) of the runtime, also known as data profiles (Moré & Wild, 2009) . Such curves report the fraction of problems solved as a function of the number of objective evaluations. For a given precision ∆, a problem is said to be solved if the best function evaluation made so far is smaller than f (x * ) + ∆. We create 200 problems, equally spaced on a log-scale from ∆ = 10 2 to ∆ = 10 -5 and, as in the COCO framework, aggregate them over 15 function instances. Results are presented in Figure 5 for the two benchmark functions and in dimensions d = 2, 5, 10.

Multimodal landscapes:

We now compare the performances of the different algorithms on a collection of three multimodal objectives: the Rastrigin function, the Griewank-Rosenbrock function and the Schwefel function. Extensive details about these objectives can be found in Hansen et al. (2010) . When using ES algorithms to optimize multimodal functions, it is usual to augment them with restart strategies (Hansen, 2016) . When convergence is detected, the search distribution is re-initialized in order to search another part of the landscape, and often the population size is increased. This allows to fairly compared algorithms that converge fast to potentially bad local minima, and algorithms that converges slower to better minima. Their exist a large variety of restart strategies (Loshchilov et al., 2012; Auger & Hansen, 2005) ; as the official implementation of xNES is not equipped with a default one, we trigger a restart whenever the algorithm makes no progress for more than 30 × d iterations.

The standard deviation of the search distribution is set back to 1, and its mean sampled uniformly within the compact X of interest (defined by the COCO framework). At each restart, the population size of the algorithm is multiplied by 2, as in Auger & Hansen (2005) . This restart strategy is used for both xNES and GNN-xNES.

We measure performance as the number of functions evaluations to find an objective value smaller than f (x * ) + 10 -5 within a budget of d × 10 5 function evaluations, averaged over 15 function instances. When an algorithm is not able to discover the global minimum within the given budget, we use the maximum number of evaluations as its performance. For visualization purposes, this measure of performance is divided by d 2 . Results are reported in Figure 6 . On all objectives, and for all dimensions, GNN-xNES discovers (in average) the global minimum faster than xNES. Additional results on others multimodal functions are presented in Appendix E. The goal of this section is to present additional comparison between xNES and GNN-xNES on RL-based objective functions -less synthetic than the previously considered BBOB functions. ES algorithms have recently been used for direct policy search in Reinforcement Learning (RL) and shown to reach performances comparable with state-of-the-art MDP-based techniques (Liu et al., 2019; Salimans et al., 2017) . Direct Policy Search ignores the MDP structure of the RL environment and rather considers it as a black-box. The search for the optimal policy is performed directly in parameter space to maximize the average reward per trajectory:

EQUATION

where p x is the distribution of trajectories induced by the policy (the state-conditional distribution over actions) parametrized by x, and r the rewards generated by the environment. The objective (17) can readily be approximated from samples by simply rolling out M trajectories, and optimized using ES. In our experiments2 , we set M = 10 and optimize deterministic linear policies (as in Rajeswaran et al. (2017) ).

In Figures 7a and 7b we report results of the GNN-xNES algorithm compared to xNES, when run on the Mujoco locomotion tasks Swimmer and InvertedDoublePendulum, both from the OpenAI Gym (Brockman et al., 2016) . Performance is measured by the average reward per trajectory as a function of the number of evaluations of the objective f . Results are averaged over 5 random seeds (ruling the initialization of the environment and the initial distribution over the policy parameters x). In all three environments, GNN-xNES discovers behaviors of high rewards faster than xNES.

6. CONCLUSION:

In this work, we motivate the use of GNNs for improving Evolutionary Strategies by pinpointing the limitations of classical search distributions, commonly used by standard ES algorithms. We propose a new algorithm that leverages the high flexibility of distributions generated by bijective GNNs with an ES objective. We highlight that this algorithm can be seen as a plug-in extension to existing ES algorithms, and therefore can virtually incorporate any of them. Finally, we show its empirical advantages across a diversity of synthetic objective functions, as well as from objectives coming from Reinforcement Learning. Beyond the proposal of this algorithm, we believe that our work highlights the role of expressiveness in exploration for optimization tasks. This idea could be leverage in other settings where exploration is crucial, such a MDP-based policy search methods. An interesting line of future work could focus on optimizing GNN-based conditional distribution for RL tasks -an idea already developed in Ward et al. (2019) ; Mazoure et al. (2019) . Other possible extensions to our work could focus on investigating first-order and mixed oracles, such as in Grathwohl et al. (2017) ; Faury et al. (2018) .

A COMPUTING THE MODE OF THE SEARCH DISTRIBUTION:

We prove here the fact that if µ denotes the location of the mode of the latent distribution ν ω , then g η (µ) is a mode for π ω,η . Indeed, under reasonable smoothness assumptions, one has that y is a mode for π ω,η if and only if:

EQUATION

Since π ω,η (x) = ν ω (h η (x)), this is therefore equivalent to:

EQUATION

In the NICE model, we have that ∂hη(x) ∂x

= 1 for all x hence the matrix ∂hη(x) ∂x x=y is invertible and its kernel is reduced to the null vector. Therefore:

EQUATION

and therefore µ = h η (y) by definition of µ (the only critical point of ν ω ). Hence since h -1 η = g η , we have that y = g η (µ) which concludes the proof.

B ALGORITHM PSEUDO-CODE:

We provide below the pseudo-code for the generic algorithm GNN-A-ES, where A is a generic ES algorithm operating on a parametric distribution ν ω . The additional hyper-parameters are the horizon T as well as the clipping constant ε. The function clip(x, lb, ub) clips the input x between a lower-bound lb and an upper-bound ub. Initialize NICE MLPs weights and biases with η 0 . Let H be a circular buffer of length T × λ while not terminate do (Sampling)

Sample Z = {z 1 , . . . z λ } i.i.d ∼ ν ωt Apply f ηt to Z, obtain X = {x 1 , . . . x λ } i.i.d ∼ π ωt,ηt Evaluate F = {f (x 1 ), . . . , f (x λ )}. Let H ← H + {X, F, π ωt,ηt } (ES update)

//One step ES-based optimization of the latent space Apply A to the latent distribution

ω t+1 ← A (ν ωt , (Z, F )) (GNN update) * //Many-steps gradient based optimization of the GNN η t+1 η x,f ∈H f • clip πω t+1 ,η (x) π∈H π(x) , πω t+1 ,η t (x) π∈H π(x) • (1 -ε), πω t+1 ,η t (x) π∈H π(x) • (1 + ε) end *

The (GNN iteration) step can be performed with virtually any gradient descent solver. In all our experiments, we used Adam (Kingma & Ba, 2014) with learning rate 1e-4 for 500 epochs. 

E ADDITIONAL RESULTS:

We present here some additional results on some unimodal and multimodal synthetic functions. Figure 9 present ECDFs curve obtained from the Attractive Sector function, a highly asymmetrical function around its global minimum. On such a function, GNN-xNES seems to accelerate xNES in small dimensions, however this speed-up disappears in higher dimensions. Figure 10 As detailed in Section 4, one can apply Algorithm 2 as a plug-in to any ES method. So far, we empirically evaluated the benefits of our approach by comparing xNES against its GNN extension (GNN-xNES). We present in Figure 12 additional evaluations obtained by comparing CMA-ES and its GNN extension (denoted GNN-CMA-ES) on the Rosenbrock function in dimension 2,5 and 10. CMA-ES is considered to be the state-of-the-art ES algorithm, and improving its performances is a non-trivial task. On the considered example GNN-CMA-ES improves CMA-ES, highlighting the empirical benefit of our approach for a large class of ES algorithm. One can however observe that the performance boost brought by the GNN extension is milder for GNN-CMA-ES then for GNN-xNES. We suspect that this is due to the use of cumulation via an evolution path in CMA-ES3 , which basically introduces a momentum-like update when optimizing the latent distribution. While using an evolution path makes a lot of sense when optimizing a stationary objective, it can be quite harmful for non-stationary ones. We therefore believe that the cumulation step in CMA-ES (for the latent distribution) and the GNN optimization (making the objective optimized by CMA-ES in the 

F ABLATION STUDY:

We present here an ablation study for two additional tools that we introduced after the alternating optimization view: the mode preserving (16) extension as well as the history augmentation (15). Figure 13 presents ECDFs curves on the Rosenbrock, Rotated Rosenbrock and Bent Cigar functions in 2D, for a version of GNN-xNES that doesn't use history but only the current population. Using history and therefore exposing the GNN to larger datasets improves the procedure. Figure 14 present similar results on a version of GNN-xNES without the mode preserving property ( 16). Again, one can notice that ensuring that the GNN training is mode-preserving is crucial to improve experimental results. 

:

available in the PyBrain library(Schaul et al.,

2010)

We used the rllab library(Duan et al., 2016) for our experiments.

We used the PyCMA library(Hansen et al., 2019) for these experiments.

Figures and Tables:

Figure 2: Example of an undesirable behavior of a Gaussian search distribution. The dashed lines represent density level lines of the search distribution. Because the latter cannot have a curved profile, it is forced to drastically reduce its entropy until it reaches the straight part of the valley.

Figure 3: Rosenbrock

Figure 5: ECDFs curves comparing GNN-xNES and xNES on the Rotated Rosenbrock and Bent Cigar functions, in dimensions d=2,5,10.

Figure 6: Scaling comparison of GNN-xNES and xNES on the Rastrigin, Griewank-Rosenbrock and Schwefel functions, d=2,5,10.

Figure 7: Direct Policy Search experiments

GNN-A-ES (ex: GNN-xNES, GNN-CMA-ES) inputs : objective function f , distribution ν ω and its related ES algorithm A hyper-parameters: clipping constant ε, NICE model architecture, initial parameters ω 0 , initial weights η 0 , horizon T , population size λ (Initialization)

Figure 8: Two-dimensional visualizations. The black dotted lines represent the isolines of the level curves of a NICE search distribution trained with GNN-xNES.

Figure 9: ECDFs curve for the Attractive Sector function, d=2,5,10

Figure 11: Scaling comparison of GNN-xNES and xNES on the Gallagher's Gaussian 101 Peaks and Gallagher's Gaussian 21 Peaks functions, d=2,5,10.

Mean number of restarts needed to discover the global minimum on the Rastrigin function.efficiently fit each optimization landscapes, without having to reduce its entropy like a multivariate normal would.

Table content: <table><tr><td>Algorithm</td><td colspan="3">mean(# restarts), d=2 mean(# restarts), d=5 mean(# restarts), d=10</td></tr><tr><td>xNES</td><td>2.3</td><td>2.7</td><td>3.4</td></tr><tr><td>GNN-xNES</td><td>1.3</td><td>2.5</td><td>2.9</td></tr></table>

presents results on the Rosenbrock function (without random rotations). Again, GNN-xNES accelerates the xNES algorithm. Figure 11 present results on the multimodal functions Gallagher's Gaussian 101 Peaks and Gallagher's Gaussian 21 Peaks. Again, GNN-xNES discovers the global minimum faster (on average) than xNES.In our multimodal experiments, we used simulated restarts as a fair mean of comparing different algorithm (this is common practice in order to fairly compare algorithms that converge fast to potentially bad local minima to algorithms that converge slowly to the global minimum). If the empirical results prove that GNN-xNES accelerate xNES in the discovery of the global minimum, it does not prove that GNN-xNES leverages the flexibility of the GNN to detect the global minimum when xNES misses it. In an attempt to prove that it is indeed the case, we report in Table2the number of restarts needed by both GNN-xNES and xNES to discover the global minimum on the Rastrigin function (averaged over the 15 randomly initialized run). For this instance, GNN-xNES consistently discovers the global minimum with less restarts than xNES.

Table content: <table/>

Peer review (iclr2020):
Summary: In ES the goal is to find a distribution pi_theta(x) such that the expected value of f(x) under this distribution is high. This can be optimized with REINFORCE or with more sophisticated methods based on the natural gradient. The functional form of pi_theta is almost always a Gaussian, but this isn*t sufficiently flexible (e.g. multi-modal) to provide a good optimization algorithm. In response, the authors advocate for using a flexible family of generative neural networks for pi_theta. Using NICE as a generative model is desirable because it maintains volumes. This means that we can adjust volumes in latent space and this directly corresponds to volumes in x space. Doing so is useful to be able to tune how concentrated the search distribution is and to explicitly reason about the mode of the search distribution. Overall, I found that there were a number of technical details that were well motivated, such as how to leverage the *mode preservation* of NICE, how to use importance sampling to be able to use samples from multiple rounds of optimization when updating theta and the fact that any existing ES algorithm can be used to do the optimization in the latent space. *We found that the PGES algorithm (naive stochastic gradient descent of (8) with the score-function estimator) applied to the NICE distribution suffers from the same limitations as when applied to the Gaussian; it is inable to precisely locate any local minimum.* I don*t understand this. Can*t the Gaussian become very concentrated? You write: *ES implicitly balance the need for exploration and exploitation of the optimization landscape. The exploitation phase consists in updating the search distribution, and exploration happens when samples are drawn from the search distribution’s tails.* This is a weak form of exploration, since there is no explicit mechanism that encourages f(x) to be evaluated at regions that it has never been evaluated on before. The search distribution*s tails will have low probability mass, so exploration unlikely. Your proposed method uses a pi(x) that is flexible enough to represent multi-modal distributions. However, how can you ensure that your search procedure actually uses this flexibility? In other words, how is your proposed method any better at exploration that the baseline ES method? It would be great to have a slightly more detailed alg. box in the main text for your proposed method, instead of having it in the appendix. Some details I found confusing, such as whether you perform one step of alternating optimization per call to f(x) or if you perform many steps of alternating optimization. Your *Mode preserving properties* trick is cool. However, I don*t fully understand how it is used. Surely you need to be able to be able to change the mode of the distribution some time? Do you only use the mode preservation trick at certain optimization steps? Again, incorporating this in the alg box would be helpful. In the results, I was disappointed that you required restart strategies. I thought that one of the key advantages of using NICE was that you could capture a multi-modal search distribution. Can you explain? What do you mean by the *global volume of the distribution?* What is the volume of a distribution? I understand what concept you*re trying to convey, but can you be more precise?