ARY STRATEGIES WITH PAST DESCENT DIRECTIONS
Year: 2019

Abstract:
Evolutionary Strategies (ES) are known to be an effective black-box optimization technique for deep neural networks when the true gradients cannot be computed, such as in Reinforcement Learning. We continue a recent line of research that uses surrogate gradients to improve the gradient estimation of ES. We propose a novel method to optimally incorporate surrogate gradient information. Our approach, unlike previous work, needs no information about the quality of the surrogate gradients and is always guaranteed to find a descent direction that is better than the surrogate gradient. This allows to iteratively use the previous gradient estimate as surrogate gradient for the current search point. We theoretically prove that this yields fast convergence to the true gradient for linear functions and show under simplifying assumptions that it significantly improves gradient estimates for general functions. Finally, we evaluate our approach empirically on MNIST and reinforcement learning tasks and show that it considerably improves the gradient estimation of ES at no extra computational cost.

heading: 1 INTRODUCTION
text: Evolutionary Strategies (ES) (1; 2; 3) are a black-box optimization technique, that estimate the gradient of some objective function with respect to the parameters by evaluating parameter perturbations in random directions. The benefits of using ES in Reinforcement Learning (RL) were exhibited in (4). ES approaches are highly parallelizable and account for robust learning, while having decent data-efficiency. Moreover, black-box optimization techniques like ES do not require propagation of gradients, are tolerant to long time horizons, and do not suffer from sparse reward distributions (4). This lead to a successful application of ES in variety of different RL settings (5; 6; 7; 8). Applications of ES outside RL include for example meta learning (9).
In many scenarios, the true gradient is impossible to compute, however surrogate gradients are available. Here, we use the term surrogate gradients for directions that are correlated but usually not equal to the true gradient, e.g. they might be biased or unbiased approximations of the gradient. Such scenarios include models with discrete stochastic variables (10), learned models in RL like Q-learning (11), truncated backpropagation through time (12) and feedback alignment (13), see (14) for a detailed exhibition. If surrogate gradients are available, it is beneficial to preferentially sample parameter perturbations from the subspace defined by these directions (14). The proposed algorithm (14) requires knowing in advance the quality of the surrogate gradient, does not always provide a descent direction that is better than the surrogate gradient, and it remains open how to obtain such surrogate gradients in general settings.
In deep learning in general, experimental evidence has established that higher order derivatives are usually "well behaved", in which case gradients of consecutive parameter updates correlate and applying momentum speeds up convergence (15; 16; 17). These observations suggest that past update directions are promising candidates for surrogate gradients.
In this work, we extend the line of research of (14). Our contribution is threefold:
• First, we show theoretically how to optimally combine the surrogate gradient directions with random search directions. More precisely, our approach computes the direction of the subspace spanned by the evaluated search directions that is most aligned with the true gradient. Our gradient estimator does not need to know the quality of the surrogate gradients
and always provides a descent direction that is more aligned with the true gradient than the surrogate gradient. • Second, above properties of our gradient estimator allow us to iteratively use the last update
direction as a surrogate gradient for our gradient estimator. Repeatedly using the last update direction as a surrogate gradient will aggregate information about the gradient over time and results in improved gradient estimates. In order to demonstrate how the gradient estimate improves over time, we prove fast convergence to the true gradient for linear functions and show, that under simplifying assumptions, it offers an improvement over ES that depends on the Hessian for general functions. • Third, we validate experimentally that these results transfer to practice, that is, the proposed
approach computes more accurate gradients than standard ES. We observe that our algorithm considerably improves gradient estimation on the MNIST task compared to standard ES and that it improves convergence speed and performance on the tested Roboschool reinforcement learning environments.

heading: 2 RELATED WORK
text: Evolutionary strategies (1; 2; 3) are black box optimization techniques that approximate the gradient by sampling finite differences in random directions in parameter space. Promising potential of ES for the optimization of neural networks used for RL was demonstrated in (4). They showed that ES gives rise to efficient training despite the noisy gradient estimates that are generated from a much smaller number of samples than the dimensionality of parameter space. This placed ES on a prominent spot in the RL tool kit (5; 6; 7; 8).
The history of descent directions was previously used to adapt the search distribution in covariance matrix adaptation ES (CMA-ES) (18). CMA-ES constructs a second-order model of the underlying objective function and samples search directions and adapts step size according to it. However, maintaining the full covariance matrix makes the algorithm quadratic in the number of parameters, and thus impractical for high-dimensional spaces. Linear time approximations of CMA-ES like diagonal approximations of the covariance matrix (19) often do not work well, in the sense that their gradient estimates do not converge to the true gradient even if the true gradient does not change over time. Our approach differs as we simply improve the gradient estimation and then feed the gradient estimate to a first-order optimization algorithm.
Our work is inspired by the line of research of (14), where surrogate gradient directions are used to improve gradient estimations by ’elongating’ the search space along these directions. That approach has two shortcomings. First, the bias of the surrogate gradients needs to be known to adapt the covariance matrix. Second, once the bias of the surrogate gradient is too small, the algorithm will not find a better descent direction than the surrogate gradient.
Another related area of research investigates how to use momentum for the optimization of deep neural networks. Applying different kinds of momentum has become one of the standard tools in current deep learning and it has been shown to speed-up learning in a very wide range of tasks (20; 16; 17). This hints, that for many problems the higher-order terms in deep learning models are "well-behaved" and thus, the gradients do not change too drastically after parameter updates. While these approaches use momentum for parameter updates, our approach can be seen as a form of momentum when sampling directions from the search space of ES.

