Title: SAMPLE EFFICIENT POLICY GRADIENT METHODS WITH RECURSIVE VARIANCE REDUCTION

Abstract: Improving the sample efficiency in reinforcement learning has been a longstanding research problem. In this work, we aim to reduce the sample complexity of existing policy gradient methods. We propose a novel policy gradient algorithm called SRVR-PG, which only requires O(1/ 3/2 ) 1 episodes to find anapproximate stationary point of the nonconcave performance function J(θ) (i.e., θ such that ∇J(θ) 2 2 ≤ ). This sample complexity improves the existing result O(1/ 5/3 ) for stochastic variance reduced policy gradient algorithms by a factor of O(1/ 1/6 ). In addition, we also propose a variant of SRVR-PG with parameter exploration, which explores the initial policy parameter from a prior probability distribution. We conduct numerical experiments on classic control problems in reinforcement learning to validate the performance of our proposed algorithms.

1. INTRODUCTION:

Reinforcement learning (RL) (Sutton & Barto, 2018) has received significant success in solving various complex problems such as learning robotic motion skills (Levine et al., 2015) , autonomous driving (Shalev-Shwartz et al., 2016) and Go game (Silver et al., 2017) , where the agent progressively interacts with the environment in order to learn a good policy to solve the task. In RL, the agent makes its decision by choosing the action based on the current state and the historical rewards it has received so far. After performing the chosen action, the agent's state will change according to some transition probability model and a new reward would be revealed to the agent by the environment based on the action and new state. Then the agent continues to choose the next action until it reaches a terminal state. The aim of the agent is to maximize its expected cumulative rewards. Therefore, the pivotal problem in RL is to find a good policy which is a function that maps the state space to the action space and thus informs the agent which action to take at each state. To optimize the agent's policy in the high dimensional continuous action space, the most popular approach is the policy gradient method (Sutton et al., 2000) that parameterizes the policy by an unknown parameter θ ∈ R d and directly optimizes the policy by finding the optimal θ. The objective function J(θ) is chosen to be the performance function, which is the expected return under a specific policy and is usually non-concave. Our goal is to maximize the value of J(θ) by finding a stationary point θ * such that ∇J(θ * ) 2 = 0 using gradient based algorithms.

Due to the expectation in the definition of J(θ), it is usually infeasible to compute the gradient exactly. In practice, one often uses stochastic gradient estimators such as REINFORCE (Williams, 1992) , PGT (Sutton et al., 2000) and GPOMDP (Baxter & Bartlett, 2001) to approximate the gradient of the expected return based on a batch of sampled trajectories. However, this approximation will introduce additional variance and slow down the convergence of policy gradient, which thus requires a huge amount of trajectories to find a good policy. Theoretically, these stochastic gradient (SG) based algorithms require O(1/ 2 ) trajectories (Robbins & Monro, 1951) to find an -approximate stationary point such that E[ ∇J(θ) 2 2 ] ≤ . In order to reduce the variance of policy gradient algorithms, Papini et al. (2018) proposed a stochastic variance-reduced policy gradient (SVRPG) O(1/ 2 ) PGT (Sutton et al., 2000) O(1/ 2 ) GPOMDP (Baxter & Bartlett, 2001) O(1/ 2 ) SVRPG (Papini et al., 2018) O(1/ 2 ) SVRPG (Xu et al., 2019) O(1/ 5/3 ) SRVR-PG (This paper)

O(1/ 3/2 ) algorithm by borrowing the idea from the stochastic variance reduced gradient (SVRG) (Johnson & Zhang, 2013; Allen-Zhu & Hazan, 2016; Reddi et al., 2016a) in stochastic optimization. The key idea is to use a so-called semi-stochastic gradient to replace the stochastic gradient used in SG methods. The semi-stochastic gradient combines the stochastic gradient in the current iterate with a snapshot of stochastic gradient stored in an early iterate which is called a reference iterate. In practice, SVRPG saves computation on trajectories and improves the performance of SG based policy gradient methods. Papini et al. (2018) also proved that SVRPG converges to an -approximate stationary point θ of the nonconcave performance function J(θ) with E[ ∇J(θ) 2 2 ] ≤ after O(1/ 2 ) trajectories, which seems to have the same sample complexity as SG based methods. Recently, the sample complexity of SVRPG has been improved to O(1/ 5/3 ) by a refined analysis (Xu et al., 2019) , which theoretically justifies the advantage of SVRPG over SG based methods.

This paper continues on this line of research. We propose a Stochastic Recursive Variance Reduced Policy Gradient algorithm (SRVR-PG), which provably improves the sample complexity of SVRPG. At the core of our proposed algorithm is a recursive semi-stochastic policy gradient inspired from the stochastic path-integrated differential estimator (Fang et al., 2018) , which accumulates all the stochastic gradients from different iterates to reduce the variance. We prove that SRVR-PG only takes O(1/ 3/2 ) trajectories to converge to an -approximate stationary point θ of the performance function, i.e., E[ ∇J(θ) 2 2 ] ≤ . We summarize the comparison of SRVR-PG with existing policy gradient methods in terms of sample complexity in Table 1 . Evidently, the sample complexity of SRVR-PG is lower than that of REINFORCE, PGT and GPOMDP by a factor of O(1/ 1/2 ), and is lower than that of SVRPG (Xu et al., 2019) by a factor of O(1/ 1/6 ).

In addition, we integrate our algorithm with parameter-based exploration (PGPE) method (Sehnke et al., 2008; 2010) , and propose a SRVR-PG-PE algorithm which directly optimizes the prior probability distribution of the policy parameter θ instead of finding the best value. The proposed SRVR-PG-PE enjoys the same trajectory complexity as SRVR-PG and performs even better in some applications due to its additional exploration over the parameter space. Our experimental results on classical control tasks in reinforcement learning demonstrate the superior performance of the proposed SRVR-PG and SRVR-PG-PE algorithms and verify our theoretical analysis.

1.1. ADDITIONAL RELATED WORK:

We briefly review additional relevant work to ours with a focus on policy gradient based methods. For other RL methods such as value based (Watkins & Dayan, 1992; Mnih et al., 2015) and actorcritic (Konda & Tsitsiklis, 2000; Peters & Schaal, 2008a; Silver et al., 2014 ) methods, we refer the reader to Peters & Schaal (2008b) ; Kober et al. (2013) ; Sutton & Barto (2018) for a complete review.

To reduce the variance of policy gradient methods, early works have introduced unbiased baseline functions (Baxter & Bartlett, 2001; Greensmith et al., 2004; Peters & Schaal, 2008b) to reduce the variance, which can be constant, time-dependent or state-dependent. Schulman et al. (2015b) proposed the generalized advantage estimation (GAE) to explore the trade-off between bias and variance of policy gradient. Recently, action-dependent baselines are also used in Tucker et al. (2018) ; Wu et al. (2018) which introduces bias but reduces variance at the same time. Sehnke et al. (2008; 2010) proposed policy gradient with parameter-based exploration (PGPE) that explores in the parameter space. It has been shown that PGPE enjoys a much smaller variance (Zhao et al., 2011) . The Stein variational policy gradient method is proposed in Liu et al. (2017) . See Peters & Schaal (2008b); Deisenroth et al. (2013) ; Li (2017) for a more detailed survey on policy gradient.

Stochastic variance reduced gradient techniques such as SVRG (Johnson & Zhang, 2013; Xiao & Zhang, 2014) , batching SVRG (Harikandeh et al., 2015) , SAGA (Defazio et al., 2014) and SARAH (Nguyen et al., 2017) were first developed in stochastic convex optimization. When the objective function is nonconvex (or nonconcave for maximization problems), nonconvex SVRG (Allen-Zhu & Hazan, 2016; Reddi et al., 2016a) and SCSG (Lei et al., 2017; Li & Li, 2018) were proposed and proved to converge to a first-order stationary point faster than vanilla SGD (Robbins & Monro, 1951) with no variance reduction. The state-of-the-art stochastic variance reduced gradient methods for nonconvex functions are the SNVRG (Zhou et al., 2018) and SPIDER (Fang et al., 2018) algorithms, which have been proved to achieve near optimal convergence rate for smooth functions.

