Accelerating reinforcement learning with a Directional-Gaussian-Smoothing evolution strategy

1.   Introduction

Reinforcement learning is a class of problems which aim to find, through trial and error, a feedback policy that prescribes how an agent should act in an uncertain, complex environment to maximize some notion of cumulative reward [30]. RL can be viewed as a generalization of stochastic optimal control. Traditionally, RL algorithms have mainly been employed for small input and action spaces, and suffered difficulties when scaling to high-dimensional problems. With the recent emergence of deep learning, powerful non-linear function approximators such as deep neural networks (DNN) can be integrated into RL and extend the capability of RL in a variety of challenging tasks which would otherwise be infeasible, ranging from playing Atari from pixels [18], playing expert-level Go [28] to robotic control [2]. Among the most popular current deep RL algorithms are Q-learning methods, policy gradient methods, and evolution strategies. Deep Q-learning algorithms [18] use a DNN to approximate the optimal Q function, yielding policies that, for a given state, choose the action that maximizes the Q-value. Policy gradient methods [26] improve the policies with a gradient estimator obtained from sample trajectories in action space, examples of which are A3C [17], TRPO [24] and PPO [25].

This work concerns RL techniques based on evolution strategies. ES refers to a family of blackbox optimization algorithms inspired by ideas of natural evolution, often used to optimize functions when gradient information is inaccessible. This is exactly the prominent challenge in a typical RL problem, where the environment and policy are usually nonsmooth or can only be accessed via noisy sampling. It is not totally surprising ES has become a convincing competitor to Q-learning and policy gradient in deep RL in recent years. Unlike policy gradient, ES perturbs and performs policy search directly in the parameter space to find an effective policy, which is now generally considered to be superior to action perturbation [27]. The policy search can be guided by a surrogate gradient [23], completely population-based and gradient-free [29], and hybridized with other exploration strategies such as novelty search and quality diversity [7]. It has been shown in those works that ES is easy to parallelize, and requires low communication overhead in a distributed setting. More importantly, being able to achieve competitive performance on many RL tasks, these methods are advantageous over other RL approaches in their high scalability and substantially lower wall-clock time. Given wide availability of distributed computing resources, all environment simulations at each iteration of training can be conducted totally in parallel. Thus, a more reasonable metric for the performance of a training algorithm is the number of non-parallelizable iterations, as opposed to the sample complexity. In this metric, ES is truly an appealing choice.

Nevertheless, there are several challenges that need to be addressed in order to further improve the performance of ES in training complex policies. First, most ES methods cope with the non-smoothness of the objective function by considering a Gaussian-smoothed version of the expected total reward. The gradient of this function is intractable and must be estimated to provide the policy parameter update. In the pioneering work [23], a gradient estimator is proposed based on random parameter sampling. Developing efficient sampling strategies for gradient estimates has become an interest in ES research since then, and several improvements have been proposed, based on imposing structures on parameter perturbation [5], or reusing past evaluations, [6,15,16]. Yet, most of these gradient estimators are of Monte Carlo type, therefore arguably affected by the low accuracy of Monte Carlo methods. For faster convergence of training (i.e., reducing the number of iterations), more accurate gradient estimators are desired, particularly in RL tasks where the policy has a large number of parameters to learn. Another prominent challenge is that the landscape of the objective function is complex and possesses plentiful local maxima. There is a risk for any optimization algorithm to get trapped in some of those points and unable to explore the parameter space effectively. The Gaussian smoothing, with its ability to smooth out function and damps out small, insignificant fluctuations, is a strong candidate in this very challenge. Specifically, with a moderately large smoothing parameter (i.e., strong smoothing effect), we can expect the gradient of the smoothed objective function will be able to look outside unimportant variation in the adjacent area and detect the general trend of the function from a distance, therefore an efficient nonlocal search direction. This potential of Gaussian smoothing, however, has not been explored in reinforcement learning.