heading: 3 GRADIENT ESTIMATION
text: We aim at minimizing a function f : Rn → R by steepest descent. In scenarios where the gradient 5f does not exist or is inefficient to compute, we are interested in obtaining some estimate of the (smoothed) gradient of f that provides a good parameter update direction.

heading: 3.1 THE ES GRADIENT ESTIMATOR
text: ES considers the function fσ that is obtained by Gaussian smoothing
fσ(θ) = E ∼N (0,I)[f(θ + σ )] ,
where σ is a parameter modulating the size of the smoothing area and N (0, I) is the n-dimensional Gaussian distribution with 0 being the all 0 vector and I being the n-dimensional identity matrix. The gradient of fσ with respect to parameters θ is given by
5fσ = 1
σ E ∼N (0,I)[f(θ + σ ) ],
which can be sampled by a Monte Carlo estimator, see (5). Often antithetic sampling is used, as it reduces variance (5). The antithetic ES gradient estimator using P samples is given by
gES = 1
2σP P∑ i=1 ( f(θ + σ i)− f(θ − σ i) ) i , (1)
where i are independently sampled from N (0, I) for i ∈ {1, . . . , P}. This gradient estimator has been shown to be effective in RL settings (4).

heading: 3.2 OUR ONE STEP GRADIENT ESTIMATOR
text: We first give some intuition before presenting our gradient estimator formally. Given one surrogate gradient direction ζ, our one step gradient estimator applies the following sampling strategy. First, it estimates how much the gradient points into the direction of ζ by antithetically evaluating f in the direction of ζ. Second, it estimates the part of5f that is orthogonal to ζ by evaluating random, pairwise orthogonal search directions that are orthogonal to ζ. In this way, our estimator detects the optimal lengths of the parameter update step into both the surrogate direction and the evaluated orthogonal directions (e.g. if ζ and 5f are parallel, the update step is parallel to ζ, and if they are orthogonal the step into direction ζ has length 0). Additionally, if the surrogate direction and the gradient are not perfectly aligned, then the gradient estimate almost surely improves over the surrogate direction due to the contribution from the evaluated directions orthogonal to ζ. In the following we define our estimator formally and prove that the estimated direction possesses best possible alignment with the gradient that can be achieved with our sampling scheme.
We assume that k pairwise orthogonal surrogate gradient directions ζ1, . . . , ζk are given to our estimator. Denote by Rζ the subspace of Rn that is spanned by the ζi, and by R⊥ζ the subspace that is orthogonal to Rζ . Further, for vectors v and5f , we denote by v̂ and 5̂f the normalized vector v ‖v‖ and 5f ‖5f‖ , respectively. Let ̂
1, . . . , ̂P be random orthogonal unit vectors from R⊥ζ . Then, our estimator is defined as
g⊥ = k∑ i=1 f(θ + σζ̂i)− f(θ − σζ̂i) 2σ ζ̂i + P∑ i=1 f(θ + σ̂i)− f(θ − σ̂i) 2σ ̂i . (2)
We write5f = 5f‖ζ +5f⊥ζ , where5f‖ζ and5f⊥ζ are the projections of5f on Rζ and R⊥ζ , respectively. In essence, the first sum in (2) computes5f‖ζ by assessing the quality of each surrogate gradient direction, and the second sum estimates5f⊥ζ similar to an orthogonalized antithetic ES gradient estimator, that samples directions from R⊥ζ , see (5). We remark that we require pairwise orthogonal unit directions ̂i for the optimality proof. Due to the orthogonality of the directions, no normalization factor like the 1/P factor in (1) is required in (2). In practice, the dimensionality n is often much larger than P . Then, sampling pairwise orthogonal unit vecotrs i is nearly identical to sampling the is from a N (0, I) distribution, because in high-dimensional space the norm of i ∼ N (0, I) is highly concentrated around 1 and the cosine of two such random vectors is highly concentrated around 0.
For the sake of analysis, we assume that f is differentiable and that the second order approximation f(x + ) ≈ f(x) + 〈 ,5f(x)〉 + TH(x) , where H(x) denotes the Hessian matrix of f at x, is exact. This assumption implies that
f(θ + σ̂)− f(θ − σ̂) 2σ = 〈5f(θ), ̂〉 , (3)
because the even terms cancel for antithetic sampling. The following proposition and theorems provide theoretical understanding, how our gradient estimation scheme improves gradient estimation
in this smooth, noise-free setting. In the following, we will omit the θ in5f(θ). Our first proposition states that g⊥ computes the direction in the subspace spanned by ζ1, . . . , ζk, 1, . . . , P that is most aligned with5f . Proposition 1 (Optimality of g⊥). Let ζ1, . . . , ζk, 1, . . . , P be pairwise orthogonal vectors in Rn. Then, g⊥ = ∑k i=1〈5f, ζ̂i〉ζ̂i + ∑P i=1〈5f, ̂i〉̂i computes the projection of 5f on the subspace spanned by ζ1, . . . , ζk, 1, . . . , P . Especially, = g⊥ is the vector of that subspace that maximizes the cosine 〈5̂f, ̂〉 between5f and . Moreover, the squared cosine between g⊥ and5f is given by
〈5̂f, ĝour〉2 = k∑ i=1 〈5̂f, ζ̂i〉2 + P∑ i=1 〈5̂f, ̂i〉2 . (4)
We remark that when evaluating 〈5f, vi〉 for arbitrary directions vi, no information about search directions orthogonal to the subspace spanned by the vis is obtained. Therefore, one can only hope for finding the best approximation of 5f lying within the subspace spanned by the vis, which is accomplished by g⊥. The proof of Proposition 1 follows easily from the Cauchy-Schwarz inequality and is given in the appendix.