There are yet not many papers studying variance reduced gradient techniques in RL. Du et al. (2017) first applied SVRG in policy evaluation for a fixed policy. Xu et al. (2017) introduced SVRG into trust region policy optimization for model-free policy gradient and showed that the resulting algorithm SVRPO is more sample efficient than TRPO. Yuan et al. (2019) further applied the techniques in SARAH (Nguyen et al., 2017) and SPIDER (Fang et al., 2018) to TRPO (Schulman et al., 2015a) . However, no analysis on sample complexity (i.e., number of trajectories required) was provided in the aforementioned papers (Xu et al., 2017; Yuan et al., 2019) . We note that a recent work by Shen et al. (2019) proposed a Hessian aided policy gradient (HAPG) algorithm that converges to the stationary point of the performance function within O(H 2 / 3/2 ) trajectories, which is worse than our result by a factor of O(H 2 ) where H is the horizon length of the environment. Moreover, they need additional samples to approximate the Hessian vector product, and cannot handle the policy in a constrained parameter space. Another related work pointed out by the anonymous reviewer is Yang & Zhang (2019) , which extended the stochastic mirror descent algorithm (Ghadimi et al., 2016) in the optimization field to policy gradient methods and achieved O(H 2 / 2 ) sample complexity. After the ICLR conference submission deadline, Yang & Zhang (2019) revised their paper by adding a new variance reduction algorithm that achieves O(H 2 / 3/2 ) sample complexity, which is also worse than our result by a factor of O(H 2 ).

Apart from the convergence analysis of the general nonconcave performance functions, there has emerged a line of work (Cai et al., 2019; Liu et al., 2019; Yang et al., 2019; Wang et al., 2019) that studies the global convergence of (proximal/trust-region) policy optimization with neural network function approximation, which applies the theory of overparameterized neural networks (Du et al., 2019b; a; Allen-Zhu et al., 2019; Zou et al., 2019; Cao & Gu, 2019) to reinforcement learning.

Notation v 2 denotes the Euclidean norm of a vector v ∈ R d and A 2 denotes the spectral norm of a matrix A ∈ R d×d . We write a n = O(b n ) if a n ≤ Cb n for some constant C > 0. The Dirac delta function δ(x) satisfies δ(0) = +∞ and δ(x) = 0 if x = 0. Note that δ(x) satisfies +∞ -∞ δ(x)dx = 1. For any α > 0, we define the Rényi divergence (Rényi et al., 1961) between distributions P and Q as

D α (P ||Q) = 1 α -1 log 2 x P (x) P (x) Q(x) α-1 dx,

which is non-negative for all α > 0. The exponentiated Rényi divergence is d α (P ||Q) = 2 Dα(P ||Q) .

2. BACKGROUNDS ON POLICY GRADIENT:

Markov Decision Process: A discrete-time Markov Decision Process (MDP) is a tuple M = {S, A, P, r, γ, ρ}. S and A are the state and action spaces respectively. P(s |s, a) is the transition probability of transiting to state s after taking action a at state s. Function r(s, a) : S × A → [-R, R] emits a bounded reward after the agent takes action a at state s, where R > 0 is a constant. γ ∈ (0, 1) is the discount factor. ρ is the distribution of the starting state. A policy at state s is a probability function π(a|s) over action space A. In episodic tasks, following any stationary policy, the agent can observe and collect a sequence of state-action pairs τ = {s 0 , a 0 , s 1 , a 1 , . . . , s H-1 , a H-1 , s H }, which is called a trajectory or episode. H is called the trajectory horizon or episode length. In practice, we can set H to be the maximum value among all the actual trajectory horizons we have collected. The sample return over one trajectory τ is defined as the discounted cumulative reward R(τ ) =

H-1 h=0 γ h r(s h , a h ). Policy Gradient: Suppose the policy, denoted by π θ , is parameterized by an unknown parameter θ ∈ R d . We denote the trajectory distribution induced by policy π θ as p(τ |θ). Then

p(τ |θ) = ρ(s 0 ) H-1 h=0 π θ (a h |s h )P (s h+1 |s h , a h ).

(2.1)

We define the expected return under policy π θ as

J(θ) = E τ ∼p(•|θ) [R(τ )|M],

which is also called the performance function. To maximize the performance function, we can update the policy parameter θ by iteratively running gradient ascent based algorithms, i.e., θ k+1 = θ k + η∇ θ J(θ k ), where η > 0 is the step size and the gradient ∇ θ J(θ) is derived as follows:

∇ θ J(θ) = τ R(τ )∇ θ p(τ |θ)dτ = τ R(τ )(∇ θ p(τ |θ)/p(τ |θ))p(τ |θ)dτ = E τ ∼p(•|θ) [∇ θ log p(τ |θ)R(τ )|M]. (2.2)

However, it is intractable to calculate the exact gradient in (2.2) since the trajectory distribution p(τ |θ) is unknown. In practice, policy gradient algorithm samples a batch of trajectories {τ i } N i=1 to approximate the exact gradient based on the sample average over all sampled trajectories:

∇ θ J(θ) = 1 N N i=1 ∇ θ log p(τ i |θ)R(τ i ).

(2.3)

At the k-th iteration, the policy is then updated by θ k+1 = θ k + η ∇ θ J(θ k ). According to (2.1), we know that ∇ θ log p(τ i |θ) is independent of the transition probability matrix P . Recall the definition of R(τ ), we can rewrite the approximate gradient as follows

EQUATION

where τ i = {s i 0 , a i 0 , s i 1 , a i 1 , . . . , s i H-1 , a i H-1 , s i H } for all i = 1, . . . , N and g(τ i |θ) is an unbiased gradient estimator computed based on the i-th trajectory τ i . The gradient estimator in (2.4) is based on the likelihood ratio methods and is often referred to as the REINFORCE gradient estimator (Williams, 1992) . Since E[∇ θ log π θ (a|s)] = 0, we can add any constant baseline b t to the reward that is independent of the current action and the gradient estimator still remains unbiased. With the observation that future actions do not depend on past rewards, another famous policy gradient theorem (PGT) estimator (Sutton et al., 2000) removes the rewards from previous states:

g(τ i |θ) = H-1 h=0 ∇ θ log π θ (a i h |s i h ) H-1 t=h γ t r(s i t , a i t ) -b t , (2.5)

where b t is a constant baseline. It has been shown (Peters & Schaal, 2008b ) that the PGT estimator is equivalent to the commonly used GPOMDP estimator (Baxter & Bartlett, 2001 ) defined as follows:

g(τ i |θ) = H-1 h=0 h t=0 ∇ θ log π θ (a i t |s i t ) γ h r(s i h , a i h ) -b h . (2.6)

All the three gradient estimators mentioned above are unbiased (Peters & Schaal, 2008b) . It has been proved that the variance of the PGT/GPOMDP estimator is independent of horizon H while the variance of REINFORCE depends on H polynomially (Zhao et al., 2011; Pirotta et al., 2013) . Therefore, we will focus on the PGT/GPOMDP estimator in this paper and refer to them interchangeably due to their equivalence.

3. THE PROPOSED ALGORITHM:

The approximation in (2.3) using a batch of trajectories often causes a high variance in practice.

In this section, we propose a novel variance reduced policy gradient algorithm called stochastic recursive variance reduced policy gradient (SRVR-PG), which is displayed in Algorithm 1. Our SRVR-PG algorithm consists of S epochs. In the initialization, we set the parameter of a reference policy to be θ 0 = θ 0 . At the beginning of the s-th epoch, where s = 0, . . . , S -1, we set the initial policy parameter θ s+1 0 to be the same as that of the reference policy θ s . The algorithm then samples N episodes {τ i } N i=1 from the reference policy π θ s to compute a gradient estimator

v s 0 = 1/N N i=1 g(τ i | θ s ), where g(τ i | θ s )

is the PGT/GPOMDP estimator. Then the policy is immediately update as in Line 6 of Algorithm 1.

Within the epoch, at the t-th iteration, SRVR-PG samples B episodes {τ j } B j=1 based on the current policy π θ s+1 t . We define the following recursive semi-stochastic gradient estimator:

v s+1 t = 1 B B j=1 g(τ j |θ s+1 t ) - 1 B B j=1 g ω (τ j |θ s+1 t-1 ) + v s+1 t-1 , (3.1)

where the first term is a stochastic gradient based on B episodes sampled from the current policy, and the second term is a stochastic gradient defined based on the step-wise important weight between the current policy π θ s+1 t and the reference policy π θ s . Take the GPOMDP estimator for example, for a behavior policy π θ1 and a target policy π θ2 , the step-wise importance weighted estimator is defined as follows

EQUATION

where

ω 0:h (τ |θ 2 , θ 1 ) = h h =0 π θ2 (a h |s h )/π θ1 (a h |s h ) is the importance weight from p(τ h |θ s+1 t ) to p(τ h |θ s+1 t-1

) and τ h is a truncated trajectory {(a t , s t )} h t=0 from the full trajectory τ . It is easy to

