I. INTRODUCTION

Reinforcement Learning (RL) is a form of implicit stochastic adaptive control where the optimal control policy is estimated without directly estimating the underlying model. This paper considers reinforcement learning for an infinite horizon discounted cost continuous state Markov decision process. In a MDP, actions affect the Markovian state dynamics and result in rewards for the agent. The objective is to find a map from the states to actions, also known as policy, that results in the accumulation of largest expected return. There are many approaches to estimate a policy: policy iteration, Q-learning [1], [2] (which operates in "value" space [3]), policy-gradients [4], [5] (that operate in policy space); see [6], [7].

Recently, policy-gradient algorithms have gained popularity due to their ability to address complex real-world RL problems with continuous state-action spaces. Given a parametrized policy space, usually designed to incorporate domain knowledge, policy-gradient algorithms update policy parameters along an estimated ascent direction of the expected return. Depending on whether the expected reward or the value function is convex or non-convex, the parameters converge to a minimum or a stationary point; for a comprehensive survey see [8], [9].

Typically, to compute the ascent direction in policy search [10], one employs the Policy Gradient Theorem [7] to write the gradient as the product of two factors: the Q-function 1 and the score function (a likelihood ratio). This score function approach has yielded numerous policy search techniques [11], [12], [13], [7], although the resulting gradient estimates are afflicted with high variance: the score function is a martingale and so for a Markov process its variance is O(N ) for N measurements. In pursuit of reducing the variance, we propose replacing the score function with the weak derivatives; see [14] for a textbook treatment. 2 Weak and measured-valued derivatives have been used for real-time reinforcement learning of constrained average cost MDPs (with finite action spaces) in [16], [17], [18], [19], [20]. These papers derive constant step size policy gradient algorithms and show analytically and via numerical examples that substantial variance reduction can be achieved compared to the score function method; moreover the optimal (randomized) policy can be tracked over time when the unknown constrained MDP parameters evolve.

In comparison to [16], [20], this paper considers off-line (decreasing step size) reinforcement learning for continuous state-continuous action infinite horizon discounted cost MDPs when the underlying system can be simulated using statistically independent trials with different policies. To estimate the Q-function in the policy gradient [7], we use Monte Carlo roll-outs with random path lengths akin to [21], motivated by the fact that obtaining unbiased estimates of continuous stateaction Q-function in the infinite horizon case is otherwise challenging. The product of these terms yields a valid estimate of the overall policy gradient, as in [7].

Our main results are: 1) A decreasing step size policy gradient algorithm using Jordan decomposition for the policy gradient. We establish that the resulting policy gradient algorithm, named Policy Gradient with Jordan Decomposition (PG-JD), yields unbiased estimates of the gradient of the reward function. 2) to establish that the PG-JD algorithm converges to a stationary point of the parametrized value function almost surely under decreasing step-sizes. 3) to derive the iteration (and samplefoot_2 ) complexity as

where k is the time step. This shows that the convergence rate is similar to stochastic gradient method for non-convex settings. 4) to upper-bound the expected variance of the gradient estimates obtained using the PG-JD algorithm, which isshown to be lower than those generated by score function methods using Monte Carlo roll-outs with random path lengths, for common policy parametrizations. The setup and problem formulation are discussed in Sec. II. The new policy gradient theorem using weak derivatives (Jordan decomposition) is derived in Sec. III. The algorithm to compute the stochastic gradient and the policy parameter update is given in Sec. IV. Convergence analysis of the stochastic gradient ascent algorithm and its statistical properties are derived in Sec. V. Numerical studies on OpenAI gym using the pendulum environment is discussed in Sec. VI.

II. PROBLEM FORMULATION AND POLICY SEARCH

The problem of reinforcement learning is considered in the framework of Markov Decision Process, which is defined as a tuple (X , A, T , r, γ) consisting of the state space X ⊆ R p , a subset of Euclidean space with elements x ∈ X ; the action space A ⊆ R q , a subset of Euclidean space with elements a ∈ A; the transition law T , a probability density function T (•|a, x) ∈ P(X ) that assigns a next-state upon taking action a in state x, where P(X ) denotes the set of all probability measures on X ; the reward function r(x, a), a real valued function on the product space X × A; the discount γ ∈ (0, 1), a parameter that scales the importance of future rewards.

A stochastic Markov policy µ = {µ k } is defined as a sequence of transition probabilities from X to A such that µ k (D(x)|x) = 1 for each x ∈ X and k = 0, 1, • • • . Here D maps each x ∈ X to the set of all available actions D(x). Let Σ denote the class of stochastic Markov policies.