heading: 3.3 ITERATIVE GRADIENT ESTIMATION USING PAST DESCENT DIRECTIONS
text: Our gradient estimation algorithm iteratively applies the one step gradient estimator g⊥ by using the gradient estimate of the last time step as surrogate direction for the current time step. Therefore, our algorithm relies on the assumption that gradients are correlated from one time step to the next. This assumption is justified since it is one of the reasons momentum-based optimizers (15; 16; 17) are successful in deep learning. We also explicitly test this assumption experimentally, see Figure 1a. The quality of the gradient estimate at any time step, depends on the quality of the surrogate gradient, which might be restrained because (a) the previous gradient estimate might be very noisy, and (b) the gradient changes with parameter updates. Our algorithm efficiently tackles problem (a), since it improves the gradient estimate over the surrogate gradient at any time step. In order to theoretically quantify this, we first analyse the case where the gradient does not change with parameter updates, i.e. if a linear function is optimized. In this setting, we show that our algorithm needs a small factor more than n (dimension) samples to align with the true gradient, see Theorem 1, which is close to optimal, because at least n samples are required to determine the true gradient. Next, we incorporate problem (b) into our analysis by considering general, non-linear functions, for which the gradient changes with parameter updates. We show under some simplifying assumptions that, also in this case, our algorithm builds up an improved gradient estimate over time, see Theorem 2.
We first need some notation. Denote by θt the search point, by 5ft = 5f(θt) the gradient and by ζt the parameter update step at time t, that is, θt+1 = θt + ζt. The iterative gradient estimation algorithm obtains the gradient estimate ζt by computing g⊥ with the last update direction ζt−1 as surrogate gradient and P new random directions ̂i. Formally, let ̂1t , . . . , ̂ P t be pairwise orthogonal unit directions chosen uniformly from the unit sphere, that is, they are conditioned to be pairwise orthogonal and are marginally uniformly distributed. By defining t = ∑P i=1〈5ft, ̂it〉̂it, and setting ζt = g⊥, we obtain
ζt = 〈5ft, ζ̂t−1〉ζ̂t−1 + P∑ i=1 〈5ft, ̂it〉̂it = 〈5ft, ζ̂t−1〉ζ̂t−1 + 〈5ft, ̂t〉̂t , (5)
where we used t = ‖ t‖2/‖ t‖ · ̂t = 〈5ft, t〉/‖ t‖ · ̂t = 〈5f, ̂t〉̂t. Then, Equation 4 of Proposition 1 turns into
〈5̂ft, ζt〉2 = 〈5̂ft, ζ̂t−1〉2 + 〈5̂ft, ̂t〉2 . (6)
The next theorem quantifies how fast the cosine between ζt and5ft converges to 1, if5ft does not change over time. Theorem 1 (Convergence rate for linear functions). Let ζt be iteratively computed, as in Equation (5), using the past update direction and P pairwise orthogonal random directions and let Xt = 〈5̂f, ζ̂t〉 be the random variable that denotes the cosine between ζt and 5ft at time t. Then, the expected
drift of X2t is E[X2t −X2t−1|Xt−1 = xt−1] = (1− x2t−1) PN−1 . Moreover, let > 0 and define T to be the first point in time t with X2t ≥ 1− δ. It holds
E[T ] ≤ N − 1 P min{(1− δ)/δ, 1 + ln(1/δ)} .
The first bound E[T ] ≤ N−1P 1−δ δ is tight for δ close to 1 and follows by an additive drift theorem, while the second bound E[T ] ≤ N−1P (1 + ln(1/δ)) is tight for δ close to 0 and follows by a variable drift theorem, see appendix. We remark that in a smooth and noise-free setting, where one can sample the true directional derivative with Equation (3), a cosine squared of 1−δ can be reached by sampling (1 − δ)N random orthogonal directions, see Proposition 2 in the appendix. Since our algorithm evaluates P + 1 directions per time step, it requires approximately min{1/δ, 1+ln(1/δ)1−δ } times more samples to reach the same alignment.
Naturally, the linear case is not the most interesting one. However, it is hard to rigorously analyse the case of general f , because it is unpredictable how the gradient5ft differs from5ft−1. Note that 5ft−5ft−1 ≈ Hζt−1, where H is the Hessian matrix of f at θt−1. We define αt = 〈5̂ft, 5̂ft−1〉 and write 5̂ft = αt5̂ft−1 +5f⊥ where5f⊥ is orthogonal to 5̂ft and has squared norm 1− α2t . Then, the first term of (6) is equal to
〈5̂ft, ζ̂t−1〉2 = 〈αt5̂ft−1 +5f⊥, ζ̂t−1〉2 = ( αt〈5̂ft−1, ζ̂t−1〉+ 〈5f⊥, ζ̂t−1〉 )2 .
In the following, we assume that5f⊥ is a direction orthogonal to5ft−1 chosen uniformly at random. Though, this assumption is not entirely true, it allows to get a grasp on the approximate cosine that our estimator is going to converge to. Theorem 2. Let ζt be iteratively computed using the past update direction and P pairwise orthogonal random directions, see Equation (5), and let Xt = 〈5̂f, ζ̂t〉 be the random variable that denotes the cosine between ζt and5ft at time t. Further, write 5̂ft = αt5̂ft−1+5f⊥, where 1 ≥ αt ≥ 0 and assume that5f⊥ has a random direction orthogonal to 5̂ft−1. Choose ζt according to Equation (5) and define Xt to be the cosine between 5̂ft and ζ̂t. Then,
E[X2t |Xt−1 = xt−1] = ( α2tx 2 t−1 + (1− α2t )(1− x2t−1) 1
N − 1
)( 1− P
N − 1
) + P
N − 1 .
The last theorem implies that the evolution of the cosine depends heavily on the cosine αt between consecutive gradients. Let A = (1−α2t ) 1N−1 (1− P N−1 )+ P N−1
1−(α2t+(1−α2t ) 1N−1 )(1− P N−1 ) . Then, the theorem implies that the drift
E[X2t −X2t−1|Xt−1 = xt−1] is positive if xt−1 ≤ A and negative otherwise. Thus, if αt would not change over time, we would expect Xt to converge to A.
TODO:explain comparison to practise plot: The predicted alignement of the gradient estimate with the true gradient of Theorem 2 is plotted in Figure 1b. Though the assumptions of Theorem 2 are technically not true, the close fit of the theoretical prediction gives empirical evidence that the analysis captures the general behaviour of the gradient estimation process.