verify that E τ ∼p(τ |θ1) [g ω (τ j |θ 1 )] = E τ ∼p(τ |θ2) [g(τ |θ 2 )].

The difference between the last two terms in (3.1) can be viewed as a control variate to reduce the variance of the stochastic gradient. In many practical applications, the policy parameter space is a subset of R d , i.e., θ ∈ Θ with Θ ⊆ R d being a convex set. In this case, we need to project the updated policy parameter onto the constraint set. Base on the semi-stochastic gradient (3.1), we can update the policy parameter using projected gradient ascent along the direction of v s+1 t :

θ s+1 t+1 = P Θ (θ s+1 t + ηv s+1 t

), where η > 0 is the step size and the projection operator associated with Θ is defined as

P Θ (θ) = argmin u∈Θ θ -u 2 2 = argmin u∈R d 1 Θ (u) + 1 2η θ -u 2 2 , (3.3) where 1 Θ (u) is the set indicator function on Θ, i.e., 1 Θ (u) = 0 if u ∈ Θ and 1 Θ (u) = +∞ otherwise.

η > 0 is any finite real value and is chosen as the step size in our paper. It is easy to see that 1 Θ (•) is nonsmooth. At the end of the s-th epoch, we update the reference policy as θ s+1 = θ s+1 m , where θ s+1 m is the last iterate of this epoch. The goal of our algorithm is to find a point θ ∈ Θ that maximizes the performance function J(θ) subject to the constraint, namely, max θ∈Θ J(θ) = max θ∈R d {J(θ) -1 Θ (θ)}. The gradient norm ∇J(θ) 2 is not sufficient to characterize the convergence of the algorithm due to additional the constraint. Following the literature on nonsmooth optimization (Reddi et al., 2016b; Ghadimi et al., 2016; Nguyen et al., 2017; Li & Li, 2018; Wang et al., 2018) , we use the generalized first-order stationary condition: G η (θ) = 0, where the gradient mapping G η is defined as follows

G η (θ) = 1 η (P Θ (θ + η∇J(θ)) -θ). (3.4)

We can view G η as a generalized projected gradient at θ. By definition if Θ = R d , we have G η (θ) ≡ ∇J(θ). Therefore, the policy is update is displayed in Line 10 in Algorithm 1, where

Algorithm 1 Stochastic Recursive Variance Reduced Policy Gradient (SRVR-PG) 1: Input: number of epochs S, epoch size m, step size η, batch size N , mini-batch size B, gradient estimator g, initial parameter θ 0 = θ 0 ∈ Θ 2: for s = 0, . . . , S -1 do 3:

θ s+1 0 = θ s 4: Sample N trajectories {τ i } from p(•| θ s ) 5: v s+1 0 = ∇ θ J( θ s ) := 1/N N i=1 g(τ i | θ s ) 6: θ s+1 1 = P Θ (θ s+1 0 + ηv s+1 0 ) 7: for t = 1, . . . , m -1 do 8: Sample B trajectories {τ j } from p(•|θ s+1 t ) 9: v s+1 t = v s+1 t-1 + 1 B B j=1 g τ j |θ s+1 t -g ω τ j |θ s+1 t-1

10:

θ s+1 t+1 = P Θ (θ s+1 t + ηv s+1 t ) 11:

end for 12: 2019) due to the gradient estimator g ω (τ j |θ s+1 t-1 ) defined in (3.2) that is equipped with step-wise importance weights. This term is essential to deal with the non-stationarity of the distribution of the trajectory τ . Specifically, {τ j } B j=1 are sampled from policy π θ s+1 t while the PGT/GPOMDP estimator g(•|θ s+1 t-1 ) is defined based on policy π θ s+1 t-1 according to (2.6). This inconsistency introduces extra challenges in the convergence analysis of SRVR-PG. Using importance weighting, we can obtain

θ s+1 =

E τ ∼p(τ |θ s+1 t ) [g ω (τ |θ s+1 t-1 )] = E τ ∼p(τ |θ s+1 t-1 ) [g(τ |θ s+1 t-1

)], which eliminates the inconsistency caused by the varying trajectory distribution.

It is worth noting that the semi-stochastic gradient in (3.1) also differs from the one used in SVRPG (Papini et al., 2018) because we recursively update v s+1 t using v s+1 t-1 from the previous iteration, while SVRPG uses a reference gradient that is only updated at the beginning of each epoch. Moreover, SVRPG wastes N trajectories without updating the policy at the beginning of each epoch, while Algorithm 1 updates the policy immediately after this sampling process (Line 6), which saves computation in practice.

We notice that very recently another algorithm called SARAPO (Yuan et al., 2019) is proposed which also uses a recursive gradient update in trust region policy optimization (Schulman et al., 2015a) . Our Algorithm 1 differs from their algorithm at least in the following ways: (1) our recursive gradient v s t defined in (3.1) has an importance weight from the snapshot gradient while SARAPO does not;

(2) we are optimizing the expected return while Yuan et al. (2019) optimizes the total advantage over state visitation distribution and actions under KullbackLeibler divergence constraint; and most importantly (3) there is no convergence or sample complexity analysis for SARAPO.

4. MAIN THEORY:

In this section, we present the theoretical analysis of Algorithm 1. We first introduce some common assumptions used in the convergence analysis of policy gradient methods. Assumption 4.1. Let π θ (a|s) be the policy parameterized by θ. There exist constants G, M > 0 such that the gradient and Hessian matrix of log π θ (a|s) with respect to θ satisfy

∇ θ log π θ (a|s) ≤ G, ∇ 2 θ log π θ (a|s) 2 ≤ M, for all a ∈ A and s ∈ S.

The above boundedness assumption is reasonable since we usually require the policy function to be twice differentiable and easy to optimize in practice. Similarly, in Papini et al. (2018) , the authors assume that ∂ ∂θi log π θ (a|s) and ∂ 2 ∂θi∂θj log π θ (a|s) are upper bounded elementwisely, which is actually stronger than our Assumption 4.1.

In the following proposition, we show that Assumption4.1 directly implies that the Hessian matrix of the performance function ∇ 2 J(θ) is bounded, which is often referred to as the smoothness assumption and is crucial in analyzing the convergence of nonconvex optimization (Reddi et al., 2016a; Allen-Zhu & Hazan, 2016) . Proposition 4.2. Let g(τ |θ) be the PGT estimator defined in (2.5). Assumption 4.1 implies:

(1). g(τ |θ 1 ) -g(τ |θ 2 ) 2 ≤ L θ 1 -θ 2 2 , ∀θ 1 , θ 2 ∈ R d , where the smoothness parameter is L = M R/(1 -γ) 2 + 2G 2 R/(1 -γ) 3 ; (2). J(θ) is L-smooth, namely ∇ 2 θ J(θ) 2 ≤ L; (3). g(τ |θ) 2 ≤ C g for all θ ∈ R d , with C g = GR/(1 -γ) 2 .

Similar properties are also proved in Xu et al. (2019) . However, in contrast to their results, the smoothness parameter L and the bound on the gradient norm here do not rely on horizon H. When H ≈ 1/(1-γ) and γ is sufficiently close to 1, we can see that the order of the smoothness parameter is O(1/(1 -γ) 3 ), which matches the order O(H 2 /(1 -γ)) in Xu et al. (2019) . The next assumption requires the variance of the gradient estimator is bounded. Assumption 4.3. There exists a constant ξ > 0 such that Var g(τ |θ) ≤ ξ 2 , for all policy π θ .

In Algorithm 1, we have used importance sampling to connect the trajectories between two different iterations. The following assumption ensures that the variance of the importance weight is bounded, which is also made in Papini et al. (2018) ; Xu et al. (2019) .

Assumption 4.4. Let ω(•|θ 1 , θ 2 ) = p(•|θ 1 )/p(•|θ 2 ).

There is a constant W < ∞ such that for each policy pairs encountered in Algorithm 1,

Var(ω(τ |θ 1 , θ 2 )) ≤ W, ∀θ 1 , θ 2 ∈ R d , τ ∼ p(•|θ 2 ).

4.1. CONVERGENCE RATE AND SAMPLE COMPLEXITY OF SRVR-PG:

Now we are ready to present the convergence result of SRVR-PG to a stationary point: Theorem 4.5. Suppose that Assumptions 4.1, 4.3 and 4.4 hold. In Algorithm 1, we choose the step size η ≤ 1/(4L) and epoch size m and mini-batch size B such that

B ≥ 72mηG 2 (2G 2 /M + 1)(W + 1)γ (1 -γ) 2 .

Then the generalized projected gradient of the output of Algorithm 1 satisfies