For an initial state x 0 and a stochastic Markov policy µ ∈ Σ, define the expected reward function

For an initial state x 0 and a Markov policy µ ∈ Σ, using Ionescu Tulcea theorem [22], [23], define P x0 µ as

Here µ 0 ∈ P(X ) is an atomic measure with µ 0 (x 0 ) = 1. The expectation E x0 µ in (1) is with respect to P x0 µ in (2). Our goal is to find the policy µ that maximizes the long-term reward accumulation, or value:

For the infinite horizon problem (3), it is sufficient [24], [25], [23], [26] to restrict the class Σ of policies to the class Σ s ⊂ Σ of stationary stochastic Markov policies. A stationary stochastic Markov policy µ(= {µ}) ∈ Σ s is defined as the transition probability from X to A such that µ(D(x)|x) = 1 for each x ∈ X . In order to solve (3) we resort to direct policy search over the space of continuous stationary policies. It is convenient to parametrize the stationary policy µ(•|•) as µ θ (•|•) for θ ∈ Θ ⊆ R d , for d ∈ N, and search over the space of θ. For example, consider Gaussian policy µ θ (•|x) = N (θ φ(x), σ 2 ). Here the function φ(•) is commonly referred to as the feature map and σ denotes the standard deviation. With a slight abuse of notation, the problem (3) can be reformulated in terms of the finding a parameter vector θ to satisfy:

Here E x0 µ θ is the expectation with respect to the measure induced by the probability measure as in (2) with the policy µ θ = {µ θ } and initial state x 0 .

III. POLICY GRADIENT THEOREM VIA HAHN-JORDAN

The foundation of any valid policy search technique is a valid ascent direction on the value function with respect to the policy parameters. Classically, one may derive that the policy gradient decomposes into two factors: the action-value (Q) function and the score function [4]. Here we establish that one may obviate the need for the log trick that gives rise to the score function through measure-valued differentiation by employing the Jordan decomposition of signed measures [15]. To begin doing so, define the Q-function as

The weak derivative of the signed measure ∇µ θ (•|x) using Jordan decomposition 4 is given as

Here the decomposed positive and negative component measures µ ⊕ θ (•|x) and µ θ (•|x) are orthogonal in L 2 (see Example 1 below). The ergodic measure associated with the transition kernel T (•|x 0 , a 0 ) and policy

. Using this measure (weak) derivative representation of the policy, we can write the gradient of the value function with respect to policy parameters θ in an unusual way which is given in the following theorem. [15] [Hahn Decomposition] Let µ be a finite signed measure on the measurable space (Ω, F ). There exists a disjoint partition of the set Ω into Ω + and Ω -such that Ω = Ω + ∪ Ω -, µ(A) ≥ 0, ∀A ⊂ Ω + , and µ(B) ≤ 0, ∀B ⊂ Ω -.

Result 2. [15]

[Jordan Decomposition] Every finite signed measure µ has a unique decomposition into a difference µ = µ +µ -of two finite non-negative measures µ + and µ -such that for any Hahn decomposition

Theorem 1. (Jordan Decomposition for Policy Gradients)

The policy gradient using Jordan decomposition takes the form

)

where g(θ, x) is a normalizing constant to ensure µ ⊕ and µ are valid measures.

Discussion: The proof can be found in [27]. Theorem 1 is the policy gradient theorem using weak derivatives, specifically Jordan decomposition. In Theorem 1, note that the Q functions in the expectations are the same, indicating that the model is unaffected by the measure decomposition; only the induced measures are different. The expression for the gradient in (7) contains a difference of two expectations. Unlike, the method of score functions, the expectation obviates the need for a score function term. Intuitively, this allows us to avoid computing the logarithm of the policy which may amplify useless parts of the state-action space and cause variance to needlessly be increased, and instead yield a sharp "perceptron-like" behavior. In subsequent sections, we indeed establish that this representation may reduce variance but this reduction intrinsically depends on the policy parameterization. Note that g(θ, x) for a given parameter θ and state x, is a constant, which makes the stochastic gradient easier to compute in Algorithm 2. Before continuing, we present a representative example; see [14] for several examples.

Example 1.

Consider a gaussian policy µ θ (•|x) = N (θ φ(x), σ 2 ), where the mean of the gaussian distribution is modulated by the optimization parameter. The Jordan decomposition of the gaussian policy can be derived as follows:

Here we may glean the normalizing constant g(θ, x) = φ(x)

