Accelerating Reinforcement Learning with a Directional-Gaussian-Smoothing Evolution Strategy

ES belongs to the family of blackbox optimization that employs random search techniques to optimize an objective function (Rechenberg & Eigen, 1973; Schwefel, 1977). The search can be guided via either covariance matrix adaptation (Hansen & Ostermeier, 2001; Hansen, 2016), or search gradients to be fed to first order optimization schemes, (Wierstra et al., 2014). The DGS-ES method is in line with the second approach. Gaussian smoothed as a method for approximating search direction has been introduced in (Flaxman et al., 2005; Nesterov & Spokoiny, 2015). In recent years, ES has been revived as a popular strategy in machine learning. In addition to those mentioned above, applications of ES to RL can also be found in Fazel et al. (2018); Ha & Schmidhuber (2018); Houthooft et al. (2018); Müller & Glasmachers (2018); Liu et al. (2019); Pourchot & Sigaud (2019). Applications of ES beyond RL include meta-learning (Metz et al., 2018), optimizing neural network architecture (Real et al., 2017; Miikkulainen et al., 2017; Cui et al., 2018), and direct training of neural network (Morse & Stanley, 2016).

Effective exploration is also critical for deep RL on high dimensional action and state spaces, and a wide variety of exploration strategies have been developed in recent years, to name a few, count-based exploration (Ostrovski et al., 2017; Tang et al., 2017), intrinsic motivation (Bellemare et al., 2016), curiosity (Pathak et al., 2017) and variational information maximization (Houthooft et al., 2016). For exploration techniques which add noise directly to the parameter space of policy, see (Fortunato et al., 2017; Plappert et al., 2017). Also, exploration can be combine with deep RL methods to improve sample efficiency (Colas et al., 2018). Finally, for a recent survey on policy search, we refer to (Sigaud & Stulp, 2019).

Reinforcement learning considers agents that are able to take a sequence of actions in an environment, denoted by $\mathcal{E}$, over a number of discrete time steps $t\in\{1,\ldots,T\}$. At each time step $t$, the agent receives a state $\bm{s}_{t}\in\mathcal{S}$ and produces a follow-up action $\bm{a}_{t}\in\mathcal{A}$. The agent will then observe another state $\bm{s}_{t+1}$ and a scalar reward $r_{t}$. The goal of the agents is to learn a policy $\pi(\bm{a}_{t}|\bm{s}_{t})$ that maximizes the objective function, i.e., the expected cumulative return of the form

where $0\leq\gamma\leq 1$ is a discount rate, $\bm{a}_{t}$ is drawn from the policy $\pi(\bm{a}_{t}|\bm{s}_{t})$, and $\bm{s}_{t+1}=\mathcal{E}(\bm{s}_{t},\bm{a}_{t})$ is generated by running the environment dynamics.

In policy-based reinforcement learning approaches, the policy $\pi(\bm{a}|\bm{s};\bm{\theta})$ is parameterized by $\bm{\theta}\in\mathbb{R}^{d}$, where the vector $\bm{\theta}:=(\theta_{1},\ldots,\theta_{d})^{\top}$ represents the parameters of the policy, e.g., weights of a neural network. Then, the task of learning a good policy $\pi$ becomes iterative updating the parameter $\bm{\theta}$ to solve the following optimization problem

where we denote $J(\bm{\theta}):=J(\pi(\bm{a}|\bm{s};\bm{\theta}))$ by an abuse of notation. In reinforcement learning, it is usually true that the gradient of the environment $\mathcal{S}$ is inaccessible, so automatic differentiation cannot be used to obtain the gradient of $J(\bm{\theta})$. Thus, much of the innovation in reinforcement learning algorithms is focused on addressing the lack of access to or the existence of gradients of the environment/policy. In this work, we focus on the evolution strategy and other types of training algorithms for reinforcement learning are summarized in §1.