E G η θ out 2 2 ≤ 8[J(θ * ) -J(θ 0 ) -1 Θ (θ * ) + 1 Θ (θ 0 )] ηSm + 6ξ 2 N ,

where θ * = argmax θ∈Θ J(θ). 2019) are for ∇ θ J(θ)2 2 ≤ , while our result in Corollary 4.7 is more general. In particular, when the policy parameter θ is defined on the whole space R d instead of Θ, our result reduces to the case for Xu et al. (2019) , the authors improved the sample complexity of SVRPG (Papini et al., 2018) from O(1/ 2 ) to O(1/ 5/3 ) by a sharper analysis. According to Corollary 4.7, SRVR-PG only needs O(1/ 3/2 ) number of trajectories to achieve ∇ θ J(θ) 2 2 ≤ , which is lower than the sample complexity of SVRPG by a factor of O(1/ 1/6 ). This improvement is more pronounced when the required precision is very small.

∇ θ J(θ) 2 2 ≤ since Θ = R d and G η (θ) = ∇ θ J(θ). In

4.2. IMPLICATION FOR GAUSSIAN POLICY:

Now, we consider the Gaussian policy model and present the sample complexity of SRVR-PG in this setting. For bounded action space A ⊂ R, a Gaussian policy parameterized by θ is defined as

π θ (a|s) = 1 √ 2π exp - (θ φ(s) -a) 2 2σ 2 , (4.1)

where σ 2 is a fixed standard deviation parameter and φ : S → R d is a mapping from the state space to the feature space. For Gaussian policy, under the mild condition that the actions and the state feature vectors are bounded, we can verify that Assumptions 4.1 and 4.3 hold, which can be found in Appendix D. It is worth noting that Assumption 4.4 does not hold trivially for all Gaussian distributions. In particular, Cortes et al. (2010) showed that for two Gaussian distributions π θ1 (a|s) ∼ N (µ 1 , σ 2 1 ) and π θ2 (a|s) ∼ N (µ 2 , σ 2 2 ), if σ 2 > √ 2/2σ 1 , then the variance of ω(τ |θ 1 , θ 2 ) is bounded. For our Gaussian policy defined in (4.1) where the standard deviation σ 2 is fixed, we have σ > √ 2/2σ trivially hold, and therefore Assumption 4.4 holds for some finite constant W > 0 according to (2.1).

Recall that Theorem 4.5 holds for any general models under Assumptions 4.1, 4.3 and 4.4. Based on the above arguments, we know that the convergence analysis in Theorem 4.5 applies to Gaussian policy. In the following corollary, we present the sample complexity of Algorithm 1 for Gaussian policy with detailed dependency on precision parameter , horizon size H and the discount factor γ. Corollary 4.8. Given the Gaussian policy defined in (4.1), suppose Assumption 4.4 holds and we have |a| ≤ C a for all a ∈ A and φ(s) 2 ≤ M φ for all s ∈ S, where C a , M φ > 0 are constants. If we set step size as η = O((1-γ) 3 ), the mini-batch sizes and epoch length as the experiments, we use the Gaussian policy defined in (4.1). In addition, we found that the proposed algorithm works well without the extra projection step. Therefore, we did not use projection in our experiments. For baselines, we compare the proposed SRVR-PG algorithm with the most relevant methods: GPOMDP (Baxter & Bartlett, 2001) and SVRPG (Papini et al., 2018) . For the learning rates η in all of our experiments, we use grid search to directly tune η. For instance, we searched η for the Cartpole problem by evenly dividing the interval [10 -5 , 10 -1 ] into 20 points in the logspace. For the batch size parameters N and B and the epoch length m, according to Corollary 4.7, we choose N = O(1/ ), B = O(1/ 1/2 ) and thus m = O(1/ 1/2 ), where > 0 is a user-defined precision parameter. In our experiments, we set N = C 0 / , B = C 1 / 1/2 and m = C 2 / 1/2 and tune the constant parameters C 0 , C 1 , C 2 using grid search. The detailed parameters used in the experiments are presented in Appendix E.