and the positive and negative component measures are

Observe that µ ⊕ θ (•|x) and µ θ (•|x) define the Rayleighfoot_4 policy. They are orthogonal in the sense that µ ⊕ θ (•|x) is defined on 6χ(a > θ φ(x)) and µ θ (•|x) is defined over χ(a < θ φ(x)).

IV. POLICY SEARCH VIA JORDAN DECOMPOSITION

In order to develop a policy search method based on Theorem 1, we need samples of both factors inside the expectation in (7) which are unbiased. We first focus on the later factor, the Q-function.

A. Estimating the Action-Value

The estimation of the Q-function is carried out using Monte Carlo roll-outs of random path lengths, similar to [21].

Here the random length is a geometric random variable with parameter γ, the discount factor in the reinforcement learning problem. Specifically, we simulate T ∼ Geom (1γ) and then simulate state-action pairs according to the positive and negative induced policies π ⊕ and π . For this time horizon, we collect rewards for the two different trajectories.

More specifically, from a given starting state x 0 , a (real) trajectory is simulated to update the policy parameters θ. At each epoch k of the parameter update θ k , the simulator (modeled as (S(= X ), A, T , r, γ)) is called two times to simulate two different (phantom 7 ) trajectories. These trajectories correspond to the random Monte-Carlo roll-outs used to estimate the Q-functions with two different policies, the positive and negative policy measure, and hence the stochastic gradient of the expected reward function. Let T denote a geometrically distributed random variable: T ∼ Geom(1γ) where γ is the discount factor. Let the path-wise cost be defined by

. Discussion: Algorithm 2 with Algorithm 1 is the stochastic gradient algorithm that is used to update the policy parameters. The simulation consists of a single simulation (real trajectory) to update the parameters and multiple phantom simulations to estimate the gradient of the expected reward function. The two phantom trajectories correspond to different polices and not different models, starting from the system's state represented by the state corresponding to the real trajectory. The stochastic gradient computation is summarized in three steps: For a fixed initial state-(i) Simulate two phantom initial actions from the measures obtained using Jordan decomposition, i.e, µ θ k (•|s 0 ) and µ ⊕ θ k (•|s ⊕ 0 ). (ii) Simulate a geometric random variable T k , and (iii) Perform Monte Carlo roll-outs of length T k -1 (i.e, simulate and feed actions to the simulator and collect the rewards) using the system policy derived from old parameters, i.e using

}. The merit of using these random horizons for estimation of the Q function, as summarized in Algorithm 1, is that one may establish that it is an unbiased estimate in the infinitehorizon discounted case, as we summarize in the following theorem.

Theorem 2. For a geometric r.v T , let the approximate state-action value function (Q-function) be defined by

T denote a geometrically distributed random variable. Then,

The proof can be found in [27]. Now that we may obtain unbiased samples of the action-value function, we shift focus to how to compute the stochastic gradients needed for policy search based on Jordan decomposition (Theorem 1).

B. Stochastic Gradient Algorithm

With the estimation of the action-value function addressed, we now discuss how we can sample the former factor: the signed measure gradients. Specifically, Theorem 1 can be used to effectively compute the gradient given access to an oracle/simulator that may generate state-action-reward triples. It is well known that one only needs to compute estimates of the gradient that are unbiased in expectation to ensure convergence of the iterates to a stationary point [7]. This results in a modification of the gradient expression as in REINFORCE algorithm [11], [7], which is a stochastic gradient, for computing the optimal policy of the reinforcement learning problem. Let E T denote the expectation with respect to the geometric distribution.

Using Theorem 2 and Fubini's Theorem [28], the gradient in (7) can be rewritten to make it implementable on a simulator:

We have from Theorem 2 and ( 14),

Here the initial state simulated from the ergodic measure is x 0 ∼ π µ θ (x), and the policies that simulate the two trajectories are:

Here the initial actions are simulated from the decomposed measures and the parametrized policy is used for the remainder of the trajectory simulation. Here (15) is the (stochastic) gradient estimate for a random path length T and ( 16) is the (stochastic) gradient estimate Algorithm 2 Policy Gradient with Jordan Decomposition (PG-JD)

Input: System state x k+1 , parameter vector θ k , and continuous random policy µ θ k .

Output: Parameter θ k+1 and next system input a k+1 ∼ µ θ k+1 .

Step 1. Simulate T k ∼ Geom(1γ), i.e., P (T k = t) = (1γ)γ t . Define the initial conditions: s ⊕ 0 , s 0 = x k+1 . Define: μ⊕

} as the policy for trajectory 1.