We briefly recall the evolution strategy methods, e.g., (Hansen & Ostermeier, 2001; Salimans et al., 2017), which use a multivariate Gaussian distribution to generate the population around the current parameter value $\bm{\theta}_{t}$ at the $t$-iteration. When the Gaussian distribution can be factorized to $d$ independent one-dimensional Gaussian distributions, the standard ES method can be mathematically interpreted based on the Gaussian smoothing technique (Flaxman et al., 2005; Nesterov & Spokoiny, 2015). Specifically, a smoothed version of $J(\bm{\theta})$ in Eq. (1), denoted by $J_{\bm{\sigma}}(\bm{\theta})$, is defined by

where $\mathcal{N}(0,\mathbf{I}_{d})$ is a $d$-dimensional standard Gaussian distribution, the notation “$\circ$” represents the element-wise product, and the vector $\bm{\sigma}=(\sigma_{1},\ldots,\sigma_{d})$ controls the smoothing effect. It is well known that $J_{\bm{\sigma}}$ is always differentiable even if $J$ is not. In addition, most of the characteristics of the original objective function $J(\bm{\theta})$, e.g., convexity, the Lipschitz constant, are inherited by $J_{\bm{\sigma}}(\bm{\theta})$. When $J(\bm{\theta})$ is differentiable, the difference $\nabla J-\nabla J_{\bm{\theta}}$ can be bounded by its Lipschitz constant (see (Nesterov & Spokoiny, 2015), Lemma 3 for details). Thus, the original optimization problem in Eq. (2) can be replaced by a smoothed version, i.e.,

where the gradient of $J_{\bm{\sigma}}(\bm{\theta})$ is given by

The standard ES method (Salimans et al., 2017) uses Monte Carlo sampling to estimate the gradient $\nabla J_{\bm{\sigma}}(\bm{\theta})$ and update the state $\bm{\theta}$ from iteration $n$ to $n+1$ by

where $\lambda$ is the learning rate, $\bm{u}_{m}$ are sampled from the Gaussian distribution $\mathcal{N}(0,\mathbf{I}_{d})$.

One drawback of the ES method and its variants is the slow convergence of the training process, due to the low accuracy of the MC-based gradient estimator (see (Berahas et al., 2019)), also (Zhang et al., 2020), for extended discussions on the accuracy of gradient approximations using Eq. (6) and related methods). On the other hand, the evaluations of $J(\bm{\theta}_{n}+\bm{\sigma}\circ\bm{u}_{m})$ for $m=1,\ldots,M$ at the $n$-th iteration can be generated totally in parallel, which makes it well suited to be scaled up to a large number of parallel workers on modern supercomputers. Therefore, the motivation of this work is to develop a new gradient operator to replace the one in Eq. (5), such that the new gradient can be approximated in a much more accurate way and the embarrassing parallelism feature can be retained.

For a given direction $\bm{\xi}\in\mathbb{R}^{d}$, the restriction of the objective function $J(\bm{\theta})$ along $\bm{\xi}$ can be represented by

where $\bm{\theta}$ is the current state of the agent’s parameters. Then, we can define the one-dimensional Gaussian smoothing of $G(y)$, denoted by $G_{\sigma}(y)$, by

which is also the Gaussian smoothing of $J(\bm{\theta})$ along $\bm{\xi}$ in the neighbourhood of $\bm{\theta}$. The derivative of $G_{\sigma}(y|\bm{\theta},\bm{\xi})$ at $y=0$ is given by

where $\mathscr{D}$ denotes the differential operator. It is easy to see that $\mathscr{D}[G_{\sigma}(0\,|\,\bm{x},\bm{\xi})]$ only involves the directionally smoothed objective function given in Eq. (8).

We can assemble a new gradient, i.e., the DGS gradient, by putting together the derivatives in Eq. (9) along $d$ orthogonal directions, i.e.,

where $\bm{\Xi}:=(\bm{\xi}_{1},\ldots,\bm{\xi}_{d})^{\top}$ represents the matrix consisting of $d$ orthonormal vectors. It is important to notice that

for any $\bm{\sigma}>0$, because of the directional Gaussian smoothing used in Eq. (10). However, there is consistency between the two quantities as $\bm{\sigma}\rightarrow 0$, i.e.,