In this paper, we propose a new strategy to accelerate the time-to-solution of reinforcement learning by exploiting the Directional Gaussian Smoothing Evolution Strategy (DGS-ES) method, developed in our recently work [33]. The DGS-ES method introduced a new directional Gaussian smoothing (DGS) gradient operator, that smooths the original objective function only along $d$-orthogonal directions in the parameter space. In other words, the DGS gradient requires $d$ one-dimensional Gaussian convolutions, instead of one $d$-dimensional convolution in the existing ES methods. There are several advantages of using the DGS gradient operator in reinforcement learning. First, each component of the DGS gradient, represented as a one-dimensional integral, can be accurately approximated with various classic numerical integration techniques. When having Gaussian kernels, we use Gauss-Hermite quadrature rule which can provide spectral accuracy (see [1]) in the DGS gradient approximation. Second, the use of Gauss-Hermite quadrature also features embarrassing parallelism as the random sampling used in existing ES methods. Since the communication cost between computing processors/cores is neglectable, the total computing time for each iteration of training does not increase with the number of environment simulations given sufficient computing resources. Third, the directional smoothing approach enables nonlocal exploration which takes into account large variation of the objective function and disregards local fluctuations. This property will greatly help skipping local optima or saddle points during the training. It is demonstrated in §4 that the proposed strategy can significantly reduce the number of iterations in training several benchmark reinforcement learning problems, compared to the state-of-the-art baselines.

The rest of the paper is organized as follows: the RL problem under consideration and a brief review about the classic Gaussian smoothing technique is given in Section 2, the DGS gradient and the corresponding DGS-ES algorithm is introduced in Section 3, extensive numerical experiments including tests on benchmark optimization problems and on several benchmark RL problems are provided in Section 4, some concluding remarks are given in Section 5.

2.   Problem setting

We study the continuous-time, continuous-state stochastic control problem via RL. The evolution of the state ${\boldsymbol{s}}_t \in \mathcal{S}$ (the state space) depends on the control (i.e., the action in RL) ${\boldsymbol{a}}_t \in \mathcal{A}$ (the action space), by a stochastic differential equation, i.e., the state dynamics,

where $b$ is the local drift, $\sigma$ is the local diffusion, and $w_t$ is the standard Brownian motion. Let $\pi({\boldsymbol{a}}_t | {\boldsymbol{s}}_t)$ denote a control strategy (i.e., the policy in RL) that is a conditional probability distribution of the action ${\boldsymbol{a}}_t$ given the current state ${\boldsymbol{s}}_t$. Then the objective of RL is to maximize a functional $\mathcal{J}$ of the policy $\pi$, i.e.,

where $T$ is the terminal time, $r({\boldsymbol{s}}_t, {\boldsymbol{a}}_t)$ is the reward and $0 \le \gamma \le 1$ is a discount rate. After applying a temporal discretization scheme to the state process in (1), the continuous objective functional $\mathcal{J}(\pi)$ can be approximated by a discretized functional,

where $t_n = n \Delta t$ for $n = 0, \ldots, N$.

In policy-based RL, the policy $\pi({\boldsymbol{a}}|{\boldsymbol{s}}; \boldsymbol \theta)$ is parameterized by $\boldsymbol \theta \in \mathbb{R}^d$, where the vector $\boldsymbol \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 $\boldsymbol \theta$ to solve the following optimization problem