heading: 4 EXPERIMENTS
text: In this section, we will empirically evaluate the performance of our gradient estimation scheme when combined with deep neural networks. In Section 4.1, we show that it significantly improves gradient estimation for digit classifiers on MNIST. In Section 4.2, we suggest how to overcome issues that arise from function evaluation noise. Finally, in Section 4.3, we evaluate our gradient estimation scheme on RL environments and investigate further issues arising in this setting.

heading: 4.1 GRADIENT ESTIMATION AND PERFORMANCE ON MNIST
text: We observe that our approach significantly improves gradient estimation compared to standard ES. Figure 1a shows that the key requisite of our iterative gradient estimation scheme is satisfied during
training on MNIST, that is, that gradients between consecutive parameter update steps are correlated. Figure 1b shows that our approach improves gradient estimation compared to ES during the whole training process and strongly improves it in the beginning of training, where consecutive gradients are most correlated, see Figure 1a. We observe that our approach strongly outperforms ES in convergence speed and reaches better final performance for all hyperparameters we tested, see Figure 2 and Table 1.
Implementation details: For these experiments, we used a fully connected neural network with two hidden layers with a tanh non-linearity and 1000 units each, to have a high dimensional model (∼ 1.8 million parameters) . For standard ES 128 random search directions are evaluated at each step. For our algorithm the previous gradient estimate and 126 random search directions are evaluated. We evaluated all directions on the same batch of images in order eliminate function evaluation noise and we resampled after every update step. We used small parameter perturbations (σ = 0.001). This is possible because no function evaluation noise is present and because the objective function is already differentiable and therefore no smoothing is required. We test both SGD and Adam optimizers with learning rates in the range 100.5, 100, . . . , 10−3.