for fixed $\bm{\theta}$ and $\bm{\Xi}$. If $\nabla J(\bm{\theta})$ exists, then ${\nabla}_{\bm{\sigma},\bm{\Xi}}[J](\bm{\theta})$ will also converge to $\nabla J(\bm{\theta})$ as $\bm{\sigma}\rightarrow 0$. Such consistency naturally led to the idea of replacing $\nabla J_{\bm{\sigma}}(\bm{\theta})$ with ${\nabla}_{\bm{\sigma},\bm{\Xi}}[J](\bm{\theta})$ in the ES framework.

Since each component of ${\nabla}_{\bm{\sigma},\bm{\Xi}}[J](\bm{\theta})$ in Eq. (10) only involves a one-dimensional integral, we can use Gaussian quadrature rules (Quarteroni et al., 2007) to obtain spectral convergence. In the case of Gaussian smoothing, a natural choice is the Gauss-Hermite (GH) rule, which is used to approximate integrals of the form $\int_{\mathbb{R}}g(x){\rm e}^{-x^{2}}dx$. By doing a simple change of variable in Eq. (9), the GH rule can be directly used to obtain the following estimator:

where $w_{m}$ are the GH quadrature weights defined by

$v_{m}$ are the roots of the Hermite polynomial of degree $M$

and $M$ is the number of function evaluations, i.e., environment simulations, needed to compute the quadrature in Eq. (12). The weights $\{w_{m}\}_{m=1}^{M}$ and the roots $\{v_{m}\}_{m=1}^{M}$ can be found in Abramowitz & Stegun (1972). The approximation error of the GH formula can be bounded by $\widetilde{\mathscr{D}}^{M}[G_{\sigma}]$ is

where $M!$ is the factorial of $M$, and the constant $C_{0}>0$ is independent of $M$ and $\sigma$.

Applying the GH quadrature rule $\widetilde{\mathscr{D}}^{M}$ to each component of ${\nabla}_{\bm{\sigma},\bm{\Xi}}[J](\bm{\theta})$ in Eq. (10), we define the following estimator:

which requires a total of $M\times d$ parallelizable environment simulations at each iteration of training. The error bound in Eq. (15) indicates that $\widetilde{\nabla}^{M}_{\bm{\sigma},\bm{\Xi}}[J](\bm{\theta})$ is an accurate estimator of the DGS gradient for a small $M$, regardless of the value of $\bm{\sigma}$. This enables the use of relatively big values of $\bm{\sigma}$ in Eq. (16) to exploit the nonlocal features of the landscape of $J(\bm{\theta})$ in the training process. The nonlocal exploitation ability of $\widetilde{\nabla}^{M}_{\bm{\sigma},\bm{\Xi}}[J]$ is demonstrated in §4 to be effective in reducing the necessary number of iterations to achieve a prescribed reward score.

On the other hand, as the quadrature weights $w_{m}$ and $v_{m}$ defined in Eq. (13) and Eq. (14), are deterministic values, the DGS estimator $\widetilde{\nabla}^{M}_{\bm{\sigma},\bm{\Xi}}[J](\bm{x})$ in Eq. (16) is also a deterministic for fixed $\bm{\Xi}$ and $\bm{\sigma}$. To introduce random exploration ability to our approach, we add random perturbations to both $\bm{\Xi}$ and $\bm{\sigma}$. For the orthonormal matrix $\bm{\Xi}$, we add a small random rotation, denoted by $\Delta\bm{\Xi}$, to the current matrix $\bm{\Xi}$. The matrix $\Delta\bm{\Xi}$ is generated as a random skew-symmetric matrix, of which the magnitude of each entry is smaller than a prescribed threshold $\alpha>0$. The perturbation of $\bm{\sigma}$ is conducted by drawing random samples from a uniform distribution $U(r-\beta,r+\beta)$ with two hyper-parameters $r$ and $\beta$ with $r-\beta>0$. The random perturbation of $\bm{\Xi}$ and $\bm{\sigma}$ can be triggered by various types of indicators e.g., the magnitude of the DGS-gradient, the number of iteration done since last perturbation.