N = O((1-γ) -3 -1 ), B = O((1 -γ) -1 -1/2 ) and m = O((1 -γ) -2 -1/2 ), then the output of Algorithm 1 satisfies E[ G η (θ out ) 2

We evaluate the performance of different algorithms in terms of the total number of trajectories they require to achieve a certain threshold of cumulative rewards. We run each experiment repeatedly for 10 times and plot the averaged returns with standard deviation. For a given environment, all experiments are initialized from the same random initialization. Figures 1(a) , 1(b) and 1(c) show the results on the comparison of GPOMDP, SVRPG, and our proposed SRVR-PG algorithm across three different RL environments. It is evident that, for all environments, GPOMDP is overshadowed by the variance reduced algorithms SVRPG and SRVR-PG significantly. Furthermore, SRVR-PG outperforms SVRPG in all experiments, which is consistent with the comparison on the sample complexity of GPOMDP, SVRPG and SRVR-PG in Table 1 .

Corollaries 4.7 and 4.8 suggest that when the mini-batch size B is in the order of O( √ N ), SRVR-PG achieves the best performance. Here N is the number of episodes sampled in the outer loop of Algorithm 1 and B is the number of episodes sampled at each inner loop iteration. To validate our theoretical result, we conduct a sensitivity study to demonstrate the effectiveness of different batch sizes within each epoch of SRVR-PG on its performance. The results on different environments are displayed in Figures 1(d) , 1(e) and 1(f) respectively. To interpret these results, we take the Pendulum problem as an example. In this setting, we choose outer loop batch size N of Algorithm 1 to be N = 250. By Corollary 4.8, the optimal choice of batch size in the inner loop of Algorithm 1 is B = C √ N , where C > 1 is a constant depending on horizon H and discount factor γ. shows that B = 50 ≈ 3 √ N yields the best convergence results for SRVR-PG on Pendulum, which validates our theoretical analysis and implies that a larger batch size B does not necessarily result in an improvement in sample complexity, as each update requires more trajectories, but a smaller batch size B pushes SRVR-PG to behave more similar to GPOMDP. Moreover, by comparing with the outer loop batch size N presented in Table 2 for SRVR-PG in Cartpole and Mountain Car environments, we found that the results in Figures 1(d ) and 1(e) are again in alignment with our theory. Due to the space limit, additional experiment results are included in Appendix E.

6. CONCLUSIONS:

We propose a novel policy gradient method called SRVR-PG, which is built on a recursively updated stochastic policy gradient estimator. We prove that the sample complexity of SRVR-PG is lower than the sample complexity of the state-of-the-art SVRPG (Papini et al., 2018; Xu et al., 2019) algorithm. We also extend the new variance reduction technique to policy gradient with parameter-based exploration and propose the SRVR-PG-PE algorithm, which outperforms the original PGPE algorithm both in theory and practice. Experiments on the classic reinforcement learning benchmarks validate the advantage of our proposed algorithms.

A EXTENSION TO PARAMETER-BASED EXPLORATION:

Although SRVR-PG is proposed for action-based policy gradient, it can be easily extended to the policy gradient algorithm with parameter-based exploration (PGPE) (Sehnke et al., 2008) . Unlike action-based policy gradient in previous sections, PGPE does not directly optimize the policy parameter θ but instead assumes that it follows a prior distribution with hyper-parameter ρ: θ ∼ p(θ|ρ). The expected return under the policy induced by the hyper-parameter ρ is formulated as follows2 = p(θ|ρ)p(τ |θ)R(τ )dτ dθ. (A.1) PGPE aims to find the hyper-parameter ρ * that maximizes the performance function J(ρ). Since p(θ|ρ) is stochastic and can provide sufficient exploration, we can choose π θ (a|s) = δ(a -µ θ (s)) to be a deterministic policy, where δ is the Dirac delta function and µ θ (•) is a deterministic function.

For instance, a linear deterministic policy is defined as π θ (a|s) = δ(a -θ s) (Zhao et al., 2011; Metelli et al., 2018) . Given the policy parameter θ, a trajectory τ is only decided by the initial state distribution and the transition probability. Therefore, PGPE is called a parameter-based exploration approach. Similar to the action-based policy gradient methods, we can apply gradient ascent to find ρ * . In the k-th iteration, we update ρ k by ρ k+1 = ρ k + η∇ ρ J(ρ). The exact gradient of J(ρ) with respect to ρ is given by

∇ ρ J(ρ) = p(θ|ρ)p(τ |θ)∇ ρ log p(θ|ρ)R(τ )dτ dθ.

To approximate ∇ ρ J(ρ), we first sample N policy parameters {θ i } from p(θ|ρ). Then we sample one trajectory τ i for each θ i and use the following empirical average to approximate ∇ ρ J(ρ)

∇ ρ J(ρ) = 1 N N i=1 ∇ ρ log p(θ i |ρ) H h=0 γ h r(s i h , a i h ) := 1 N N i=1 g(τ i |ρ), (A.2)

where γ ∈ [0, 1) is the discount factor. Compared with the PGT/GPOMDP estimator in Section 2, the likelihood term ∇ ρ log p(θ i |ρ) in (A.2) for PGPE is independent of horizon H.

Algorithm 1 can be directly applied to the PGPE setting, where we replace the policy parameter θ with the hyper-parameter ρ. When we need to sample N trajectories, we first sample N policy parameters {θ i } from p(θ|ρ). Since the policy is deterministic with given θ i , we sample one trajectory τ i from each policy p(τ |θ i ). The recursive semi-stochastic gradient is given by

v s+1 t = 1 B B j=1 g(τ j |ρ s+1 t ) - 1 B B j=1 g ω (τ j |ρ s+1 t-1 ) + v s+1 t-1 , (A.3)

where g ω (τ j |ρ s+1 t-1 ) is the gradient estimator with step-wise importance weight defined in the way as in (3.2). We call this variance reduced parameter-based algorithm SRVR-PG-PE, which is displayed in Algorithm 2.

Under similar assumptions on the parameter distribution p(θ|ρ), as Assumptions 4.1, 4.3 and 4.4, we can easily prove that SRVR-PG-PE converges to a stationary point of J(ρ) with O(1/ 3/2 ) sample complexity. In particular, we assume the policy parameter θ follows the distribution p(θ|ρ) and we update our estimation of ρ based on the semi-stochastic gradient in (A.3). Recall the gradient ∇ ρ J(ρ) derived in (A.2). Since the policy in SRVR-PG-PE is deterministic, we only need to make the boundedness assumption on p(θ|ρ). In particular, we assume that 1. ∇ ρ log p(θ|ρ) 2 and ∇ 2 ρ log p(θ|ρ) 2 are bounded by constants in a similar way to Assumption 4.1; 2. the gradient estimator g(τ |ρ) = ∇ ρ log p(θ|ρ) H h=0 γ h r(s h , a h ) has bounded variance; 3. and the importance weight ω(τ j |ρ s+1 t-1 , ρ s+1 t ) = p(θ j |ρ s+1 t-1 )/p(θ j |ρ s+1 t

) has bounded variance in a similar way to Assumption 4.4.

Then the same gradient complexity O(1/ 3/2 ) for SRVR-PG-PE can be proved in the same way as the proof of Theorem 4.5 and Corollary 4.7. Since the analysis is almost the same as that of SRVR-PG, we omit the proof of the convergence of SRVR-PG-PE. In fact, according to the analysis in Zhao et al. (2011); Metelli et al. (2018) , all the three assumptions listed above can be easily verified under a Gaussian prior for θ and a linear deterministic policy. Sample one trajectory τ i from each policy π θi 6:

v s+1 0 = ∇ ρ J(ρ s ) := 1 N N i=1 g(τ i | ρ s ) 7: ρ s+1 1 = ρ s+1 0 + ηv s+1 0 8: for t = 1, . . . , m -1 do 9: Sample B policy parameters {θ j } from p(•|ρ s+1 t ) 10:

Sample one trajectory τ j from each policy π θj 11:

v s+1 t = v s+1 t-1 + 1 B B j=1 g τ j |ρ s+1 t -g ω τ j |ρ s+1 t-1 12: ρ s+1 t+1 = ρ s+1 t + ηv s+1 t 13:

end for 14: end for 15: return ρ out , which is uniformly picked from {ρ s t } t=0,...,m;s=0,...,S

B PROOF OF THE MAIN THEORY:

In this section, we provide the proofs of the theoretical results for SRVR-PG (Algorithm 1). Before we start the proof of Theorem 4.5, we first lay down the following key lemma that controls the variance of the importance sampling weight ω. Lemma B.1. For any

θ 1 , θ 2 ∈ R d , let ω 0:h τ |θ 1 , θ 2 = p(τ h |θ 1 )/p(τ h |θ 2 )

, where τ h is a truncated trajectory of τ up to step h. Under Assumptions 4.1 and 4.4, it holds that

Var ω 0:h τ |θ 1 , θ 2 ≤ C ω θ 1 -θ 2 2 2 ,

where

C ω = h(2hG 2 + M )(W + 1).

Recall that in Assumption 4.4 we assume the variance of the importance weight is upper bounded by a constant W . Based on this assumption, Lemma B.1 further bounds the variance of the importance weight via the distance between the behavioral and the target policies. As the algorithm converges, these two policies will be very close and the bound in Lemma B.1 could be much tighter than the constant bound.

Proof of Theorem 4.5. By plugging the definition of the projection operator in (3.3) into the update rule θ s+1 t+1 = P Θ θ s+1 t + ηv s+1 t , we have

θ s+1 t+1 = argmin u∈R d 1 Θ (u) + 1/(2η) u -θ s+1 t 2 2 -v s+1 t , u . (B.1)

Similar to the generalized projected gradient G η (θ) defined in (3.4), we define G s+1 t to be a (stochastic) gradient mapping based on the recursive gradient estimator v s+1 t :

G s+1 t = 1 η θ s+1 t+1 -θ s+1 t = 1 η P Θ θ s+1 t + ηv s+1 t -θ s+1 t . (B.2)

The definition of

G s+1 t differs from G η (θ s+1 t

) only in the semi-stochastic gradient term v s+1 t , while the latter one uses the full gradient ∇J(θ s+1 t ). Note that 1 Θ (•) is convex but not smooth. We assume that p ∈ ∂ 1 Θ (θ s+1 t+1 ) is a sub-gradient of 1 Θ (•). According to the optimality condition of (B.1), we have p + 1/η(θ s+1 t+1 -θ s+1 t ) -v s+1 t = 0. Further by the convexity of 1 Θ (•), we have

1 Θ (θ s+1 t+1 ) ≤ 1 Θ (θ s+1 t ) + p, θ s+1 t+1 -θ s+1 t = 1 Θ (θ s+1 t ) -1/η(θ s+1 t+1 -θ s+1 t ) -v s+1 t , θ s+1 t+1 -θ s+1 t . (B.3)

By Proposition 4.2, J(θ) is L-smooth, which by definition directly implies

J θ s+1 t+1 ≥ J θ s+1 t + ∇J θ s+1 t , θ s+1 t+1 -θ s+1 t - L 2 θ s+1 t+1 -θ s+1 t 2 2 .

For the simplification of presentation, let us define the notation Φ(θ) = J(θ) -1 Θ (θ). Then according to the definition of 1 Θ we have argmax θ∈R d Φ(θ) = argmax θ∈Θ J(θ) := θ * . Combining the above inequality with (B.3), we have

Φ θ s+1 t+1 ≥ Φ θ s+1 t + ∇J θ s+1 t -v s+1 t , θ s+1 t+1 -θ s+1 t + 1 η - L 2 θ s+1 t+1 -θ s+1 t 2 2 = Φ θ s+1 t + ∇J θ s+1 t -v s+1 t , η G s+1 t + η G s+1 t 2 2 - L 2 θ s+1 t+1 -θ s+1 t 2 2 ≥ Φ θ s+1 t - η 2 ∇J θ s+1 t -v s+1 t 2 2 + η 2 G s+1 t 2 2 - L 2 θ s+1 t+1 -θ s+1 t 2 2 = Φ θ s+1 t - η 2 ∇J θ s+1 t -v s+1 t 2 2 + η 4 G s+1 t 2 2 + 1 4η - L 2 θ s+1 t+1 -θ s+1 t 2 2 ≥ Φ θ s+1 t - η 2 ∇J θ s+1 t -v s+1 t 2 2 + η 8 G η θ s+1 t 2 2 - η 4 G η (θ s+1 t ) -G s+1 t 2 2 + 1 4η - L 2 θ s+1 t+1 -θ s+1 t 2 2 , (B.4)

where the second inequality holds due to Young's inequality and the third inequality holds due to the fact that G η (θ s+1 t

) 2 2 ≤ 2 G s+1 t 2 2 + 2 G η (θ s+1 t ) -G s+1 t 2 2 . Denote θs+1 t+1 = prox η 1 Θ (θ s+1 t + η∇J(θ s+1 t ))

. By similar argument in (B.3) we have