heading: 4.2 ROBUSTNESS TO FUNCTION EVALUATION NOISE
text: In practice, our iterative gradient estimation scheme may suffer from function evaluation noise because it builds up good gradient estimates over several parameter update steps. Suppose that the past update direction is a good descent direction but it performs poorly on the current batch used for evaluation due to randomness in the batch selection or network evaluation process. Then, this direction is weighted lightly when computing the new update direction, see Equation 5, and therefore the information about this direction will be discarded. We empirically show, that our approach suffers heavily from this issue when artificially injecting noise in the function evaluation process, see Figure 3b . Figure 3c shows, that this issue can be resolved by using the last k update directions for our gradient estimator (see Equation 2). In this case, a good direction is only discarded, when it performs poorly in k consecutive evaluation steps, which is very unlikely. We remark that the magnitude of the parameter updates naturally limits k, because the k-th last update direction is only useful if it is still correlated with the current gradient. Concretely, we found that using the last 4 parameters updates was extremely helpful for smaller learning rates, even in the absence of noise (see Figure 3c). However, it did not offer an advantage for larger learning rates.

heading: 4.3 ROBOTIC RL ENVIRONMENTS
text: For the next set of experiments, we evaluate our algorithm on three robotics tasks of the Roboschool environment: RoboschoolInvertedPendulum-v1, RoboschoolHalfCheetah-v1 and RoboschoolAnt-v1. Our approach outperforms ES in the pendulum task, and offers a small improvement over ES in the other two tasks, see Figure 4. The improvement of our approach over standard ES is smaller on RL tasks than on the MNIST task. Therefore, we first empirically confirmed that past updated direction are also in RL correlated with the gradient. To test this, we kept track of the average difference between random perturbations and the direction given by our algorithm, after normalizing the rewards.
We found that, the direction of our algorithm had an average weight of 1.11 versus the 0.65 of a random direction.
RL robotics tasks bring two additional major challenges compared to the MNIST task. First, exploration is crucial to escape local optima and find new solutions, and second, the function evaluation noise is huge due to each perturbation being tested only on a single trajectory. Our proposed solution of robustness against function evaluation noise intertwines with the exploration issue. A rather small step size is necessary in order to use more past directions as surrogate gradients. However, exploration in ES is driven by large perturbation sizes and noisy optimization trajectories. We did not observe improvements when combining the approach of using several past directions with standard hyperparameter settings. We believe that an exhaustive empirical study can shed light onto the effect of our approach on exploration and may further improve the performance on RL tasks. However, running extensive experiments for complex RL environment is computationally expensive.
Implementation details: We use most of the hyper-parameters from the OpenAI implementation 1 . That is, two hidden layers of 256 units each and with a tanh non-linearity. Further, we use a learning rate of 0.01, a perturbation standard deviation of σ = 0.02 and the Adam optimizer, and we also apply fitness shaping (19). For standard ES 128 random perturbations are evaluated at each step. For our algorithm the previous gradient estimate and 126 random perturbations are evaluated. For the Ant and Cheetah environments, we observed with this setup, that agents often get stuck in a local optima where they stay completely still, instead of running forward. As this happens for both, ES and our algorithm, we tweaked the environments in order to ensure that a true solution to the task is learned and not some some degenerate optima, we tweaked the environments in the following way. We remove the penalty for using electricity and finish the episode if the agent does not make any progress in a given amount of time. In this way, agents consistently escape the local minima. We use a tanh non-linearity on the output of the network, which increased stability of training, as otherwise the output of the network would become very large without an electricity penalty.