We compared Algorithm 1 against several RL baselines, including ES, PPO and TRPO, as well as the state-of-the-art algorithms such as ASEBO, DDPG, and TD3. Below is the information of the packages used in this effort.

• •

  ES: The Evolution Strategy proposed in Salimans et al. (2017). We used the implementation of ES from the open-source code https://github.com/hardmaru/estool.

• •

  ASEBO: Adaptive ES-Active Subspaces for Blackbox Optimization, which was recently developed by Choromanski et al. (2019a). We used the implementation released by the authors at https://github.com/jparkerholder/ASEBO.

• •

  PPO: Proximal Policy Optimization in Schulman et al. (2017), which is available in OpenAI’s baselines repository at https://github.com/openai/baselines (Dhariwal et al., 2017).

• •

  TRPO: Trust Region Policy Optimization, developed by Schulman et al. (2015). We also used the OpenAI’s baselines implementation (Dhariwal et al., 2017).

• •

  DDPG: Deep Deterministic Policy Gradient, proposed by Lillicrap et al. (2016). We used the implementation from https://github.com/georgesung/TD3 where the benchmark DDPG in PyBullet is provided.

• •

  TD3: Twin Delayed Deep Deterministic Policy Gradient (Fujimoto et al., 2018), which was built upon the DDPG. The original results were reported for the MuJoCo version environments using the implementation from https://github.com/sfujim/TD3, but we used the PyBullet implementation from https://github.com/georgesung/TD3.

The hyper-parameters for all algorithms above were set to match the original papers without further tuning to improve performance on the testing benchmark examples.

Figure 1 shows the comparison results of CartPole, Pendulum and MountainCar problems from the OpenAI Gym. We compared the DGS-ES with classical ES and the improved ASEBO method. In general, the DGS-ES method features faster convergence than the baselines. For the simplest CartPole problem, the three methods perform equally well. Discrepancy appears in the Pendulum test, where the DGS-ES method not only converges faster than the baselines, but also achieves a higher average return. There is a much bigger discrepancy between the DGS-ES and the baselines appear in the MountainCar test. According to the guideline provided in the OpenAI Gym, the success threshold is to achieve an average return of 90. The DGS-ES method achieves the threshold within 500 iterations, while the average returns of the ES and ASEBO methods are around zero. It is well known that the challenge of this problem is that the surface of the objective function $J(\bm{\theta})$ is very flat at most locations in the parameter space, which makes it difficult to capture the peak of $J(\bm{\theta})$. This test is a good demonstration of the nonlocal exploration ability of the DGS-ES method. Since the mean of the smoothing factor $\bm{\sigma}$ is set to 1.0, DGS-ES can capture the peak of $J(\bm{\theta})$ much faster than the baselines. In fact, we can see that it took around 50 iterations for the DGS-ES to find the peak region, while ES and ASEBO needed more than 3000 iterations (not plotted) to move out of the flat region.

PyBullet library. We compared the DGS-ES with six baselines including ES, ASEBO, PPO, TRPO, DDPG and TD3. As expected, the DGS-ES method shows better performance in terms of the convergence speed. For the Hopper-v0 problem, the DGS-ES achieves the highest return of all the methods within 400 iterations. It is worth pointing out that some of the baselines will eventually catch up with the DGS-ES given a sufficiently large number of iterations. For example, DDPG and TD3 do not provide much improvement within 400 iterations, but DDPG could reach 1650 average return with over 3000 iterations, and TD3 could reach even higher, around 2200 average return, with over 6000 iterations, according to the baselines provided by OpenAI (Dhariwal et al., 2017) and Fujimoto et al. (2018). This phenomenon illustrates the fast convergence feature of DGS-ES. For the InvertedPendulum-v0 problem, DGS-ES can achieve the maximum return 1000 (default value in the PyBellet library) around 30 iterations. In comparison, ES and ASEBO can reach the maximum return but with more iterations than DGS-ES. According to the benchmark results for PyBullet environments in Fujimoto et al. (2018), DDPG and TD3 cannot converge to the maximum return even with a large number of iterations. For the Reacher-v0 problem, DGS-ES and ASEBO are still the top performers, and the advantage of DGS-ES is, again, faster convergence.