1 Θ (θ s+1 t+1 ) ≤ 1 Θ ( θs+1 t+1 ) -1/η(θ s+1 t+1 -θ s+1 t ) -v s+1 t , θ s+1 t+1 -θs+1 t+1 , 1 Θ ( θs+1 t+1 ) ≤ 1 Θ (θ s+1 t+1 ) -1/η(θ s+1 t+1 -θ s+1 t ) -∇J θ s+1 t , θs+1 t+1 -θ s+1 t+1 . Adding the above two inequalities immediately yields θs+1 t+1 -θ s+1 t+1 2 ≤ η ∇J(θ s+1 t ) -v s+1 t 2 , which further implies G η (θ s+1 t ) -G s+1 t 2 ≤ ∇J(θ s+1 t ) -v s+1 t 2 . Submitting this result into (B.4), we obtain Φ θ s+1 t+1 ≥ Φ θ s+1 t - 3η 4 ∇J θ s+1 t -v s+1 t 2 2 + η 8 G η θ s+1 t 2 2 + 1 4η - L 2 θ s+1 t+1 -θ s+1 t 2 2 . (B.5)

We denote the index set of {τ j } B j=1 in the t-th inner iteration by B t . Note that

∇J θ s+1 t -v s+1 t 2 2 = ∇J θ s+1 t -v s+1 t-1 + 1 B j∈Bt g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t 2 2 = ∇J θ s+1 t -∇J(θ s+1 t-1 ) + 1 B j∈Bt g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t + ∇J(θ s+1 t-1 ) -v s+1 t-1 2 2 = ∇J θ s+1 t -∇J(θ s+1 t-1 ) + 1 B j∈Bt g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t 2 + 2 B j∈Bt ∇J θ s+1 t -∇J(θ s+1 t-1 ) + g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t , ∇J(θ s+1 t-1 ) -v s+1 t-1 + ∇J(θ s+1 t-1 ) -v s+1 t-1 2 2 . (B.6)

Conditional on θ s+1 t , taking the expectation over B t yields

E ∇J θ s+1 t -g τ j |θ s+1 t , ∇J(θ s+1 t-1 ) -v s+1 t-1 = 0.

Similarly, taking the expectation over θ s+1 t and the choice of B t yields

E ∇J(θ s+1 t-1 ) -g ω τ j |θ s+1 t-1 , ∇J(θ s+1 t-1 ) -v s+1 t-1 = 0.

Combining the above equations with (B.6), we obtain

E ∇J θ s+1 t -v s+1 t 2 2 = E ∇J θ s+1 t -∇J(θ s+1 t-1 ) + 1 B j∈Bt g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t 2 2 + E ∇J(θ s+1 t-1 ) -v s+1 t-1 2 2 = 1 B 2 j∈Bt E ∇J θ s+1 t -∇J(θ s+1 t-1 ) + g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t 2 2 + E ∇J(θ s+1 t-1 ) -v s+1 t-1 2 2 , (B.7) ≤ 1 B 2 j∈Bt E g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t 2 2 + ∇J(θ s+1 t-1 ) -v s+1 t-1 2 2 , (B.8)

where (B.7) is due to the fact that E x 1 + . . . + x n 2 2 = E x 1 2 + . . . + E x n 2 for independent zero-mean random variables, and (B.8) holds due to the fact that x 1 , . . . ,

x n is due to E x -Ex 2 2 ≤ E x 2 2 .

For the first term, we have

g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t 2 ≤ g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t-1 2 + L θ s+1 t-1 -θ s+1 t

2 by triangle inequality and Proposition 4.2.

E g ω τ j |θ s+1 t-1 -g τ j |θ s+1 t-1 2 2 = E H-1 h=0 (ω 0:h -1) h t=0 ∇ θ log π θ (a i t |s i t ) γ h r(s i h , a i h ) 2 2 = H-1 h=0 E (ω 0:h -1) h t=0 ∇ θ log π θ (a i t |s i t ) γ h r(s i h , a i h ) 2 2 ≤ H-1 h=0 h 2 (2G 2 + M )(W + 1) θ s+1 t-1 -θ s+1 t 2 2 • h 2 G 2 γ h R ≤ 24RG 2 (2G 2 + M )(W + 1)γ (1 -γ) 5 θ s+1 t-1 -θ s+1 t 2 2 , (B.9)

where in the second equality we used the fact that E[∇ log π θ (a|s)] = 0, the first inequality is due to Lemma B.1 and in the last inequality we use the fact that ∞ h=0 h 4 γ h = γ(γ 3 + 11γ 2 + 11γ + 1)/(1 -γ) 5 for |γ| < 1. Combining the results in (B.8) and (B.9), we get Submitting the above result into (B.5) yields

E ∇J θ s+1 t -v s+1 t 2 2 ≤ C γ B θ s+1 t -θ s+1 t-1 2 2 + ∇J(θ s+1 t-1 ) -v s+1 t-1 2 2 ≤ C γ B t l=1 θ s+1 l -θ s+1 l-1 2 2 + ∇J(θ s+1 0 ) -v s+1

E N,B Φ θ s+1 t+1 ≥ E N,B Φ θ s+1 t + η 8 G η θ s+1 t 2 2 + 1 4η - L 2 θ s+1 t+1 -θ s+1 t 2 2 - 3ηC γ 4B E N,B t l=1 θ s+1 l -θ s+1 l-1 2 2 - 3ηξ 2 4N , (B.12)

for t = 1, . . . , m -1.Recall Line 6 in Algorithm 1, where we update θ t+1 1 with the average of a mini-batch of gradients v s 0 = 1/N N i=1 g(τ i | θ s ). Similar to (B.5), by smoothness of J(θ), we have

Φ θ s+1 1 ≥ Φ θ s+1 0 - 3η 4 ∇J θ s+1 0 -v s+1 0 2 2 + η 8 G η θ s+1 0 2 2 + 1 4η - L 2 θ s+1 1 -θ s+1 0 2 2 .

Further by (B.11), it holds that

E Φ θ s+1 1 ≥ E Φ θ s+1 0 - 3ηξ 2 4N + η 8 G η θ s+1 0 2 2 + 1 4η - L 2 θ s+1 1 -θ s+1 0 2 2 . (B.13)

Telescoping inequality (B.12) from t = 1 to m -1 and combining the result with (B.13), we obtain 

E N,B Φ θ s+1 m ≥ E N,B Φ θ s+1 0 + η 8 m-1 t=0 E N G η θ s+1 t 2 2 - 3mηξ 2 4N + 1 4η - L 2 m-1 t=0 θ s+1 t+1 -θ s+1 t 2 2 - 3ηC γ 2B E N,B m-1 t=0 t l=1 θ s+1 l -θ s+1 l-1 2 2 ≥ E N,B Φ θ s+1 0 + η 8 m-1 t=0 E N G η θ s+1 t 2 2 - 3mηξ 2 4N + 1 4η - L 2 - 3mηC γ 2B m-1 t=0 θ s+1 t+1 -θ s+1

η ≤ 1 4L , B m ≥ 3ηC γ L = 72ηG 2 (2G 2 + M )(W + 1)γ M (1 -γ) 2 , (B.15)

and note that θ s+1 0 = θ s , θ s+1 m = θ s+1 , then (B.14) leads to

E N Φ θ s+1 ≥ E N Φ θ s + η 8 m-1 t=0 E N G η θ s+1 t 2 2 - 3mηξ 2 4N . (B.16)

Summing up the above inequality over s = 0, . . . , S -1 yields

η 8 S-1 s=0 m-1 t=0 E G η θ s+1 t 2 2 ≤ E Φ θ S -E Φ θ 0 + 3Smηξ 2 4N , which immediately implies E G η θ out 2 2 ≤ 8 E Φ θ S -E Φ θ 0 ηSm + 6ξ 2 N ≤ 8(Φ(θ * ) -Φ(θ 0 )) ηSm + 6ξ 2 N .

This completes the proof.

Proof of Corollary 4.7. Based on the convergence results in Theorem 4.5, in order to ensure

E ∇J θ out 2 2 ≤ , we can choose S, m and N such that 8(J(θ * ) -J(θ 0 )) ηSm = 2 , 6ξ 2 N = 2 ,

which implies Sm = O(1/ ) and N = O(1/ ). Note that we have set m = O(B). The total number of stochastic gradient evaluations T g we need is

T g = SN + SmB = O N B + B = O 1 3/2 ,

where we set B = 1/ 1/2 .

C PROOF OF TECHNICAL LEMMAS:

In this section, we provide the proofs of the technical lemmas. We first prove the smoothness of the performance function J(θ).

Proof of Proposition 4.2. Recall the definition of PGT in (2.5). We first show the Lipschitzness of g(τ |θ) with baseline b = 0 as follows:

