1 Introduction
Blackbox optimization (also called zeroorder optimization) methods have received a great deal of attention in recent years due to their broad applicability to real world problems [18, 15, 10, 7, 9, 11]. When the structure of the objective function is unknown or too complex to model directly, or when gradient information is unavailable or unreliable, such methods are often seen as the last resort because all they require is that the objective function can be evaluated at specific points.
In continuous blackbox optimization problems, the stateoftheart algorithms, such as xNES [4] and CMAES [8], are all based the same principle [4, 1]: a Gaussian search distribution is repeatedly updated based on the objective function values of sampled points. Usually the full covariance matrix of the distribution is updated, allowing the algorithm to adapt the size and shape of the Gaussian to the local characteristics of the objective function. Full parameterization also provides invariance under affine transformations of the coordinate system, so that illshaped, highly nonseparable problems can be tackled. However, this generality comes with a price. The number of parameters scales quadratically in the number of dimensions, and the computational cost per sample is at least quadratic in the number of dimensions [13], and sometimes cubic [4, 16].
This cost is often justified since evaluating the objective function can dominate the computation, and thus the main focus is on the improvement of sampling efficiency. However, there are many problems, such as optimizing weights in neural networks where the dimensionality of the parameter space can be very large (e.g. many thousands of weights), the quadratic cost of updating the search distribution can become the computational bottleneck. One possible remedy is to restrict the covariance matrix to be diagonal
[14], which reduces the computation per function evaluation to , linear in the number of dimensions. Unfortunately, this “diagonal” approach performs poorly when the problem is nonseparable because the search distribution cannot follow directions that are not parallel to the current coordinate axes.In this paper, we propose a new variant of the natural evolution strategy family [17], termed Rank One NES (R1NES). This algorithm stays within the general NES framework in that the search distribution is adjusted according to the natural gradient [2], but it uses a novel parameterization of the covariance matrix,
where and are the parameters to be adjusted. This parameterization allows the predominant eigen direction, , of to be aligned in any direction, enabling the algorithm to tackle highly nonseparable problems while maintaining only parameters. We show through rigorous derivation that the natural gradient can also be effectively computed in per sample. R1NES scales well to high dimensions, and dramatically outperforms diagonal covariance matrix algorithms on nonseparable objective functions. As an example, R1NES reliably solves the the nonconvex Rosenbrock function up to dimensions.
2 The NES framework
Natural evolution strategies (NES) are a class of evolutionary algorithms for realvalued optimization that maintain a search distribution, and adapt the distribution parameters by following the
natural gradient of the expected function value. The success of such algorithms is largely attributed to the use of natural gradient, which has the advantage of always pointing in the direction of the steepest ascent, even if the parameter space is not Euclidean. Moreover, compared to regular gradient, natural gradient reduces the weights of gradient components with higher uncertainty, therefore making the update more reliable. As a consequence, NES algorithms can effectively cope with objective functions with illshaped landscapes, especially preventing premature convergence on plateaus and avoiding overaggressive steps on ridges [16].The general framework of NES is given as follows: At each time step, the algorithm samples new samples , with being the search distribution parameterized by . Let be the objective function to maximize. The expected function value under the search distribution is
Using the loglikelihood trick, the gradient w.r.t. the parameters can be written as
from which we obtain the MonteCarlo estimate
of the search gradient. The key step of NES then consists of replacing this gradient by the natural gradient
where
is the Fisher information matrix (See Fig.1 for an illustration). This leads to a straightforward scheme of natural gradient ascent for iteratively updating the parameters
(1) 
The sequence of 1) sampling an offspring population, 2) computing the corresponding Monte Carlo estimate of the gradient, 3) transforming it into the natural gradient, and 4) updating the search distribution, constitutes one iteration of NES.
The most difficult step in NES is the computation of the Fisher information matrix with respect to the parameterization. For full Gaussian distribution, the Fisher can be derived analytically
[4, 16]. However, for arbitrary parameterization of , the Fisher matrix can be highly nontrivial.3 Natural gradient of the rankone covariance matrix approximation
In this paper, we consider a special formulation of the covariance matrix
with parameter set
. The special part of the parameterization is the vector
, which corresponds to the predominant direction of . This allows the search distribution to be aligned in any direction by adjusting , enabling the algorithm to follow valleys not aligned with the current coordinate axes, which is essential for solving nonseparable problems.Since should always be positive, following the same procedure in [4], we parameterize , so that can be adjusted freely using gradient descent without worrying about becoming negative. The parameter set is adjusted to accordingly.
From the derivation of [16], the natural gradient on the sample mean is given by
(2) 
In the subsequent discussion we always assume for simplicity. It is straightforward to sample from ^{1}^{1}1For succinctness, we always assume the mean of the search distribution is . This can be achieved easily by shifting the coordinates. by lettting , , then
The inverse of can also be computed easily as
where . Using the relation , the determinant of is
Knowing and allows the loglikelihood to be written explicitly as
The regular gradient with respect to and can then be computed as:
(3)  
(4) 
Replacing with , then the Fisher can be computed by marginalizing out i.i.d. standard Gaussian variables and , namely,
Since elements in are essentially polynomials of and , their expectations can be computed analytically^{2}^{2}2The derivation is tedious, thus omitted here. All derivations are numerically verified using MonteCarlo simulation., which gives the exact Fisher information matrix
with
Let , then
The inverse of is thus given by
We apply the formula for block matrix inverse in [12]
where , and are the Schur complements. Let be partitioned as above, then
and the Shur complements are
and
whose inverse is given by
Combining the results gives the analytical solution of the inverse Fisher:
Multiplying with the regular gradient in Eq.3 and Eq.4 gives the natural gradient for and :
(5) 
and
(6) 
Note that computing both and requires only the inner products and , therefore can be done storage and time.
3.1 Reparameterization
The natural gradient above is obtained with respect to . However, direct gradient update on has an unpleasant property when is in the opposite direction of , which is illustrated in Fig. 2(a). In this case, the gradient tends to shrink . However, if is large, adding the gradient will flip the direction of , and the length of might even grow. This causes numerical problems, especially when is small. A remedy is to separate the length and direction of , namely, reparameterize , where and is the length of . Then the gradient update on will never flip , and thus avoid the problem.
Note that for small change , the update on and can be obtained from
and
The natural gradient on and is given by letting , thanks to the invariance property:
(7)  
(8) 
Note that computing and involves only inner products between vectors, which can also be done linearly in the number of dimensions.
Using the parameterization introduces another problem. When is small, tends to be large, and thus directly updating causes to grow exponentially, resulting in numerical instability, as shown in Fig. 2(b). In this case, the additive update on , rather than the update on is more stable. In our implementation, the additive update on is used if , otherwise the update is on . This solution proved to be numerially stable in all our tests. Algorithm 1 shows the complete R1NES algorithm in pseudocode.
4 Experiments
The R1NES algorithm was evaluated on the twelve noisefree unimodal functions [6] in the ‘BlackBox Optimization Benchmarking’ collection (BBOB) from the 2010 GECCO Workshop for RealParameter Optimization. In order to make the results comparable those of other methods, the setup in [5] was used, which transforms the pure benchmark functions to make the parameters nonseparable (for some) and avoid trivial optima at the origin.
R1NES was compared to xNES [3], SNES [14] on each benchmark with problem dimensions , (20 runs for each setup), except for xNES, which was only run up . Note that xNES serves as a proper baseline since it is stateoftheart, achieving performance on par with the popular CMAES. The reference machine is an Intel Core i7 processor with 1.6GHz and 4GB of RAM.
Fig. 3 shows the results for the eight benchmarks on which R1NES performs at least as well as the other methods, and often much better. For dimensionality under , R1NES is comparable to xNES in terms of the number of fitness evaluations, indicating that the rankone parameterization of the search distribution effectively captures the local curvature of the fitness function (see Fig.5 for example). However, the time required to compute the update for the two algorithms differs drastically, as depicted in Fig. 4. For example, a typical run of xNES in dimensions takes hours (hence the truncated xNES curves in all graphs), compared to minutes for R1NES. As a result, R1NES can solve these problems up to dimensions in acceptable time. In particular, the result on the dimensional Rosenbrock function is, to our knowledge, the best to date. We estimate that optimizing the 512dimensional sphere function with xNES (or any other full parameterization method, e.g. CMAES) would take over a year in computation time on the same reference hardware. It is also worth pointing out that sNES, though sharing similar, low complexity per function evaluation, can only solve separable problems (Sphere, Linear, AttractiveSector, and Ellipsoid).
Fig. 6 shows four cases (Ellipsoid, StepEllipsoid, RotatedEllipsoid, and Tablet) for which R1NES is not suited, highlighting a limitation of the algorithm. Three of the four functions are from the Ellipsoid family, where the fitness functions are variants of the type
The eigenvalues of the Hessian span several orders of magnitude, and the parameterization with a single predominant direction is not enough to approximate the Hessian, resulting in poor performance. The other function where R1NES fails is the Tablet function where all but a one eigendirection has a large eigenvalue. Since the parameterization of R1NES only allows a single direction to have a large eigenvalue, the shape of the Hessian cannot be effectively approximated.
5 Conclusion and future work
We presented a new blackbox optimization algorithm R1NES that employs a novel parameterization of the search distribution covariance matrix which allows a predominant search direction to be adjusted using the natural gradient with complexity linear in the dimensionality. The algorithm shows excellent performance in a number of highdimensional nonseparable problems that, to date, have not been solved with other parameterizations of similar complexity.
Future work will concentrate on overcoming the limitations of the algorithm (shown in Fig 6). In particular, we intend to extend the algorithm to a) incorporate multiple search directions, and b) enable each search direction to shrink as well as grow.
Acknowledgement
This research was funded in part by Swiss National Science Foundation grants 200020122124, 200020125038, and EU IMCLeVeR project (#231722).
References
 [1] Y. Akimoto, Y. Nagata, I. Ono, and S. Kobayashi. Bidirectional Relation between CMA Evolution Strategies and Natural Evolution Strategies. In Parallel Problem Solving from Nature (PPSN), 2010.
 [2] S. Amari. Natural gradient works efficiently in learning. Neural Computation, 10(2):251–276, 1998.

[3]
T. Glasmachers, T. Schaul, Y. Sun, D. Wierstra, and J. Schmidhuber.
Exponential Natural Evolution Strategies.
In
Genetic and Evolutionary Computation Conference (GECCO)
, Portland, OR, 2010.  [4] T. Glasmaches, T. Schaul, Y. Sun, and J. Schmidhuber. Exponential natural evolution strategies. In GECCO’10, 2010.
 [5] N. Hansen and A. Auger. Realparameter blackbox optimization benchmarking 2010: Experimental setup, 2010.
 [6] N. Hansen and S. Finck. Realparameter blackbox optimization benchmarking 2010: Noiseless functions definitions, 2010.
 [7] N. Hansen, A. S. P. Niederberger, L. Guzzella, and P. Koumoutsakos. A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion. Trans. Evol. Comp, 13:180–197, 2009.
 [8] N. Hansen and A. Ostermeier. Completely derandomized selfadaptation in evolution strategies. Evolutionary Computation, 9(2):159–195, 2001.
 [9] M. Hasenjäger, B. Sendhoff, T. Sonoda, and T. Arima. Three dimensional evolutionary aerodynamic design optimization with CMAES. In Proceedings of the 2005 conference on Genetic and evolutionary computation, GECCO ’05, pages 2173–2180, New York, NY, USA, 2005. ACM.
 [10] M. Jebalia, A. Auger, M. Schoenauer, F. James, and M. Postel. Identification of the isotherm function in chromatography using cmaes. In IEEE Congress on Evolutionary Computation, pages 4289–4296, 2007.
 [11] J. Klockgether and H. P. Schwefel. Twophase nozzle and hollow core jet experiments. In Proc. 11th Symp. Engineering Aspects of Magnetohydrodynamics, pages 141–148, 1970.
 [12] K. B. Petersen and M. S. Pedersen. The matrix cookbook, 2008.
 [13] R. Ros and N. Hansen. A Simple Modification in CMAES Achieving Linear Time and Space Complexity. In R. et al., editor, Parallel Problem Solving from Nature, PPSN X, pages 296–305. Springer, 2008.
 [14] T. Schaul, T. Glasmachers, and J. Schmidhuber. High Dimensions and Heavy Tails for Natural Evolution Strategies. In To appear in: Genetic and Evolutionary Computation Conference (GECCO), 2011.
 [15] O. M. Shir and T. Bäck. The second harmonic generation casestudy as a gateway for es to quantum control problems. In Proceedings of the 9th annual conference on Genetic and evolutionary computation, GECCO ’07, pages 713–721, New York, NY, USA, 2007. ACM.
 [16] Y. Sun, D. Wierstra, T. Schaul, and J. Schmidhuber. Stochastic search using the natural gradient. In ICML’09, 2009.
 [17] D. Wierstra, T. Schaul, J. Peters, and J. Schmidhuber. Natural Evolution Strategies. In Proceedings of the Congress on Evolutionary Computation (CEC08), Hongkong. IEEE Press, 2008.
 [18] S. Winter, B. Brendel, and C. Igel. Registration of bone structures in 3d ultrasound and ct data: Comparison of different optimization strategies. International Congress Series, 1281:242 – 247, 2005.
Comments
There are no comments yet.