where we denote $J(\boldsymbol \theta) : = J(\pi({\boldsymbol{a}} | {\boldsymbol{s}}; \boldsymbol \theta))$ by an abuse of notation. The main challenges for the RL problem under consideration lies in three aspects. First, it is usually true that the local gradient of the environment $\mathcal{S}$ is inaccessible, so automatic differentiation or adjoint methods cannot be used to obtain the local gradient of $J(\boldsymbol \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. Second, the landscape of $J(\boldsymbol \theta)$ are usually highly non-convex and multi-modal, such that the local gradient (if available) may lead to an unsatisfactory local maximum of $J(\boldsymbol \theta)$. Third, the dimension $d$ of the parameters are usually on the order of hundreds or thousands (e.g., when using a neural network to define the policy), which makes many existing optimization methods inapplicable. These challenges motivated us to develop a nonlocal gradient operator and the corresponding approximation for the RL problem under consideration.

2.1.   The standard Gaussian smoothing

We briefly recall the evolution strategy methods, e.g., [13,23], which use a multivariate Gaussian distribution to generate the population around the current parameter value $\boldsymbol \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 [10,20]. Specifically, a smoothed version of $J(\boldsymbol \theta)$ in Eq. (2), denoted by $J_{\boldsymbol \sigma}(\boldsymbol \theta)$, is defined by

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

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

The standard ES method [23] uses Monte Carlo sampling to estimate the gradient $\nabla J_{\boldsymbol \sigma}(\boldsymbol \theta)$ and update the state $\boldsymbol \theta$ from iteration $n$ to $n+1$ by

where $\lambda$ is the learning rate, ${\boldsymbol{u}}_m$ are sampled from the Gaussian distribution $\mathcal{N}(0,{\bf{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 [3]), also [33], for extended discussions on the accuracy of gradient approximations using Eq. (7) and related methods). On the other hand, the evaluations of $J(\boldsymbol \theta_n + \boldsymbol \sigma \circ {\boldsymbol{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. (6), such that the new gradient can be approximated in a much more accurate way and the embarrassing parallelism feature can be retained.

3.   The DGS-ES method for RL

This section introduces our main framework. We start by introducing in Section 3.1 the DGS gradient operator and its approximation using the Gauss-Hermite quadrature rule. In Section 3.2, we describe how to incorporate the DGS gradient operator into the ES for reinforcement learning.

3.1.   The DGS gradient operator

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

where $\boldsymbol \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(\boldsymbol \theta)$ along $\boldsymbol \xi$ in the neighbourhood of $\boldsymbol \theta$. The derivative of $G_{\sigma}(y|\boldsymbol \theta,\boldsymbol \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 \,|\, {\boldsymbol{x}}, \boldsymbol \xi)]$ only involves the directionally smoothed objective function given in Eq. (9).

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

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

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

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

3.2.   The DGS-ES algorithm

Since each component of ${\nabla}_{\boldsymbol \sigma, \boldsymbol \Xi}[J](\boldsymbol \theta)$ in Eq. (11) only involves a one-dimensional integral, we can use Gaussian quadrature rules [22] 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. (10), 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. (13). The weights $\{w_m\}_{m = 1}^M$ and the roots $\{v_m\}_{m = 1}^M$ can be found in [1]. 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}_{\boldsymbol \sigma, \boldsymbol \Xi}[J](\boldsymbol \theta)$ in Eq. (11), 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. (16) indicates that $\widetilde{\nabla}^M_{\boldsymbol \sigma, \boldsymbol \Xi}[J](\boldsymbol \theta)$ is an accurate estimator of the DGS gradient for a small $M$, regardless of the value of $\boldsymbol \sigma$. This enables the use of relatively big values of $\boldsymbol \sigma$ in Eq. (17) to exploit the nonlocal features of the landscape of $J(\boldsymbol \theta)$ in the training process. The nonlocal exploitation ability of $\widetilde{\nabla}^M_{\boldsymbol \sigma, \boldsymbol \Xi}[J]$ is demonstrated in §4 to be effective in reducing the necessary number of iterations to achieve a prescribed reward score. An illustration of the nonlocal exploitation is given in Figure 1 in optimizing the Ackley function (its definition is given in Section 4.1).

On the other hand, as the quadrature weights $w_m$ and $v_m$ defined in Eq. (14) and Eq. (15), are deterministic values, the DGS estimator $\widetilde{\nabla}^M_{\boldsymbol \sigma, \boldsymbol \Xi}[J]({\boldsymbol{x}})$ in Eq. (17) is also a deterministic for fixed $\boldsymbol \Xi$ and $\boldsymbol \sigma$. To introduce random exploration ability to our approach, we add random perturbations to both $\boldsymbol \Xi$ and $\boldsymbol \sigma$. For the orthonormal matrix $\boldsymbol \Xi$, we add a small random rotation, denoted by $\Delta \boldsymbol\Xi$, to the current matrix $\boldsymbol \Xi$. The matrix $\Delta \boldsymbol\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 $\boldsymbol \sigma$ is conducted by drawing random samples from a uniform distribution $U(r-\beta, r+\beta)$ with two hyperparameters $r$ and $\beta$ with $r-\beta>0$. The random perturbation of $\boldsymbol \Xi$ and $\boldsymbol \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 suggest to set the hyperparameters in Algorithm 1 as follows: $M = 7, \alpha = 2.0, r = 1.0, \beta = 0.2, \gamma = 0.01$. Note that both input and output variables need to be properly normalized before running the DGS-ES algorithm. In practice one can tune the critical hyperparameters ($M, r$ and $\ell_r$) given the following suggested range: $M \in [7,9]$, $r \in [0.5,1.0]$ and $\ell_r \in [0.01, 0.1]$. Since our method does not have many hyperparameters, we suggest to use a simple grid search for hyperparameter tuning, and special attention should be spent to tuning the smoothing radius $\sigma$.

4.   Numerical experiments

We evaluate the DGS-ES method using two sets of problems: (a) classic high-dimensional benchmark functions to test the performance of DGS-ES in solving high-dimensional multi-modal optimization problems, and (b) several reinforcement learning benchmark problems.

4.1.   Tests on high-dimensional benchmark functions

Here we investigate the DGS-ES performance on the high-dimensional functions: Sphere, Ackley, Lévy and Rastrigin functions, whose definitions are given below.

● The Sphere function $F_1({\boldsymbol{x}})$ is defined by

where $d$ is the dimension and $\mathit{\boldsymbol{x}} \in [-5.12, 5.12]^d$ is the input domain. The global minimum is $f(\mathit{\boldsymbol{x}}^*) = 0$ at $\mathit{\boldsymbol{x}}^* = (0,...,0)$. It represents convex and isotropic landscapes.

● The Ackley function $F_2({\boldsymbol{x}})$ is defined by

where $d$ is the dimension and $a = 20, b = 0.2, c = 2\pi$ are used in our experiments. The input domain $\mathit{\boldsymbol{x}} \in [-32.768, 32.768]^d$. The global minimum is $f(\mathit{\boldsymbol{x}}^*) = 0$, at $\mathit{\boldsymbol{x}}^* = (0,...,0)$. The Ackley function represents non-convex landscapes with nearly flat outer region. The function poses a risk for optimization algorithms, particularly hill-climbing algorithms, to be trapped in one of its many local minima.

● The Rastrigin function $F_3({\boldsymbol{x}})$ is defined by

where $d$ is the dimension and $\mathit{\boldsymbol{x}} \in [-5.12, 5.12]^d$ is the input domain. The global minimum is $f(\mathit{\boldsymbol{x}}^*) = 0$ at $\mathit{\boldsymbol{x}}^* = (0,...,0)$. This function represents multimodal and separable landscapes.

● The Lévy function $F_4({\boldsymbol{x}})$ is defined by

where $w_i = 1+(x_i-1)/4$ for $i = 1, \ldots, d$ and $\mathit{\boldsymbol{x}} \in [-10, 10]^d$ is the input domain. The global minimum is $f(\mathit{\boldsymbol{x}}^*) = 0$ at $\mathit{\boldsymbol{x}}^* = (0,...,0)$. This represents multi-modal and anisotropic landscapes.

We compare the DGS-ES method with the following blackbox optimization methods: (a) Evolution Strategy (ES) in [23]; (b) Covariance matrix adaptation evolution strategy (CMA-ES) in [12] (we used the implementation of CMA-ES in the pycma open-source code from https://github.com/CMA-ES/pycma) and (c) the BFGS method in the Scipy library.

Figure 2 shows the scalability of the DGS-ES algorithm with the dimension of the objective function. We perform 20 repeated independent trials for each method to show the statistical performance. For the sphere function, the convergence rate does not change when we increase the dimension from 10 to 1000. Such property empirically carries over to the non-convex Lévy and Rastrigin functions. The reason is that both Lévy and Rastrigin have a globally near-convex structure when smoothing out their local minima. In contrast, the Ackley function is highly concave, as a large part of its surface is almost flat regardless of the small local minima. In this case, it takes many iterations for DGS-ES to search for the global minimum hidden in the middle of the flat surface.

Figure 3 shows the comparison between the DGS-ES and the baselines for optimizing the four benchmark functions in Figure 2 in 2000-dimensional space. For the sphere function, BFGS is a clear winner due to its optimal performance in handling convex functions. Among the three ES-type methods, the DGS-ES has much better performance than its competitors. For the other three functions, the DGS-ES method shows significant advantages over other methods.

4.2.   Tests on RL benchmark problems

To evaluate the DGS-ES algorithm, we test its performance on two classes of reinforcement learning environments: three classical control theory problems from OpenAI Gym (https://github.com/openai/gym) [4] and three continuous control tasks simulated using PyBullet (2.6.5) (https://pybullet.org/) [8] which is an open-source library. Within OpenAI Gym, we demonstrate the proposed approach on three benchmark examples: CartPole-v0 (discrete), MountainCarContinuous-v0 (continuous), Pendulum-v0 (continuous). The maximum time steps for these examples are 200,999 and 200, respectively. More details about the environment and reward settings can be found in [4]. We also examine the DGS-ES algorithm on the challenging continuous robotic control problems in PyBullet library, namely HopperBulletEnv-v0, InvertedPendulumBulletEnv-v0 and ReacherBulletEnv-v0. In these three tasks, the maximum time steps are 1000, 1000 and 150, respectively. For the purpose of reproducible comparison, we employ the original environment settings from the OpenAI Gym and the PyBullet library without modifying the rewards or the environments.

For our implementation of DGS-ES, we define our policies as a two-layer feed-forward neural network with 16 hidden nodes and tanh activation functions. For gradient-based optimization, we use Adam to adaptively update the network parameters with a learning rate of $\ell_r = 0.1$. We choose the hyperparameters used in Algorithm 1 as follows: $M = 7, \alpha = 2.0, r = 1.0, \beta = 0.2, \gamma = 0.01$. In practice one can tune the critical hyperparameters ($M, r$ and $\ell_r$) given the following suggested range: $M \in [7,9]$, $r \in [0.5,1.0]$ and $\ell_r \in [0.01, 0.1]$. For each task, our results are performed over 5 repeated independent trials (different random seeds) of the Gym/PyBullet simulators and the network policy initialization.

The DGS-ES algorithm is specifically amenable to parallelization since it only needs to communicate scalars, allowing it to scale to over a large number of parallel workers. We implement a distributed version of Algorithm 1 to the reinforcement learning tasks. The distributed DGS-ES is implemented using PyTorch [21] combined with Ray [19] (https://github.com/ray-project/ray), which does not rely on special networking setup and is tested on large-scale high performance computing facilities with thousands of computing nodes/workers.

Comparison metric. As the motivation of this work is to accelerate time-to-solution of reinforcement training under the assumption that sufficient distributed computing resource is available, we use a different metric to evaluate the performance of DGS-ES and the baselines. Specifically, we are interested in the average return $\mathbb{E}[J]$ versus the number of iterations, i.e., $N$ in Algorithm 1, because those iterations cannot be parallelized.

Baseline methods. We compare 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.

● ES: The Evolution Strategy proposed in [23]. 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 [5]. We used the implementation released by the authors at https://github.com/jparkerholder/ASEBO.

● PPO: Proximal Policy Optimization in [25], which is available in OpenAI's baselines repository at https://github.com/openai/baselines [9].

● TRPO: Trust Region Policy Optimization, developed by [24]. We also used the OpenAI's baselines implementation [9].

● DDPG: Deep Deterministic Policy Gradient, proposed by [14]. 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 [11], 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 hyperparameters for all algorithms above are set to match the original papers without further tuning to improve performance on the testing benchmark examples.

Comparative evaluation. Figure 4 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(\boldsymbol \theta)$ is very flat at most locations in the parameter space, which makes it difficult to capture the peak of $J(\boldsymbol \theta)$. This test is a good demonstration of the nonlocal exploration ability of the DGS-ES method. Since the mean of the smoothing factor $\boldsymbol \sigma$ is set to 1.0, DGS-ES can capture the peak of $J(\boldsymbol \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.

Figure 5 shows the comparison results of Hopper-v0, InvertedPendulum-v0 and Reacher-v0 problems from the PyBullet library. We compare 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 [9] and [11]. 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 [11], 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.

Figure 6 illustrates the effect of the radius $\boldsymbol \sigma$ of $\widetilde{\nabla}^M_{\boldsymbol \sigma, \boldsymbol \Xi}[J](\boldsymbol \theta)$ on the performance of the DGS-ES method in solving the Reacher-v0 problem. All the simulations are done using the same initialization. We set the mean of, i.e., the hyperparameter $r$ in Algorithm 1, to 0.5, 0.05, 0.01, and 0.005. It is easy to see that the performance of DGS-ES deteriorates with the decrease of $\boldsymbol \sigma$. Since it is evident that the surface of $J(\boldsymbol \theta)$ is not convex and may have many local maxima, a relatively big radius $\boldsymbol \sigma$ is necessary to help DGS-ES skip the local maxima. As $\boldsymbol \sigma$ becomes smaller, the DGS gradient $\widetilde{\nabla}^M_{\boldsymbol \sigma, \boldsymbol \Xi}[J](\boldsymbol \theta)$ converges to the local gradient $\nabla J(\boldsymbol \theta)$, which may get the optimizer trapped in a local maximum.

4.3.   Aerospace hierarchical stiffened shell design

Now we demonstrate the DGS-ES method on a real-world stiffened shell design problem. Due to its high strength and stiffness, the hierarchical stiffened shell has been widely used in aerospace engineering [31,32]. However, it is challenging to fully explore its optimal buckling load-carrying capacity. The goal of this work is to improve the load-carrying capacity by optimizing the representative unit cell of hierarchical stiffened shell where the inputs are 9 size variables (widths, heights and thickness) for major and minor stiffeners, as shown in Figure 7(a), and the output is the carrying-load capability factor which is calculated by high-fidelity numerical simulation, e.g., finite element method (FEM), see Figure 7 (b). Although the high-fidelity FEM simulation is time-consuming, many open-source codes and commercial software have improved the computational efficiency via scalable parallelism and GPU acceleration. It is thus feasible to apply the DGS-ES method by implementing parallel FEM simulations in supercomputers. Figure 7 presents the comparison between DGS-ES and the other blackbox optimization methods. We observe that DGS-ES outperforms all other algorithms and achieves faster convergence using only 100 iterations. The robustness of DGS-ES also outperforms the alternatives.

5.   Conclusion

Despite the successful demonstration shown in §4, there are several limitations with the DGS-ES method for reinforcement learning. First, it requires a powerful enough cluster or distributed computing resources to show superior performance. Even though more and more parallel environments for complex RL tasks, those parallel codes might not compatible with a cluster/supercomputer with a specific architecture. This will need some extra effort to modify environment codes, in order to exploit the advantage of the DGS-ES approach. Second, asynchronization between different environment simulations may drag down the total performance. In a distributed computing system, all the parallel workers receive the same number of environment simulations. However, as different parameter values may lead to different termination times of the environment simulations, there will be a potential waste of computing resources due to such asynchronization. Thus, a better scheduling algorithm is needed to further improve the performance of DGS-ES in RL. Third, even though the performance of the DGS-ES is not very sensitive to the hyperparameters, especially the radius $\boldsymbol \sigma$ in the experiments conducted in this work, optimal or even viable hyperparameters of the DGS-ES method are still problem dependent, which means hyperparameter tuning may be needed when applying the method to another RL problem.

Acknowledgments

This material was based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under contract and award numbers ERKJ369, by the DOE SciDac FastMath project, and by the Artificial Intelligence Initiative at the Oak Ridge National Laboratory (ORNL). ORNL is operated by UT-Battelle, LLC., for the U.S. Department of Energy under Contract DE-AC05-00OR22725.