Define:

} as the policy for trajectory 2.

Step 2. Simulate

using a realization T z . Using the estimates ( 16) that are computable using Algorithm 1 to estimate the Q function with respect to the signed measures, then, we may write out an iterative stochastic gradient method to optimize θ with respect to the value function as

The overall policy search routine is summarized as Algorithm 2. Its convergence and variance properties are discussed in the following section.

V. CONVERGENCE, COMPLEXITY, & VARIANCE ANALYSIS

In this section, we discuss a few properties of the stochastic gradient ascent algorithm derived using weak derivatives, namely, convergence, the iteration complexity, sample complexity, and the variance of the resulting gradient estimates.

A. Convergence Analysis

We now analyze the convergence of the PG-JD algorithm (Algorithm 2 ), establishing that the stochastic gradient estimates obtained from the algorithm are unbiased estimates of the true gradient, and that the parameter sequence (17) converges almost surely to a stationary point of the value function (4). To do so, some assumptions are required which we state next.

1) Assumptions: (i) The reward functionfoot_7 r(x, a) is bounded Lipschitz, i.e,

ψ-irreducible, positive Harris recurrent, and geometrically ergodic. (iv) The continuous policy µ θ (a|x) is Lipschitz, i.e,

for all θ ∈ Θ, and n, m > 0. Assumptions (i) -(iii) are model assumptions, whereas Assumptions (iv) -(vi) impose restricts about how the algorithm behaves. Assumption (i) is standard, and tied to learnability of the problem. Assumption (ii) is a continuity assumption on the transition law that is easily satisfied by most physical systems. Assumption (iii) makes sure that for every policy µ θ , there exists a unique invariant (stationary) measure and the Markov chain reaches stationarity geometrically fast; see [30]. All the results hold without the transition law being geometrically ergodic. Assuming geometric ergodicity makes simulating from the ergodic measure (in Algorithm 2, Sec.IV) more meaningful. Regarding the algorithmic conditions: Assumptions (iv)-(v) are standard in stochastic gradient methods; see [31]. Assumption (vi) says that the stochastic gradient is always bounded by the true gradient, which can grow unbounded with θ. This assumption makes sure that the martingale noise of the stochastic gradient is bounded by the true gradient; see [31].

Proposition 1. Under Assumption (i), the expected cost J(θ) in the reinforcement learning problem ( 4) is a bounded realvalued function, i.e,

The following result makes sure that the stochastic gradient estimates so obtained are representative of the true gradient. Theorem 3. The stochastic gradient obtained in ( 16) is an unbiased estimate of the true gradient ∇J(θ), i.e, E ∇J(θ) = ∇J(θ).