heading: 5 CONCLUSION
text: We proposed a gradient estimator that optimally incorporates surrogate gradient directions and random search directions, in the sense that it determines the direction with maximal cosine to the true gradient from the subspace of evaluated directions. Such a method has many applications as elucidated in (14). Importantly, our estimator does not require information about the quality of surrogate directions, which allows us to iteratively use past update directions as surrogate directions for our gradient estimator. We theoretically quantified the benefits of the proposed iterative gradient estimation scheme. Finally, we showed that our approach in combination with deep neural networks considerably improves the gradient estimation capabilities of ES, at no extra computational cost. The results on MNIST indicate that the speed of the Evolutionary Strategies themselves, a key part in the current Reinforcement Learning toolbox, is greatly improved. Within Reinforcement Learning an out of the box application of our algorithm yields some improvements. The smaller improvement in RL compared to MNIST is likely due to the interaction of our approach and exploration that is essential in RL. We leave it to future work to explicitly add and study appropriate exploration strategies which might unlock the true potential of our approach in RL.
1https://github.com/openai/evolution-strategies-starter

heading: A PROOF OF THEOREMS
text: In this Section we prove the theorems from the main paper rigorously.
A.1 PROOF OF PROPOSITION 1
Note that for this theorem there is no distinction between the directions ζi and i. For ease of notation, we denote ζ1, . . . , ζk, ζ1, . . . , ζP by 1, . . . , m. The theorem is a simple application of the CauchySchwarz inequality. Denote by 5f‖ i = ∑m i=1〈5f, ̂i〉̂i the projection of 5f on the subspace
spanned by the is, and let = ∑m i=1 αî
i be a vector in that subspace. Then, the Cauchy-Schwarz inequality implies
〈5f, 〉 = 〈5f‖ i , 〉 ≤ ‖ 5 f‖ i‖‖ ‖ . (7)
Equality holds if and only if and5f‖ i have the same direction, which is equivalent to = αgour for some α > 0. In particular, in this case the cosine squared between5f and is
〈5̂f, ̂〉2 = ‖ 5 f‖ i‖2
‖ 5 f‖2 = m∑ i=1 〈5̂f, ̂i〉2 (8)
A.2 EXPECTATION OF COSINE SQUARED OF RANDOM VECTORS FROM THE UNIT SPHERE
We need the following proposition for the proofs of Theorem 1 and 2.
Proposition 2. Let û be an N -dimensional unit vector, and let ̂1, . . . , ̂P be pairwise orthogonal vectors sampled uniformly from the N -dimensional unit sphere, that is, they are marginally uniformly distributed and conditioned to be pairwise orthogonal. Then, the expected cosine squared of û and = ∑P i=1〈u, ̂i〉̂i is
E[〈û, ̂〉2] = P N
Proof. Note that
〈û, ̂〉2 = 1 ‖ ‖2 〈û, 〉2 = 1∑P
i=1〈û, ̂i〉2 ( P∑ i=1 〈û, ̂i〉2 )2 = P∑ i=1 〈û, ̂i〉2
Denote by S the N -dimensional unit sphere. Linearity of expectation implies that
E[〈û, ̂〉2] = P∑ i=1 E[〈û, ̂i〉2]
= P
V ol(S) ∫ S 〈û, v〉2 dv
By rotational invariance of the unit sphere, we can replace û by e1 = (1, 0, . . . , 0) and obtain
E[〈û, ̂〉2] = P V ol(S) ∫ S 〈e1, v〉2 dv
= P
V ol(S) ∫ S v21 dv
= P
N · V ol(S) ∫ S N∑ i=1 v2i dv
= P
N · V ol(S) ∫ S 1 dv
= P
N
A.3 PROOF OF THEOREM 1
In order to prove Theorem 2, use Equation 6 to compute how X2t = 〈5̂f, ζ̂t〉2 depends on X2t−1. Then, we apply a variable transformation to Xt in order to be able to apply the additive and variable drift theorems from (21), which are stated in Section B.
We can split the normalized gradient 5̂f = 5f⊥ζ̂t−1 +5f‖ζ̂t−1 into an orthogonal to ζ̂t−1 part and a parallel to ζ̂t−1 part. It holds ‖ 5 f⊥ζ̂t−1‖
2 = 1− 〈5̂f, ζ̂t−1〉2 and 〈5̂f‖ζ̂t−1 , ̂t〉 = 0 since ̂t is a unit vector orthogonal to ζ̂t−1. Recall that 5̂f⊥ζ = 5f⊥ζ‖5f⊥ζ‖ , then
〈5̂f, ̂t〉2 = 〈5f⊥ζ̂t−1 , ̂t〉 2 = (1− 〈5̂f, ζ̂t−1〉2)〈5̂f⊥ζ , ̂t〉2 = (1−X2t−1)〈5̂f⊥ζ , ̂t〉2 , (9)
and therefore by Equation 6
X2t = 〈5̂f, ζ̂t−1〉2 + 〈5̂f, ̂t〉2 = X2t−1 + (1−X2t−1)〈5̂f⊥ζ , ̂t〉2 . (10)
Define the random process Yt = 1−X2t . It holds
Yt = 1−X2t−1 + (1−X2t−1)〈5̂f⊥ζ , ̂t〉2 = Yt−1(1− 〈5̂f⊥ζ , ̂t〉2) , and (11) E[Yt|Yt−1 = yt−1] = yt−1 ( 1− E[〈5̂f⊥ζ , ̂t〉2] ) = yt−1 ( 1− P
N − 1
) , (12)
where we used Proposition 2 in the N − 1 dimensional subspace that is orthogonal to ζt−1. In order to derive the first bound on T , we bound the drift of Yt for Yt ≥ δ.
E[Yt|Yt−1 = yt−1, yt−1 ≥ δ] = yt−1 − yt−1 P
N − 1 ≤ yt−1 − δ
P
N − 1 ,
where we used Equation 12 and yt−1 ≥ δ. In order to apply Theorem 3, we define the auxiliary process Zt = Yt − δ. Then, T is the expected time that Zt hits 0. Since E[Zt|Zt−1 = zt−1, zt−1 ≥ 0] ≤ zt−1 − δ cN−1 and Z0 = 1− δ, Theorem 3 implies that
E[T ] ≤ 1− δ δ N − 1 P
In order to apply Theorem 4, to show the second bound on T , we need to rescale Yt such that it takes values in {0} ∪ [1,∞). Define the auxiliary process Zt by
Zt = { Yt/δ if Yt ≥ δ 0 if Yt < δ . (13)
Then, T is the expected time that Zt hits 0. The process Zt satisfies E [Zt|Zt−1 = zt−1, zt−1 ≥ 1] ≤ E [Yt/δ|Yt−1 = δzt−1, zt−1 ≥ 1] ≤ zt−1 ( 1− P
N − 1
) ,
where we used Equations 13 and 12. Since Z0 = 1/δ, Theorem 4 implies for h(z) = z cN−1 that
E[T ] ≤ N − 1 P + ∫ 1/δ 1 N − 1 Pu du = N − 1 P (1 + ln(1/δ)) .
A.4 PROOF OF THEOREM 2
In order to prove the theorem, we need to understand how Xt depends on the value of Xt−1. It holds X2t = 〈5̂ft, ζ̂t−1〉2 + 〈5̂ft, ̂t〉2. As in Equation 12, we can write 〈5̂ft, ̂t〉2 = (1 − 〈5̂ft, ζ̂t−1〉2)〈5̂f⊥ζ̂t−1 , ̂t〉 2, and note that E[〈5̂f⊥ζ̂t−1 , ̂t〉 2] = PN−1 follows by Proposition 2 and the definition of ̂t. This implies that
E[X2t |Xt−1 = xt−1] = ( 1− P
N − 1
) E[〈5̂ft, ζ̂t−1〉2|Xt−1 = xt−1] + P
N − 1 (14)
To understand how the Xt evolves we need to analyze how 〈5̂ft, ζ̂t−1〉2 relates to Xt−1 = 〈5̂ft−1, ζ̂t−1〉2. To that end, we set αt = 〈5̂ft, 5̂ft−1〉 and write 5̂ft = αt5̂ft−1 +5f⊥5ft−1 where5f⊥5ft−1 is orthogonal to5ft−1 and has norm 1− α2t . Then,
〈5̂ft, ζ̂t−1〉2 = 〈αt5̂ft−1 +5f⊥5ft−1 , ζ̂t−1〉2 = ( αt〈5̂ft−1, ζ̂t−1〉+ 〈5f⊥5ft−1 , ζ̂t−1〉 )2 .
It follows that
E[〈5̂ft, ζ̂t−1〉2|Xt−1 = xt−1] (15) = α2tx 2 t−1 + 2αtxt−1E[〈5f⊥5ft−1 , ζ̂t−1〉] + E[〈5f⊥5ft−1 , ζ̂t−1〉2] (16)
= α2tx 2 t−1 + (1− α2t )(1− x2t )
1
N − 1 , (17)
where we used E[〈5f⊥5ft−1 , ζ̂t−1〉] = 0, which follows from the assumption of5f⊥5ft−1 being a random direction orthogonal to5ft−1, and that
E[〈5f⊥5ft−1 , ζ̂t−1〉2] = E[〈5f⊥5ft−1 , ζ̂t−1⊥5ft−1〉2]
= ‖ 5 f⊥5ft−1‖2‖ζ̂t−1⊥5ft−1‖2E[〈 5f⊥5ft−1 ‖ 5 f⊥5ft−1‖ , ζ̂t−1⊥5ft−1 ‖ζ̂t−1⊥5ft−1‖ 〉2]
= (1− α2t )(1− x2t ) 1
N − 1 ,
which follows from ‖5f⊥5ft−1‖2 = 1−α2, ‖ζ̂t−1⊥5ft−1‖2 = 1−x2t and Proposition 2 for P = 1 using that5f⊥5ft−1 is a random direction orthogonal to5ft−1. Then, plugging Equation (17) into (14), implies the theorem.

heading: B DRIFT THEOREMS
text: For the proof of Theorem 1, we use two drift theorems from (21), which we restate for completeness. Theorem 3 (Additive Drift, Theorem 1 from (21)). Let (Xt)t∈N0 be a Markov chain with state space S ⊂ [0,∞) and assume X0 = n. Let T be the earliest point in time t ≥ 0 such that Xt = 0. If there exists c > 0 such that for all x ∈ S, x > 0 and for all t ≥ 0 we have
E[Xt+1|Xt = x] ≤ x− c . Then,
E[T ] ≤ n c .
Theorem 4. Variable Drift, Theorem 4 from (21)] Let (Xt)t∈N be a Markov chain with state space S ⊂ {0} ∪ [1,∞) and with X0 = n. Let T be the earliest point in time t ≥ 0 such that Xt = 0. Suppose furthermore that there is a positive, increasing function h : [1,∞)→ R>0 such that for all x ∈ S, x > 0 we have for all t ≥ 0
E[Xt+1|Xt = x] ≤ x− h(x) .
Then,
E[T ] ≤ 1 h(1) + ∫ n 1 1 h(u) du .