∇g(τ |θ) 2 = H-1 h=0 ∇ 2 θ log π θ (a h |s h ) H-1 t=h γ t r(s t , a t ) 2 ≤ H-1 t=0 γ h ∇ 2 θ log π θ (a t |s t ) 2 R 1 -γ ≤ M R (1 -γ) 2 ,

where we used the fact that 0 < γ < 1. When we have a nonzero baseline b h , we can simply scale it with γ h and the above result still holds up to a constant multiplier.

Since the PGT estimator is an unbiased estimator of the policy gradient ∇ θ J(θ), we have ∇ θ J(θ) = E τ [g(τ |θ)] and thus

∇ 2 θ J(θ) = τ p(τ |θ)∇ θ g(τ |θ)dτ + τ p(τ |θ)g(τ |θ)∇ θ log p(τ |θ)dτ = E τ [∇ θ g(τ |θ)] + E τ [g(τ |θ)∇ θ log p(τ |θ)]. (C.1)

We have already bounded the norm of the first term by M R/(1 -γ) 2 . Now we take a look at the second term. Plugging the equivalent definition of g(τ |θ) in (2.6) yields

E τ [g(τ |θ)∇ θ log p(τ |θ)] = τ H-1 h=0 h t=0 ∇ θ log π θ (a t |s t ) γ h r(s h , a h )∇ θ log p(τ |θ) • p(τ |θ)dτ = τ H-1 h=0 h t=0 ∇ θ log π θ (a t |s t ) γ h r(s h , a h ) H-1 t =0 ∇ θ log π θ (a t |s t ) • p(τ |θ)dτ = τ H-1 h=0 h t=0 ∇ θ log π θ (a t |s t ) γ h r(s h , a h ) h t =0 ∇ θ log π θ (a t |s t ) • p(τ |θ)dτ, (C.2)

where the second equality is due to ∇ θ P (s t +1 |s t , a t ) = 0, and the last equality is due to the fact that for all t > h it holds that

τ h t=0 ∇ θ log π θ (a t |s t )γ h r(s h , a h )∇ θ log π θ (a t |s t ) • p(τ |θ)dτ = 0.

Therefore, we have

E τ [g(τ |θ)∇ θ log p(τ |θ)] 2 ≤ E τ H-1 h=0 h t=0 Gγ h R × (h + 1)G = H-1 h=0 G 2 R(h + 1) 2 γ h ≤ 2G 2 R (1 -γ) 3 .

Putting the above pieces together, we can obtain

∇ 2 θ J(θ) 2 ≤ M R (1 -γ) 2 + 2G 2 R (1 -γ) 3 := L, which implies that J(θ) is L-smooth with L = M R/(1 -γ) 2 + 2G 2 R/(1 -γ) 3 .

Similarly, we can bound the norm of gradient estimator as follows

g(τ |θ) 2 ≤ H-1 h=0 ∇ θ log π θ (a h |s h ) γ h R(1 -γ H-h ) 1 -γ 2 ≤ GR (1 -γ) 2 ,

which completes the proof. 

0:h τ | θ s , θ s+1 t = d 2 p(τ h | θ s )||p(τ h |θ s+1 t ) -1.

To simplify the presentation, we denote θ 1 = θ s and θ 2 = θ s+1 t in the rest of this proof. By definition, we have

d 2 (p(τ h |θ 1 )||p(τ h |θ 2 )) = τ p(τ h |θ 1 ) p(τ h |θ 1 ) p(τ h |θ 2 ) dτ = τ p(τ h |θ 1 ) 2 p(τ h |θ 2 ) -1 dτ.

Taking the gradient of d 2 (p(τ h |θ 1 )||p(τ h |θ 2 )) with respect to θ 1 , we have

∇ θ1 d 2 (p(τ h |θ 1 )||p(τ h |θ 2 )) = 2 τ p(τ h |θ 1 )∇ θ1 p(τ h |θ 1 )p(τ h |θ 2 ) -1 dτ.

In particular, if we set the value of θ 1 to be θ 1 = θ 2 in the above formula of the gradient, we get

∇ θ1 d 2 (p(τ h |θ 1 )||p(τ h |θ 2 )) θ1=θ2 = 2 τ ∇ θ1 p(τ h |θ 1 )dτ θ1=θ2 = 0.

Applying mean value theorem with respect to the variable θ 1 , we have

d 2 (p(τ h |θ 1 )||p(τ h |θ 2 )) = 1 + 1/2(θ 1 -θ 2 ) ∇ 2 θ d 2 (p(τ h |θ)||p(τ h |θ 2 ))(θ 1 -θ 2 ), (C.3)

where θ = tθ 1 + (1 -t)θ 2 for some t ∈ [0, 1] and we used the fact that d

2 (p(τ h |θ 2 )||p(τ h |θ 2 )) = 1.

To bound the above exponentiated Rényi divergence, we need to compute the Hessian matrix. Taking the derivative of ∇ θ1 d 2 (p(τ h |θ 1 )||p(τ h |θ 2 )) with respect to θ 1 further yields

∇ 2 θ d 2 (p(τ h |θ)||p(τ h |θ 2 )) = 2 τ ∇ θ log p(τ h |θ)∇ θ log p(τ h |θ) p(τ h |θ) 2 p(τ h |θ 2 ) dτ + 2 τ ∇ 2 θ p(τ h |θ)p(τ h |θ)p(τ h |θ 2 ) -1 dτ. (C.4)

Thus we need to compute the Hessian matrix of the trajectory distribution function, i.e., ∇ 2 θ p(τ h |θ), which can further be derived from the Hessian matrix of the log-density function. where C ω = h(2hG 2 + M )(W + 1).

∇ 2 θ log p(τ h |θ) = -p(τ h |θ) -2 ∇ θ p(τ h |θ)∇ θ p(τ h |θ) + p(τ h |θ) -1 ∇ 2 θ p(τ h |θ). (C.

D PROOF OF THEORETICAL RESULTS FOR GAUSSIAN POLICY:

In this section, we prove the sample complexity for Gaussian policy. According to (4.1), we can calculate the gradient and Hessian matrix of the logarithm of the policy. The following lemma gives the variance ξ 2 of the PGT estimator, which verifies Assumption 4.3. Lemma D.1 (Lemma 5.5 in Pirotta et al. (2013) ). Given a Gaussian policy π θ (a|s) ∼ N (θ φ(s), σ 2 ), if the |r(s, a)| ≤ R and φ(s) 2 ≤ M φ for all s ∈ S, a ∈ A and R, M φ > 0 are constants, then the variance of PGT estimator defined in (2.5) can be bounded as follows:

Var(g(τ |θ)) ≤ ξ 2 = R 2 M 2 φ (1 -γ) 2 σ 2 1 -γ 2H 1 -γ 2 -Hγ 2H -2γ H 1 -γ H 1 -γ .

Proof of Corollary 4.8. The proof will be similar to that of Corollary 4.7. By Theorem 4.5, to ensure that E[ ∇J(θ out ) 2 2 ] ≤ , we can set 8(J(θ * ) -J(θ 0 ))

ηSm = 2 , 6ξ 2 N = 2 .

Plugging the value of ξ 2 in Lemma D.1 into the second equation above yields N = O( -1 (1-γ) -3 ).

For the first equation, we have S = O(1/(ηm )). Therefore, the total number of stochastic gradient evaluations T g required by Algorithm 1 is Note that C γ = 24RG 2 (2G 2 + M )(W + 1)γ/(1 -γ) 5 . Plugging the values of G, N and L into the above equations yields

T g = SN + SmB = O N ηm + B η .

m = O 1 (1 -γ) 2 √ , B = O 1 (1 -γ) 1 √ .

The corresponding sample complexity is

T g = O 1 (1 -γ) 4 3/2 .

This completes the proof for Gaussian policy.

E ADDITIONAL DETAILS ON EXPERIMENTS:

Now, we provide more details of our experiments presented in Section 5. We first present the parameters for all algorithms we used in all our experiments in Tables 2 and 3 . Among the parameters, the neural network structure and the RL environment parameters are shared across all the algorithms. As mentioned in Section 5, the order of the batch size parameters of our algorithm are chosen according to Corollary 4.7 and we multiply them by a tuning constant via grid search. Similarly, the orders of batch size parameters of SVRPG and GPOMDP are chosen based on the theoretical results suggested by Papini et al. (2018) ; Xu et al. (2019) . Moerover, the learning rates for different methods are tuned by grid search.

We then present the results of PGPE and SRVR-PG-PE on Cartpole, Mountain Car and Pendulum in Figure 2 . In all three environments, our SRVR-PG-PE algorithm shows improvement over PGPE (Sehnke et al., 2010) in terms of number of trajectories. It is worth noting that in all these environments both PGPE and SRVR-PG-PE seem to solve the problem very quickly, which is consistent with the results reported in (Zhao et al., 2011; 2013; Metelli et al., 2018) . Our primary goal in this experiment is to show that our proposed variance reduced policy gradient algorithm can be easily extended to the PGPE framework. To avoid distracting the audience's attention from the variance reduction algorithm on the sample complexity, we do not thoroughly compare the performance of the parameter based policy gradient methods such as PGPE and SRVR-PG-PE with the action based policy gradient methods. We refer interested readers to the valuable empirical studies of PGPE based algorithms presented in Zhao et al. (2011; 2013); Metelli et al. (2018) . 

:

O(•) notation hides constant factors.

] ≤ after O(1/((1 -γ) 4

3/2 )) trajectories in total. Remark 4.9. For Gaussian policy, the number of trajectories Algorithm 1 needs to find anapproximate stationary point, i.e., E[ G η (θ out ) 2 2 ] ≤ , is also in the order of O( -3/2 ), which is faster than PGT and SVRPG. Additionally, we explicitly show that the sample complexity does not depend on the horizon H, which is in sharp contrast with the results inPapini et al. (2018);Xu et al. (2019). The dependence on 1/(1 -γ) comes from the variance of PGT estimator.5 EXPERIMENTSIn this section, we provide experiment results of the proposed algorithm on benchmark reinforcement learning environments including the Cartpole, Mountain Car and Pendulum problems. In all