(

Discussion: The proof can be found in [27]. Theorem 3 says that the estimates of the stochastic gradient are unbiased 9 As in [29], K(υ, ν) denotes the Kantorovich distance between probability distributions υ and ν. It is given by:

in expectation. This is required to ensure the almost sure convergence of the iterates to a stationary point [7]. (20)

Discussion: The proof can be found in [27]. The expected cost function J(θ), under model assumptions, is continuous and L-Lipschitz; see [ Chapter 7] [32] and [29]. Theorem 4 says that the sequence of iterates {θ k } converges to θ * with probability one, and since J(θ) is a continuous function, J(θ k ) converges to J(θ * ) with probability one. The gradient (which can be unbounded) at iterates {θ k } is such that ∇J(θ * ) = 0 with probability one.

B. Sample Complexity

In this section, we consider the convergence rate analysis of the PG-JD algorithm. We choose the stepsize to be k = k -b for some parameter b ∈ (0, 1). Since the optimization of J(θ) is generally non-convex, we consider the convergence rate in terms of a metric of non-stationarity, i.e., the norm of the gradient ∇J(θ) 2 . The following theorem considers a step-size that diminishes more slowly than Assumption (v), which yields a O(1/ √ k) rate for the decrement of the expected gradient norm square E ∇J(θ k ) 2 . Theorem 5. Let θ k k≥0 be the sequence of parameters of the policy µ θ k generated by Algorithm 2. Let the stepsize be k = k -b for b ∈ (0, 1) and ∆ = min ε, η for some ε, η > 0. Let

denote the number of iteration steps for the norm of the expected cost to come within the error neighborhood. Then, under Assumptions (i) -(iv), (vi)

where optimizing the complexity bound over b, we have b = 1/2. Therefore, K ∆ = O(∆ -2 ).

Discussion: See [27] for proof. Theorem 5 characterizes the iteration complexity, which is a measure of the number of iteration steps of the algorithm are required to settle down on a stationary point of the value function. The iteration complexity is O(1/ √ k) showing that the convergence rate is similar to the stochastic gradient methods for non-convex settings. We emphasize that despite the re-invention of such proofs in machine learning, in fact much more general analysis has been developed in the literature. For example Proposition 5, pg 294 in [33] gives general L q error bounds for q/geq2 in the presence of Markovian noise (whereas our setting is i.i.d. noise). Our intention is to identify these rates in the context and language of modern RL. Corollary 6. Let γ denote the discount factor and K ∆ denote the iteration complexity. The average sample complexity M ∆ γ using Algorithm 2 is given as:

Discussion: The proof can be found in [27]. Corollary 6 characterizes the sample complexity, which is a measure of the number of the expected total number of actions and states realized. Higher the discount factor γ, longer the two (random) Monte-Carlo roll-outs (trajectories) that need to simulated, and hence higher the sample complexity. Together the complexity results, Theorem 5 and Corollary 6, provide an estimate of the duration and expected number of simulations to learn a stationary solution for the reinforcement learning task considered.

C. Variance Analysis

In this section, we provide an analysis of the variance of the stochastic gradient estimates obtained using weak derivatives and score function approaches. Since the Q-function estimation in the computation of the gradient is performed using random Monte Carlo roll-outs as in [21], the stochastic gradient obtained is a function of the geometric random variable T that characterizes the roll-out (trajectory) length. To obtain a comparison of the different methods -weak derivatives and score function -we consider the expected variance of the gradient estimates. The proof of Theorem 7 can be found in [27]. The proof of Theorem 8 is similar and hence omitted.

Theorem 7. The expected variance of the gradient estimates ∇J obtained using weak derivatives is given as:

where

Theorem 8. The expected variance of the gradient estimates ∇J, if score function is used instead of weak derivatives, is given as:

where G SF = E (x,a)∼µ θ (a|x) ∇µ θ (a|x) 2 .

Corollary 9. For the Gaussian policy µ θ (•|x) = N (θ φ(x), σ 2 ), we have

Hence, the maximum expected variance of the gradient estimates using weak derivatives is smaller than those obtained using the score function method. k=0 γ k r(x k , a k )} is evaluated over 50 trajectories with γ = 0.97. Observe that the discounted return is higher on average using Monte-Carlo PG-JD as opposed to PG-SF. It can be attributed to algorithm iterates converging to a "better" stationary point due to smaller variance in the gradient estimates.

VI. NUMERICAL STUDIES

In this section, we present a simple experiment using PG-JD algorithm on the Pendulum environment in OpenAI gym [34]. The performance is compared with Monte Carlo Policy Gradient using Score Function (PG-SF) which is akin to REINFORCE [35] with random roll-out horizons; see Fig. 1. In the simulation environment, the pendulum starts at a random position, and the goal is to swing it up so that it stays upright. The environment state is a vector of dimension three, i.e., x k = (cos(ϕ k ), sin(ϕ k ), φk ) , where ϕ k is the angle between the pendulum and the upright direction, and φk is the derivative of ϕ k . The action a k is a one-dimensional scalar modified using a tanh-function, and represents the joint effort. The received reward r(x k , a k ) is given as

which lies in [-16.2734, 0], ϕ k is normalized between [-π, π] and a k lies in [-2, 2]. The transition dynamics are specified by Newton's Second Law of Motion. We use Gaussian policy π θ , which is parameterized as π θ (•|x) = N (θ T φ(x), σ 2 ), where σ = 1.0 and φ(x)(= x) being the feature vector. The policy is a stationary policy (timehomogeneous) as it is well known [6] to be sufficient for infinite or random horizon discounted MDP problems.

Observe that the discounted return is higher on average using PG-JD as opposed to PG-SF, which may attributable to the variance-reduced properties of the policy gradient estimates using signed measures as compared with the score function.

Remark: It is noted that for common parametrizations of the mean of the Gaussian policy [12], for example like linear -θ T φ(s), the score function is unbounded with respect to θ with the expression being (a-θ T φ(s)) σ 2 φ(s). This results in convergence issues in policy gradient algorithms for unbounded θ and unbounded state spaces. However, using Jordan decomposition, even with linear parametrization and unboundedness, the convergence of the policy gradient algorithm is ensured due to the absence of explicit function of θ.