heading: C ADDITIONAL EXPERIMENTS
text: This section presents some additional experiments. We will incorporate a polished version of this section into the experiments section of the main paper, for the camera ready version of the paper.
C.1 COMPARISON OF OUR APPROACH TO SEPARABLE NES AND GUIDED ES
We consider the toy task of optimizing a quadratic function f(x) = ‖Ax‖22 as considered in (14). Figure 5 compares our approach to a diagonal approximation of CMA-ES, that is, seperable natural ES (SNES) from (19), that has the same runtime and memory complexity as our approach. While SNES works verv well for small diminsional parameter space (N = 32 , see Figrue 5a ), it is clearly outperformed by our approach for high dimensional paramater spaces ( N = 10000, see Figure 5b).
Further, we compare our one step gradient estimator g⊥ to the guided-ES gradient estimator proposed in (14). Again the goal is to optimize a quadratic function f(x) = ‖Ax‖22. Both estimators receive a surrogate gradient ζ for the optimization. The surrogate gradient is created as in (14): A normalized random bias b (drawn once at the beginning of optimization) and a normalized noise direction n (resampled at every iteration) are added to the true gradient. that is, ζ(x) = 5f(x)+(b+n)‖5f(x)‖2 Figure 6 shows that our approach outperforms guided-ES no matter how the parameter α controlling the bias variance trade-off of guided-ES is set.
Implementation details: For Figure 5b, x has dimension N = 10000 and is initalized by sampling x0 ∼ N (0, I). A is a matrix of dimension 100×10000, and its entries are samples from aN (0, 0.001) distribution, where the variance is chosen such that the initial loss is approx 1. For Figure 5a, x has dimension N = 16 and A has dimension 16× 16. The hyperparamters for our approach and standard ES are the learning rate β and perturbation size σ. The hyperparamters for SNES are learning rate
β, learning rate βσ for the step size adaptation and perturbation size initialization. We performed a hyperparameter search for all these hyperparameters and plotted the best performance.