We slightly abuse the notation by overloading J as the performance function defined on the hyperparameter ρ.

Figures and Tables:

Remark 4.6. Theorem 4.5 states that under a proper choice of step size, batch size and epoch length, the expected squared gradient norm of the performance function at the output of SRVR-PG is in the order of Recall that S is the number of epochs and m is the epoch length of SRVR-PG, so Sm is the total number of iterations of SRVR-PG. Thus the first term O(1/(Sm)) characterizes the convergence rate of SRVR-PG. The second term O(1/N ) comes from the variance of the stochastic gradient used in the outer loop, where N is the batch size used in the snapshot gradient v s+1 0 in Line 5 of SRVR-PG. Compared with the O(1/(Sm) + 1/N + 1/B) convergence rate in Papini et al. (2018), our analysis avoids the additional term O(1/B) that depends on the mini-batch size within each epoch. Compared with Xu et al. (2019), our mini-batch size B is independent of the horizon length H. This enables us to choose a smaller mini-batch size B while maintaining the same convergence rate. As we will show in the next corollary, this improvement leads to a lower sample complexity. Corollary 4.7. Suppose the same conditions as in Theorem 4.5 hold. Set step size as η = 1/(4L), the batch size parameters as N = O(1/ ) and B = O(1/ 1/2 ) respectively, epoch length as m = O(1/ 1/2 ) and the number of epochs as S = O(1/ 1/2 ). Then Algorithm 1 outputs a point θ out that satisfies E[ G η (θ out ) 2 2 ] ≤ within O(1/ 3/2 ) trajectories in total.Note that the results inPapini et al. (2018);Xu et al. (

Figure 1: (a)-(c): Comparison of different algorithms. Experimental results are averaged over 10 repetitions. (d)-(f): Comparison of different batch size B on the performance of SRVR-PG.

Figure 1(f)

Stochastic Recursive Variance Reduced Policy Gradient with Parameter-based Exploration (SRVR-PG-PE) 1: Input: number of epochs S, epoch size m, step size η, batch size N , mini-batch size B, gradient estimator g, initial parameter ρ 0 m := ρ 0 := ρ 0 2: for s = 0, . . . , Sparameters {θ i } from p(•|ρ s ) 5:

t = 1, . . . , m -1, where C γ = 24RG 2 (2G 2 + M )(W + 1)γ/(1 -γ)5 . According to Algorithm 1 and Assumption 4.3, we have E ∇J θ s+1 0

step size η and the epoch length B such that

Lemma C.1 (Lemma 1 inCortes et al. (2010)). Let ω(x) = P (x)/Q(x) be the importance weight for distributions P and Q. ThenE[ω] = 1, E[ω 2 ] = d 2 (P ||Q),where d 2 (P ||Q) = 2 D2(P ||Q) and D 2 (P ||Q) is the Rényi divergence between distributions P and Q. Note that this immediately implies Var(ω) = d 2 (P ||Q) -1. Proof of Lemma B.1. According to the property of importance weight in Lemma C.1, we know Var ω

log π θ (a|s) = (a -θ φ(s))φ(s) σ 2 , ∇ 2 log π θ (a|s) =to see that Assumption 4.1 holds with G = C a M φ /σ 2 and M = M 2 φ /σ 2 .Based on this observation, Proposition 4.2 also holds for Gaussian policy with parameters defined as followsL = RM 2 φ σ 2 (1 -γ) 3 , and C g = RC a M φ σ 2 (1 -γ) 2 . (D.2)

So a good choice of batch size B and epoch length m will lead to Bm = N . Combining this with the requirement of B in Theorem 4.5, we can set m = LN ηC γ , and B = N ηC γ L .

Figure 2: Performance of SRVR-PG-PE compared with PGPE. Experiment results are averaged over 10 runs.

Comparison on sample complexities of different algorithms to achieve ∇J(θ) 2 2 ≤ .

Table content: <table><tr><td>Algorithms</td><td>Complexity</td></tr><tr><td>REINFORCE (Williams, 1992)</td><td/></tr></table>

(p(τ h |θ)||p(τ h |θ 2 )) 2 = 4 log p(τ h |θ)∇ θ log p(τ h |θ) p(τ h |θ) 2 p(τ h |θ 2 ) dτ ∇ θ log p(τ h |θ) 2 2 + 2 ∇ 2 θ log p(τ h |θ) 2 dτ ≤ (4h 2 G 2 + 2hM )E[ω(τ |θ, θ 2 ) 2 ] ≤ 2h(2hG 2 + M )(W + 1),where the second inequality comes from Assumption 4.1 and the last inequality is due to Assumption 4.4 and Lemma C.1. Combining the above result with (C.3), we haveVarω 0:h τ | θ s , θ s+1 t = d 2 p(τ h | θ s )||p(τ h |θ s+1 t ) -1 ≤ C ω θ sθ s+1

Table content: <table><tr><td>Submitting (C.5) into (C.4) yields</td><td/><td/><td/><td/><td/></tr><tr><td colspan="3">∇ 2 θ d 2 + 2</td><td>τ</td><td colspan="2">∇ 2 θ log p(τ h |θ)</td><td>p(τ h |θ) 2 p(τ h |θ 2 )</td><td>dτ</td><td>2</td></tr><tr><td>≤</td><td>τ</td><td colspan="3">p(τ h |θ) 2 p(τ h |θ 2 )</td><td colspan="2">4 t</td><td>2 2 ,</td></tr><tr><td/><td/><td/><td/><td/><td/><td>5)</td></tr></table>

Peer review (iclr2020):
Summary: The paper proposed a policy gradient method called SRVR-PG, which based on stochastic recursive gradient estimator. It shows that the complexity is better than that of SVRPG. Some experiments on standard environments are provided to show the efficiency of the algorithm over GPOMDP and SVRPG. Comments: 1) I think you may ignore the highly-related work: Yang and Zhang, *Policy Optimization with Stochastic Mirror Descent*, June 2019. Since both papers are highly-related, I would suggest the author(s) have some discussions to differentiate two papers. 2) Could you please provide the reasons why you are only choosing two methods (GPOMDP and SVRPG) to compare? There are many policy gradient algorithms such as A3C, A2C, ... which have not mentioned or discussed in this paper. Are they completely different here? Is there no way to compare the performance among them with SRVR-PG? 3) I am not sure if you are using the right word *novel* to describe your method. Basically, you adopt the existing estimator based on SARAH/SPIDER in optimization algorithms into RL problems. I do not think the word *novel* is proper here since it is not something total new. Notice that, the complexity result achieved in this paper is also matched the one for SARAH/SPIDER/SpiderBoost in nonconvex optimization. 4) There are some discussion on choosing a batch-size in the experimental part. However, I do not see the discussion on the learning rate. The choice of parameters for GPOMDP and SVRPG may also need to be discussed. 5) In Papini et al., 2018, for the experiment part, they use a snapshot policy to sample in the inner loop, and use early stopping inner loop. Moreover, they also check variance to recover the backup policy when it is blowup. Do you apply any trick to your experiments? I wonder if your numerical experiments are totally followed on your theory. Minor comments: - Redundancy *)* in eta*S*m in Theorem 4