Review summaries:
This paper addresses the issue of noisy gradient estimation in a type of evolution strategies popularized by the open AI's reinforcement learning paper. It is a follow-up paper of reference [14], and try to analyze the optimality of the gradient estimation. The goal of the paper is well stated and well motivated. The paper itself is well-organized. However, the novelty of this work is not sufficiently high and its usefulness is questionable. 

This paper provides a new type of gradient estimator that combines an Evolutionary Strategies (ES) style estimate (using function evaluations at perturbed parameters) along with surrogate gradient estimates (gradient estimates that may be biased and/or high variance). The estimator involves computing antithetic ES estimates in two subspaces: along the set of (normalized) surrogate gradients, and along a set of randomly chosen vectors in the orthogonal complement of the span of the surrogate gradients. The paper provides a proof of the optimality of the estimate, that is, the proposed gradient estimate maximizes the cosine of the angle with the true gradient over the vectors in the subspace defined by the set of surrogate gradients and sampled directions. The paper proposes an additional mechanism for generating surrogate gradients by simply using previous gradient estimates as surrogate gradients, and derives a convergence rate for when this iterative estimator will approximate a fixed, true gradient (e.g. for linear functions). Finally, the paper applies the estimate to two tasks: MNIST classification and robotic control via reinforcement learning, demonstrating improvements on both compared to standard ES.