TRULY DETERMINISTIC POLICY OPTIMIZATION
Year: 2020

Abstract:
In this paper, we present a policy gradient method that avoids exploratory noise injection and performs policy search over the deterministic landscape. By avoiding noise injection all sources of estimation variance can be eliminated in systems with deterministic dynamics (up to the initial state distribution). Since deterministic policy regularization is impossible using traditional non-metric measures such as the KL divergence, we derive a Wasserstein-based quadratic model for our purposes. We state conditions on the system model under which it is possible to establish a monotonic policy improvement guarantee, propose a surrogate function for policy gradient estimation, and show that it is possible to compute exact advantage estimates if both the state transition model and the policy are deterministic. Finally, we describe two novel robotic control environments—one with non-local rewards in the frequency domain and the other with a long horizon (8000 time-steps)— for which our policy gradient method (TDPO) significantly outperforms existing methods (PPO, TRPO, DDPG, and TD3). Policy Gradient (PG) methods can be broadly characterized by three defining elements: the policy gradient estimator, the regularization measures, and the exploration profile. For gradient estimation, episodic (Williams, 1992), importance-sampling-based (Schulman et al., 2015a), and deterministic (Silver et al., 2014) gradients are some of the most common estimation oracles. As for regularization measures, either an Euclidean distance within the parameter space (Williams, 1992; Silver et al., 2014; Lillicrap et al., 2015), or dimensionally consistent non-metric measures (Schulman et al., 2015a; Kakade & Langford, 2002; Schulman et al., 2017; Kakade, 2002; Wu et al., 2017) have been frequently adapted. Common exploration profiles include Gaussian (Schulman et al., 2015a) and stochastic processes (Lillicrap et al., 2015). These elements form the basis of many model-free and stochastic policy optimization methods successfully capable of learning high-dimensional policy parameters. Both stochastic and deterministic policy search can be useful in applications. A stochastic policy has the effect of smoothing or filtering the policy landscape, which is desirable for optimization. Searching through stochastic policies has enabled the effective control of challenging environments under a general framework (Schulman et al., 2015a; 2017). The same method could either learn robotic movements or play basic games (1) with minimal domain-specific knowledge, (2) regardless of function approximation classes, and (3) with less human intervention (ignoring reward engineering and hyper-parameter tuning) (Duan et al., 2016). Using stochasticity for exploration, although it imposes approximations and variance, has provided a robust way to actively search for higher rewards. Despite many successes, there are practical environments which remain challenging for current policy gradient methods. For example, non-local rewards (e.g., those defined in the frequency domain), long time horizons, and naturally-resonant environments all occur in realistic robotic systems (Kuo & Golnaraghi, 2002; Meirovitch, 1975; Preumont & Seto, 2008) but can present issues for policy gradient search. To tackle challenging environments such as these, this paper considers policy gradient methods based on deterministic policies and deterministic gradient estimates, which could offer advantages by allowing the estimation of global reward gradients on long horizons without the need to inject noise into the system for exploration. To facilitate a dimensionally consistent and low-variance deterministic policy search, a compatible policy gradient estimator and a metric measure for regularization should be employed. For gradient estimation we focus on Vine estimators (Schulman et al., 2015a), which can be easily applied to deterministic policies. As a metric measure we use the Wasserstein distance, which can measure meaningful distances between deterministic policy functions that have non-

text: In this paper, we present a policy gradient method that avoids exploratory noise injection and performs policy search over the deterministic landscape. By avoiding noise injection all sources of estimation variance can be eliminated in systems with deterministic dynamics (up to the initial state distribution). Since deterministic policy regularization is impossible using traditional non-metric measures such as the KL divergence, we derive a Wasserstein-based quadratic model for our purposes. We state conditions on the system model under which it is possible to establish a monotonic policy improvement guarantee, propose a surrogate function for policy gradient estimation, and show that it is possible to compute exact advantage estimates if both the state transition model and the policy are deterministic. Finally, we describe two novel robotic control environments—one with non-local rewards in the frequency domain and the other with a long horizon (8000 time-steps)— for which our policy gradient method (TDPO) significantly outperforms existing methods (PPO, TRPO, DDPG, and TD3).
Policy Gradient (PG) methods can be broadly characterized by three defining elements: the policy gradient estimator, the regularization measures, and the exploration profile. For gradient estimation, episodic (Williams, 1992), importance-sampling-based (Schulman et al., 2015a), and deterministic (Silver et al., 2014) gradients are some of the most common estimation oracles. As for regularization measures, either an Euclidean distance within the parameter space (Williams, 1992; Silver et al., 2014; Lillicrap et al., 2015), or dimensionally consistent non-metric measures (Schulman et al., 2015a; Kakade & Langford, 2002; Schulman et al., 2017; Kakade, 2002; Wu et al., 2017) have been frequently adapted. Common exploration profiles include Gaussian (Schulman et al., 2015a) and stochastic processes (Lillicrap et al., 2015). These elements form the basis of many model-free and stochastic policy optimization methods successfully capable of learning high-dimensional policy parameters.
Both stochastic and deterministic policy search can be useful in applications. A stochastic policy has the effect of smoothing or filtering the policy landscape, which is desirable for optimization. Searching through stochastic policies has enabled the effective control of challenging environments under a general framework (Schulman et al., 2015a; 2017). The same method could either learn robotic movements or play basic games (1) with minimal domain-specific knowledge, (2) regardless of function approximation classes, and (3) with less human intervention (ignoring reward engineering and hyper-parameter tuning) (Duan et al., 2016). Using stochasticity for exploration, although it imposes approximations and variance, has provided a robust way to actively search for higher rewards. Despite many successes, there are practical environments which remain challenging for current policy gradient methods. For example, non-local rewards (e.g., those defined in the frequency domain), long time horizons, and naturally-resonant environments all occur in realistic robotic systems (Kuo & Golnaraghi, 2002; Meirovitch, 1975; Preumont & Seto, 2008) but can present issues for policy gradient search.
To tackle challenging environments such as these, this paper considers policy gradient methods based on deterministic policies and deterministic gradient estimates, which could offer advantages by allowing the estimation of global reward gradients on long horizons without the need to inject noise into the system for exploration. To facilitate a dimensionally consistent and low-variance deterministic policy search, a compatible policy gradient estimator and a metric measure for regularization should be employed. For gradient estimation we focus on Vine estimators (Schulman et al., 2015a), which can be easily applied to deterministic policies. As a metric measure we use the Wasserstein distance, which can measure meaningful distances between deterministic policy functions that have non-
overlapping supports (in contrast to the Kullback-Liebler (KL) divergence and the Total Variation (TV) distance).
The Wasserstein metric has seen substantial recent application in a variety of machine-learning domains, such as the successful stable learning of generative adversarial models (Arjovsky et al., 2017). Theoretically, this metric has been studied in the context of Lipschitz-continuous Markov decision processes in reinforcement learning (Hinderer, 2005; Ferns et al., 2012). Pirotta et al. (2015) defined a policy gradient method using the Wasserestein distance by relying on Lipschitz continuity assumptions with respect to the policy gradient itself. Furthermore, for Lipschitz-continuous Markov decision processes, Asadi et al. (2018) and Rachelson & Lagoudakis (2010) used the Wasserstein distance to derive model-based value-iteration and policy-iteration methods, respectively. On a more practical note, Pacchiano et al. (2019) utilized Wasserstein regularization for behavior-guided stochastic policy optimization. Moreover, Abdullah et al. (2019) has proposed another robust stochastic policy gradient formulation. Estimating the Wasserstein distance for general distributions is more complicated than typical KL-divergences (Villani, 2008). This fact constitutes and emphasizes the contributions of Abdullah et al. (2019) and Pacchiano et al. (2019). However, for our deterministic observation-conditional policies, closed-form computation of Wasserstein distances is possible without any approximation.
Existing deterministic policy gradient methods (e.g., DDPG and TD3) use deterministic policies (Silver et al., 2014; Lillicrap et al., 2015; Fujimoto et al., 2018), meaning that they learn a deterministic policy function from states to actions. However, such methods still use stochastic search (i.e., they add stochastic noise to their deterministic actions to force exploration during policy search). In contrast, we will be interested in a method which not only uses deterministic policies, but also uses deterministic search (i.e., without constant stochastic noise injection). We call this new method truly deterministic policy optimization (TDPO) and it may have lower estimation variances and better scalability to long horizons, as we will show in numerical examples.
Scalability to long horizons is one of the most challenging aspects for policy gradient methods that use stochastic search. This issue is sometimes referred to as the curse of horizon in reinforcement learning (Liu et al., 2018). General worst-case analyses suggests that the sample complexity of reinforcement learning is exponential with respect to the horizon length (Kakade et al., 2003; Kearns et al., 2000; 2002). Deriving polynomial lower-bounds for the sample complexity of reinforcement learning methods is still an open problem (Jiang & Agarwal, 2018). Lower-bounding the sample complexity of reinforcement learning for long horizons under different settings and simplifying assumptions has been a topic of theoretical research (Dann & Brunskill, 2015; Wang et al., 2020). Some recent work has examined the scalability of importance sampling gradient estimators to long horizons in terms of both theoretical and practical estimator variances (Liu et al., 2018; Kallus & Uehara, 2019; 2020). All in all, long horizons are challenging for all reinforcement learning methods, especially the ones suffering from excessive estimation variance due to the use of stochastic policies for exploration, and our truly deterministic method may have advantages in this respect.
In this paper we focus on continuous-domain robotic environments with reset capability to previously visited states. The main contributions of this work are: (1) we introduce a Deterministic Vine (DeVine) policy gradient estimator which avoids constant exploratory noise injection; (2) we derive a novel deterministically-compatible surrogate function and provide monotonic payoff improvement guarantees; (3) we show how to use the DeVine policy gradient estimator with the Wasserstein-based surrogate in a practical algorithm (TDPO: Truly Deterministic Policy Optimization); (4) we illustrate the robustness of the TDPO policy search process in robotic control environments with non-local rewards, long horizons, and/or resonant frequencies.

heading: 1 BACKGROUND
text: MDP preliminaries. An infinite-horizon discounted Markov decision process (MDP) is specified by (S,A, P,R, µ, γ), where S is the state space, A is the action space, P : S × A → ∆(S) is the transition dynamics, R : S × A → [0, Rmax] is the reward function, γ ∈ [0, 1) is the discount factor, and µ(s) is the initial state distribution of interest (where ∆(F) denotes the set of all probability distributions over F , otherwise known as the Credal set of F). The transition dynamics P is defined as an operator which produces a distribution over the state space for the next state s′ ∼ P (s, a). The transition dynamics can be easily generalized to take distributions of states or
actions as input (i.e., by having P defined as P : ∆(S)×A → ∆(S) or P : S ×∆(A)→ ∆(S)). We may abuse the notation and replace δs and δa by s and a, where δs and δa are the deterministic distributions concentrated at the state s and action a, respectively. A policy π : S → ∆(A) specifies a distribution over actions for each state, and induces trajectories from a given starting state s as follows: s1 = s, a1 ∼ π(s1), r1 = R(s1, a1), s2 ∼ P (s2, a2), a2 ∼ π(s2), etc. We will denote trajectories as state-action tuples τ = (s1, a1, s2, a2, . . .). One can generalize the dynamics (1) by using a policy instead of an action distribution P(µs, π) := Es∼µs [Ea∼π(s)[P (s, a)]], and (2) by introducing the t-step transition dynamics recursively as Pt(µs, π) := P(Pt−1(µs, π), π) with P0(µs, π) := µs, where µs is a distribution over S. The visitation frequency can then be defined as ρπµ := (1− γ) ∑∞ t=1 γ t−1Pt−1(µ, π). Table 2 of the appendix summarizes all MDP notation.
The value function of π is defined as V π(s) := E[ ∑∞ t=1 γ
t−1rt | s1 = s;π]. Similarly, one can define Qπ(s, a) by conditioning on the first action. The advantage function can then be defined as their difference (i.e. Aπ(s, a) := Qπ(s, a)− V π(s)). Generally, one can define the advantage/value of one policy with respect to another using Aπ(s, π′) := E[Qπ(s, a) − V π(s) | a ∼ π′(·|s)] and Qπ(s, π′) := E[Qπ(s, a) | a ∼ π′(·|s)]. Finally, the payoff of a policy ηπ := E[V π(s); s ∼ µ] is the average value over the initial states distribution of the MDP. Probabilistic and mathematical notations. Sometimes we refer to ∫ f(x)g(x)dx integrals as 〈f, g〉x Hilbert space inner products. Assuming that ζ and ν are two probabilistic densities, the Kulback-Liebler (KL) divergence is DKL(ζ|ν) := 〈ζ(x), log( ζ(x)ν(x) )〉x, the Total-Variation (TV) distance is TV(ζ, ν) =: 12 〈|ζ(x) − ν(x)|, 1〉x, and the Wasserstein distance is W (ζ, ν) = infγ∈Γ(ζ,ν)〈‖x − y‖, γ(x, y)〉x,y where Γ(ζ, ν) is the set of couplings for ζ and ν. We define Lip(f(x, y);x) := supx ‖∇xf(x, y)‖2 and assume the existence of Lip(Qπ(s, a); a) and ‖Lip(∇sQπ(s, a); a)‖2 constants. Under this notation, the Rubinstein-Kantrovich (RK) duality states that the |〈ζ(x)− ν(x), f(x)〉x| ≤W (ζ, ν) · Lip(f ;x) bound is tight for all f . For brevity, we may abuse the notation and denote supsW (π1(·|s), π2(·|s)) with W (π1, π2) (and similarly for other measures). For parameterized policies, we define∇πf(π) := ∇θf(π) where π is parameterized by the vector θ. Table 1 of the appendix summarizes all these mathematical definitions.
Policy gradient preliminaries. The advantage decomposition lemma provides insight into the relationship between payoff improvements and advantages (Kakade & Langford, 2002). That is,
ηπ2 − ηπ1 = 1
1− γ · Es∼ρπ2µ [A π1(s, π2)]. (1)
We will denote the current and the candidate next policy as π1 and π2, respectively. Taking derivatives of both sides with respect to π2 at π1 yields
∇π2ηπ2 = 1
1− γ
[ 〈∇π2ρπ2µ (·), Aπ1(·, π1)〉+ 〈ρπ1µ (·),∇π2Aπ1(·, π2)〉 ] . (2)
Since π1 does not have any advantage over itself (i.e., Aπ1(·, π1) = 0), the first term is zero. Thus, the Policy Gradient (PG) theorem is derived as
∇π2ηπ2 ∣∣∣ π2=π1 = 1 1− γ · Es∼ρπ1µ [∇π2A π1(s, π2)] ∣∣∣ π2=π1 . (3)
For policy iteration with function approximation, we assume π2 and π1 to be parameterized by θ2 and θ1, respectively. One can view the PG theorem as a Taylor expansion of the payoff at θ1.
Conservative Policy Iteration (CPI) (Kakade & Langford, 2002) was one of the early dimensionally consistent methods with a surrogate of the form Lπ1(π2) = ηπ1 + 11−γ · Es∼ρπ1µ [A
π1(s, π2)] − C 2 TV
2(π1, π2). The C coefficient guarantees non-decreasing payoffs. However, CPI is limited to linear function approximation classes due to the update rule πnew ← (1− α)πold + απ′. This lead to the design of the Trust Region Policy Optimization (TRPO) (Schulman et al., 2015a) algorithm.
TRPO exchanged the bounded squared TV distance with the KL divergence by lower bounding it using the Pinsker inequality. This made TRPO closer to the Natural Policy Gradient algorithm(Kakade, 2002), and for Gaussian policies the modified terms had similar Taylor expansions within small trust
regions. Confined trust regions are a stable way of making large updates and avoiding pessimistic coefficients. For gradient estimates, TRPO used Importance Sampling (IS) with a baseline shift:
∇θ2Es∼ρπ1µ [A π1(s, π2)] ∣∣∣ θ2=θ1
= ∇θ2E s∼ρπ1µ a∼π1(·|s)
[ Qπ1(s, a)
π2(a|s) π1(a|s) ]∣∣∣ θ2=θ1 . (4)
While empirical E[Aπ1(s, π2)] and E[Qπ1(s, π2)] estimates yield identical variances in principle, the importance sampling estimator in (4) imposes larger variances. Later, Proximal Policy Optimization (PPO) (Schulman et al., 2015b) proposed utilizing the Generalized Advantage Estimation (GAE) method for variance reduction and incorporated first-order smoothing like ADAM (Kingma & Ba, 2014). GAE employed Temporal-Difference (TD) learning (Bhatnagar et al., 2009) for variance reduction. Although TD-learning was not theoretically guaranteed to converge and could add bias, it improved the estimation accuracy.
As an alternative to importance sampling, deterministic policy gradient estimators were also utilized in an actor-critic fashion. Deep Deterministic Policy Gradient (DDPG) (Lillicrap et al., 2015) generalized deterministic gradients by employing Approximate Dynamic Programming (ADP) (Mnih et al., 2015) for variance reduction. Twin Delayed Deterministic Policy Gradient (TD3) (Fujimoto et al., 2018) improved DDPG’s approximation to build an even better policy optimization method. Although both methods used deterministic policies, they still performed stochastic search by adding stochastic noise to the deterministic policies to force exploration.
Other lines of stochastic policy optimization were later proposed. Wu et al. (2017) used a Kroneckerfactored approximation for curvatures. Haarnoja et al. (2018) proposed a maximum entropy actorcritic method for stochastic policy optimization.

heading: 2 MONOTONIC POLICY IMPROVEMENT GUARANTEE
text: We use the Wasserstein metric because it allows the effective measurement of distances between probability distributions or functions with non-overlapping support, such as deterministic policies, unlike the KL divergence or TV distance which are either unbounded or maximal in this case. The physical transition model’s smoothness enables the use of the Wasserstein distance to regularize deterministic policies. Therefore, we make the following two assumptions about the transition model:
W (P(µ, π1),P(µ, π2)) ≤ Lπ ·W (π1, π2), (5) W (P(µ1, π),P(µ2, π)) ≤ Lµ ·W (µ1, µ2). (6)
Also, we make the dynamics stability assumption sup ∑t k=1 L̂ (k−1) µ,π1,π2 ∏t−1 i=k+1 L̃ (i) µ,π1,π2 <∞, with the definitions of the new constants and further discussion of the implications deferred to Section A.5 of the appendix where we also discuss Assumptions 5 and 6 and the existence of Lip(Qπ(s, a); a).
The advantage decomposition lemma can be rewritten as
ηπ2 = ηπ1 + 1
1− γ · Es∼ρπ1µ [A
π1(s, π2)] + 1
1− γ · 〈ρπ2µ − ρπ1µ , Aπ1(·, π2)〉s. (7)
The 〈ρπ2µ − ρπ1µ , Aπ1(·, π2)〉 term has zero gradient at π2 = π1, which qualifies it to be crudely called “the second-order term”. The theory behind lower-bounding this second-order term with all derivations is left to the appendix. Next, we present the theoretical bottom line.

heading: 2.1 THE MONOTONIC PAYOFF IMPROVEMENT GUARANTEE
text: Combining the results of Inequality (30) and Theorems A.5 and A.4 of the appendix leads us to define the regularization terms and coefficients:
C1 := sup s 2 · Lip(Qπ1(s, a); a) · γ · Lπ (1− γ)(1− γLµ)
C2 := sup s ‖Lip(∇sQπ1(s, a); a)‖2 · γ · Lπ (1− γ)(1− γLµ)
LWG(π1, π2; s) := W (π2(a|s), π1(a|s)) × ∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣ s′=s ∥∥∥∥ 2 . (8)
This gives us the novel lower bound for payoff improvement:
Lsupπ1 (π2) = ηπ1 + 1
1− γ Es∼ρπ1µ [A π1(s, π2)]− C1 · sup s
[ LWG(π1, π2; s) ] − C2 · sup
s
[ W (π2(a|s), π1(a|s))2 ] . (9)
We have the inequalities ηπ2 ≥ Lsupπ1 (π2) and L sup π1 (π1) = ηπ1 . This facilitates the application of Theorem 2.1 as an instance of Minorization-Maximization algorithms (Hunter & Lange, 2004). Theorem 2.1. Successive maximization of Lsup yields non-decreasing policy payoffs.
Proof. With π2 = arg maxπ Lsupπ1 (π), we have L sup π1 (π2) ≥ L sup π1 (π1). Thus,
ηπ2 ≥ Lsupπ1 (π2) and ηπ1 = L sup π1 (π1) =⇒ ηπ2 − ηπ1 ≥ L sup π1 (π2)− L sup π1 (π1) ≥ 0. (10)

heading: 3 SURROGATE OPTIMIZATION AND PRACTICAL ALGORITHM
text: Successive optimization of Lsupπ1 (π2) generates non-decreasing payoffs. However, it is impractical due to the large number of constraints and statistical estimation of maximums. To mitigate this, we take a similar approach to TRPO and optimize for the surrogate
L̄π1(π2) = ηπ1 + 1
1− γ · Es∼ρπ1µ [A π1(s, π2)]− C1 · Es∼ρπ1µ
[ LWG(π1, π2; s) ] − C2 · Es∼ρπ1µ [ W (π2(a|s), π1(a|s))2 ] . (11)
Although first order stochastic optimization methods can be directly applied to the surrogate defined in Equation (11), second order methods can be more efficient. Since LWG(π1, π2; s) is the geometric mean of two functions of quadratic order, it is also of quadratic order. However, LWG(π1, π2; s) may not be twice continuously differentiable. For this, we lower bound LWG(π1, π2; s) further using the AM-GM inequality and replace it with
LG2(π1, π2; s) := ∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣ s′=s ∥∥∥∥2 2
(12)
to form a quadratic regularization term with a definable Hessian matrix (see Section A.8 of the appendix for detailed derivations). This induces our final surrogate which is used for second order optimization:
Lπ1(π2) = ηπ1 + 1
1− γ · Es∼ρπ1µ [A π1(s, π2)]− C ′1 · Es∼ρπ1µ
[ LG2(π1, π2; s) ] − C ′2 · Es∼ρπ1µ [ W (π2(a|s), π1(a|s))2 ] . (13)
The coefficients C1 and C2 are dynamics-dependent. For simplicity, we used constant coefficients and a trust region. This yields the Truly Deterministic Policy Optimization (TDPO) as given in Algorithm 1. See the appendix for practical notes on the choice of C1 and C2. Alternatively, one could adopt processes similar to Schulman et al. (2015a) where the IS-based advantage estimator used a line search for proper step size selection, or the adaptive penalty coefficient setting in Schulman et al. (2017). We plan to consider such approaches in future work.

heading: 3.1 ON THE INTERPRETATION OF THE SURROGATE FUNCTION
text: For deterministic policies, the squared Wasserstein distance W (π2(a|s), π1(a|s))2 degenerates to the Euclidean distance over the action space. Any policy defines a sensitivity matrix at a given state s, which is the Jacobian matrix of the policy output with respect to s. The policy sensitivity
Algorithm 1 Truly Deterministic Policy Optimization (TDPO) Require: The squared Wasserestein regularization coefficient C ′2, The secondary regularization
coefficient ratio C ′1/C ′ 2, and a trust region radius δmax.
Require: Initial policy π0. Require: Advantage estimator and sample collector oracle Aπ .
1: for k = 1, 2, . . . do 2: Collect trajectories and construct the advantage estimator oracle Aπk . 3: Compute the policy gradient g at θk : g ← ∇θ′Aπk(π′)|π′=πk 4: Construct a surrogate Hessian vector product oracle v → H · v such that for θ′ = θk + δθ,
Es∼ρπkµ
[ W (π′(a|s), πk(a|s))2 ] + C ′1 C ′2 Es∼ρπkµ [ LG2(π′, πk; s) ] = 1 2 δθTHδθ+ h.o.t., (14)
where h.o.t. denotes higher order terms in δθ. 5: Find the optimal update direction δθ∗ = H−1g using the Conjugate Gradient algorithm. 6: Determine the best step size α∗ within the trust region:
α∗ = arg max α
gT (αδθ∗)− C ′ 2
2 (αδθ∗)TH(αδθ∗)
s.t. 1
2 (α∗δθ∗)TH(α∗δθ∗) ≤ δ2max (15)
7: Update the policy parameters: θk+1 ← θk + α∗δθ∗. 8: end for
term LG2(π1, π2; s) is essentially the squared Euclidean distance over the action-to-observation Jacobian matrix elements. In other words, our surrogate prefers to step in directions where the action-to-observation sensitivity is preserved within updates.
Although our surrogate uses a metric distance instead of the traditional non-metric measures for regularization, we do not consider this sole replacement a major contribution. The squared Wasserestein distance and the KL divergence of two identically-scaled Gaussian distributions are the same up to a constant (i.e., DKL(N (m1, σ)‖N (m2, σ)) = W (N (m1, σ),N (m2, σ))2/2σ2). On the other hand, our surrogate’s compatibility with deterministic policies makes it a valuable asset for our policy gradient algorithm; both W (π2(a|s), π1(a|s))2 and LG2(π1, π2; s) can be evaluated for two deterministic policies π1 and π2 numerically without resorting to any approximations to overcome singularities.

heading: 4 MODEL-FREE ESTIMATION OF POLICY GRADIENT
text: The DeVine advantage estimator is formally defined in Algorithm 2. Unlike DDPG and TD3, the DeVine estimator allows our method to perform deterministic search by not consistently injecting noise in actions for exploration. Essentially, DeVine rolls out a trajectory and computes the values of each state. Since the transition dynamics and the policy are deterministic, these values are exact. Then, it picks a perturbation state st according to the visitation frequencies using importance sampling. A state-reset to st is made, a σ-perturbed action is applied for a single time-step, followed by π1 policy. This exactly produces Qπ1(st, at + σ). Then, Aπ1(st, at + σ) can be computed by subtracting the value baseline. Finally, Aπ1(st, at) = 0 and Aπ1(st, at + σ) define a two-point linear Aπ1(st, a) model with respect to the action. Parallelization can be used to have as many states of the first roll-out included in the estimator as desired. The parameter σ acts as an exploration parameter and a finite difference to establish derivatives. While σ ' 0 can produce exact gradients, larger σ can build stabler interpolations.
Under deterministic dynamics and policies, if the DeVine oracle samples each dimension at each time-step exactly once then in the limit of small σ it can produce exact advantages, as stated in Theorem 4.1, whose proof is deferred to the appendix. Theorem 4.1. Assume a finite horizon MDP with both deterministic transition dynamics P and initial distribution µ, with maximal horizon length of H . Define K = H · dim(A), a uniform ν, and
Algorithm 2 Deterministic Vine (DeVine) Policy Advantage Estimator Require: The number of parallel workers K Require: A policy π, an exploration policy q, discrete time-step distribution ν(t), initial state
distribution µ(s), and the discount factor γ. 1: Sample an initial state s0 from µ, and then roll out a trajectory τ = (s0, a0, s1, a1, · · · ) using π. 2: for k = 1, 2, · · · ,K do 3: Sample the integer number t = tk from ν. 4: Compute the value V π1(st) = ∑∞ i=t γ
t−iR(si, ai). 5: Reset the initial state to st, sample the first action a′t according to q(·|st), and use π for the
rest of the trajectory. This will create τ ′ = (st, a′t, s ′ t+1, a ′ t+1, · · · ). 6: Compute the value Qπ1(st, a′t) = ∑∞ i=t γ t−iR(s′i, a ′ i). 7: Compute the advantage Aπ1(st, a′t) = Q π(st, a ′ t)− V π(st). 8: end for
9: Define Aπ1(π2) := 1
K K∑ k=1 dim(A) · γtk ν(tk) · (π2(s)− atk)T (a′tk − atk) (a′tk − atk)T (a ′ tk − atk) ·Aπ1(stk , a′tk).
10: Return Aπ1(π2) and ∇π2Aπ1(π2) as unbiased estimators for Es∼ρπ1µ [A π1(s, π2)] and PG, re-
spectively.
q(s;σ) = π1(s) + σej in Algorithm 2 with ej being the jth basis element for A. If the (j, tk) pairs are sampled to exactly cover {1, · · · ,dim(A)} × {1, · · · , H}, then we have
lim σ→0 ∇π2Aπ1(π2) ∣∣ π2=π1 = ∇π2ηπ2 ∣∣ π2=π1 . (16)
Theorem 4.1 provides a guarantee for recovering the exact policy gradient if the initial state distribution was deterministic and all time-steps of the trajectory were used to branch vine trajectories. Although this theorem sets the stage for computing a fully deterministic gradient, stochastic approximation can be used in Algorithm 2 by randomly sampling a small set of states for advantage estimation. In other words, Theorem 4.1 would use ν to deterministically sample all trajectory states, whereas this is not a practical requirement for Algorithm 2 and the gradients are still unbiased if a random set of vine branches is used.
The DeVine estimator can be advantageous in at least two scenarios. First, in the case of rewards that cannot be decomposed into summations of immediate rewards. For example, overshoot penalizations or frequency-based rewards as used in robotic systems are non-local. DeVine can be robust to non-local rewards as it is insensitive to whether the rewards were applied immediately or after a long period. Second, DeVine can be an appropriate choice for systems that are sensitive to the injection of noise, such as high-bandwidth robots with natural resonant frequencies. In such cases, using white (or colored) noise for exploration can excite these resonant frequencies and cause instability, making learning difficult. DeVine avoids the need for constant noise injection.

heading: 5 EXPERIMENTS
text: The next two subsections show challenging robotic control tasks including frequency-based non-local rewards, long horizons, and sensitivity to resonant frequencies. See Section A.11 of the appendix for a comparison on traditional gym environments. TDPO works similar or slightly worse on these traditional gym environments as they seem to be well-suited for stochastic exploration.

heading: 5.1 AN ENVIRONMENT WITH NON-LOCAL REWARDS 1
text: The first environment that we consider is a simple pendulum. The transition function is standard—the states are joint angle and joint velocity, and the action is joint torque. The reward function is nonstandard—rather than define a local reward in the time domain with the goal of making the pendulum
1Non-local rewards are reward functions of the entire trajectory whose payoffs cannot be decomposed into the sum of terms such as η = ∑ t ft(st, at), where functions ft only depend on nearby states and actions. An example non-local reward is one that depends on the Fourier transform of the complete trajectory signal.
stand upright (for example), we define a non-local reward in the frequency domain with the goal of making the pendulum oscillate with a desired frequency and amplitude about a desired offset. In particular, we compute this non-local reward by taking the Fourier transform of the joint angle signal over the entire trajectory and by penalizing differences between the resulting power spectrum and a desired power spectrum. We apply this non-local reward at the last time step of the trajectory. Implementation details and similar results for more pendulum variants are left to the appendix.
Figure 1 shows training curves for TDPO (our method) as compared to TRPO, PPO, DDPG, and TD3. These results were averaged over 25 experiments in which the desired oscillation frequency was 1.7 Hz (different from the pendulum’s natural frequency of 0.5 Hz), the desired oscillation amplitude was 0.28 rad, and the desired offset was 0.52 rad. Figure 1 also shows trajectories obtained by the best agents from each method. TDPO (our method) was able to learn high-reward behavior and to achieve the desired frequency, amplitude, and offset. TRPO was able to learn the correct offset but did not produce any oscillatory behavior. TD3 also learned the correct offset, but could not produce desirable oscillations. PPO and DDPG failed to learn any desired behavior.

heading: 5.2 AN ENVIRONMENT WITH LONG HORIZON AND RESONANT FREQUENCIES2
text: The second environment that we consider is a single leg from a quadruped robot (Park et al., 2017). This leg has two joints, a “hip” and a “knee,” about which it is possible to exert torques. The hip is attached to a slider that confines motion to a vertical line above flat ground. We assume the leg is dropped from some height above the ground and the task is to recover from this drop and to stand upright at rest after impact. States given to the agent are the angle and velocity of each joint (slider position and velocity are hidden), and actions are the joint torques. The reward function penalizes difference from an upright posture, slipping or chattering at the contact between the foot and the
2Resonant frequencies are a concept from control theory. In the frequency domain, signals of certain frequencies are excited more than others when applied to a system. This is captured by the frequency-domain transfer function of the system, which may have a peak of magnitude greater than one. The resonant frequency is the frequency at which the frequency-domain transfer function has the highest amplitude. Common examples of systems with a resonant frequency include the undamped pendulum, which oscillates at its natural frequency, and RLC circuits which have characteristic frequencies at which they are most excitable. See Chapter 8 of Kuo & Golnaraghi (2002) for more information.
ground, non-zero joint velocities, and steady-state joint torque deviations. We use MuJoCo for simulation (Todorov et al., 2012), with high-fidelity models of ground compliance, motor nonlinearity, and joint friction. The control loop rate is 4000 Hz and the rollout length is 2 s, resulting in a horizon of 8000 steps. Implementation details are left to the appendix.
Figure 2 shows training curves for TDPO (our method) as compared to TRPO, PPO, DDPG and TD3. These results were averaged over 75 experiments. A discount factor of γ = 0.99975 was chosen for all methods, where (1− γ)−1 is half the trajectory length. Similarly, the GAE factors for PPO and TRPO were scaled up to 0.99875 and 0.9995, respectively, in proportion to the trajectory length. Figure 2 also shows trajectories obtained by the best agents from each method. TDPO (our method) was able to learn high-reward behavior. TRPO, PPO, DDPG, and TD3 were not.
We hypothesize that the reason for this difference in performance is that TDPO better handles the combination of two challenges presented by the leg environment—an unusually long time horizon (8000 steps) and the existence of a resonant frequency that results from compliance between the foot and the ground (note the oscillations at a frequency of about 100 Hz that appear in the trajectories after impact). Both high-speed control loops and resonance due to ground compliance are common features of real-world legged robots to which TDPO seems to be more resilient.

heading: 6 DISCUSSION
text: This article proposed a deterministic policy gradient method (TDPO: Truly Deterministic Policy Optimization) based on the use of a deterministic Vine (DeVine) gradient estimator and the Wasserstein metric. We proved monotonic payoff guarantees for our method, and proposed a novel surrogate for policy optimization. We showed numerical evidence for superior performance with non-local rewards defined in the frequency domain and a realistic long-horizon resonant environment. This method enables applications of policy gradient to customize frequency response characteristics of agents. Future work should include the addition of a regularization coefficient adaptation process for the C1 and C2 parameters in the TDPO algorithm.

heading: A APPENDIX

heading: A.1 TABLES OF NOTATION
text: The same mathematical definitions and notations used in the paper were re-introduced and summarized in two tables; Table 1 describes the mathematical functions and operators used throughout the paper, and Table 2 describes the notations needed to define the Markov Decision Process (MDP). The tables consist of two columns; one showing or defining the notation, and the other includes the name in which the same notation was called in the paper.
A.2 BOUNDING W (Pt(µ, π1),Pt(µ, π2))
To review, the dynamical smoothness assumptions were
W (P(µ, π1),P(µ, π2)) ≤ Lπ ·W (π1, π2), W (P(µ1, π),P(µ2, π)) ≤ Lµ ·W (µ1, µ2).
The following lemma states that these two assumptions are equivalent to a more concise assumption. This will be used to bound the t-step visitation distance and prove Lemma A.2. Lemma A.1. Assumptions (5) and (6) are equivalent to having
W (P(µ1, π1),P(µ2, π2)) ≤ Lµ ·W (µ1, µ2) + Lπ ·W (π1, π2). (17)
Proof. To prove the (5), (6)⇒ (17) direction, the triangle inequality for the Wasserstein distance gives
W (P(µ1, π1),P(µ2, π2)) ≤W (P(µ1, π1),P(µ2, π1)) +W (P(µ2, π1),P(µ2, π2)) (18)
and using (5), (6), and (18) then implies
W (P(µ1, π1),P(µ2, π2)) ≤ Lµ ·W (µ1, µ2) + Lπ ·W (π1, π2). (19)
The other direction is trivial.
Lemma A.2. Under Assumptions (5) and (6) we have the bound
W (Pt(µ, π1),Pt(µ, π2)) ≤ Lπ · (1 + Lµ + · · ·+ Lt−1µ ) ·W (π1, π2), (20)
where Pt(µ, π) denotes the state distribution after running the MDP for t time-steps with the initial state distribution µ and policy π.
Proof. For t = 1, the lemma is equivalent to Assumption (5). This paves the way for the lemma to be proved using induction. The hypothesis is
W (Pt−1(µ, π1),Pt−1(µ, π2)) ≤ Lπ · (1 + Lµ + · · ·+ Lt−2µ ) ·W (π1, π2), (21)
and for the induction step we write
W (Pt(µ, π1),Pt(µ, π2)) = W (P(Pt−1(µ, π1), π1),P(Pt−1(µ, π2), π2)). (22)
Using Assumption (17), which according to Lemma A.1 is equivalent to Assumptions (5) and (6), we can combine (21) and (22) into
W (Pt(µ, π1),Pt(µ, π2)) ≤ Lπ ·W (π1, π2) + Lµ ·W (Pt−1(µ, π1),Pt−1(µ, π2)). (23)
Thus, by applying the induction Hypothesis (21), we have
W (Pt(µ, π1),Pt(µ, π2)) ≤ Lπ ·W (π1, π2) +Lµ ·Lπ · (1 +Lµ + · · ·+Lt−2µ ) ·W (π1, π2), (24)
which can be simplified into the lemma statement (i.e., Inequality (20)).
A.3 BOUNDING W (ρπ1µ , ρ π2 µ )
Lemma A.2 suggests making the γLµ < 1 assumption and paves the way for Theorem A.4. The γLµ < 1 assumption is overly restrictive and unnecessary, but makes the rest of the proof easier to follow. This assumption can be relaxed by a general transition dynamics stability assumption which is discussed in more detail later at section A.5.3, and an equivalent γL̄µ < 1 assumption is introduced to replace γLµ < 1.
First, we need to introduce Lemma A.3 first, which will be used in the proof of Theorem A.4. Lemma A.3. The Wasserstein distance between linear combinations of distributions can be bounded as W (β · µ1 + (1− β) · ν1, β · µ2 + (1− β) · ν2) ≤ β ·W (µ1, µ2) + (1− β) ·W (ν1, ν2).
Proof. Plugging γ = β ·γ(µ1,µ2) +(1−β) ·γ(ν1,ν2) in the Wasserstein definition yields the result.
Theorem A.4. Assuming (5), (6), and γLµ < 1, we have the inequality
W (ρπ1µ , ρ π2 µ ) ≤ γLπ 1− γLµ ·W (π1, π2). (25)
Proof. Using Lemma A.3 and the definition of ρπµ, we can write
W (ρπ1µ , ρ π2 µ ) ≤ (1− γ) ∞∑ t=0 γt ·W (Pt(µ, π1),Pt(µ, π2)). (26)
Using Lemma A.2, we can take another step to relax the inequality (26) and write
W (ρπ1µ , ρ π2 µ ) ≤ Lπ(1− γ)W (π1, π2) (Lµ − 1) ∞∑ t=0 ((γLµ) t − γt). (27)
Due to the γLµ < 1 assumption, the right-hand summation in (27) is convergent. This leads us to
W (ρπ1µ , ρ π2 µ ) ≤ Lπ(1− γ)W (π1, π2) (Lµ − 1) ( 1 1− γLµ − 1 1− γ ). (28)
Inequality (28) can then be simplified to give the result.

heading: A.4 STEPS TO BOUND THE SECOND-ORDER TERM
text: The RK duality yields the following bound:
|〈ρπ2µ − ρπ1µ , Aπ1(·, π2)〉s| ≤W (ρ π1 µ , ρ π2 µ ) · sup s ‖∇sAπ1(s, π2)‖2. (29)
To facilitate the further application of the RK duality, any advantage can be rewritten as the following inner product: Aπ1(s, π2) = 〈π2(a|s)− π1(a|s), Qπ1(s, a)〉a. Taking derivatives of both sides with respect to the state variable and applying the triangle inequality produces the bound
sup s ‖∇sAπ1(s, π2)‖2 ≤ sup s ‖〈∇s(π2(a|s)− π1(a|s)), Qπ1(s, a)〉a‖2
+ sup s ‖〈π2(a|s)− π1(a|s),∇sQπ1(s, a)〉a‖2. (30)
The second term of the RHS in (30) is compatible with the RK duality. However, the form of the first term does not warrant an easy application of RK. For this, we introduce Theorem A.5. Theorem A.5. Assuming the existence of Lip(Qπ1(s, a); a), we have the bound∥∥∥∥〈∇s(π2(a|s)− π1(a|s)), Qπ1(s, a)〉a∥∥∥∥
2
(31)
≤ 2 · Lip(Qπ1(s, a); a) · ∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|bs′)2 )∣∣∣∣ s′=s ∥∥∥∥ 2 .
Proof. By definition, we have∥∥〈∇s(π2(a|s)− π1(a|s)), Qπ1(s, a)〉a∥∥2 =
√√√√dim(S)∑ j=1 (〈 ∂ ∂s(j) (π2(a|s)− π1(a|s)), Qπ1(s, a) 〉 a )2 . (32)
For better insight, we will write the derivative using finite differences:〈 ∂ ∂s(j) (π2(a|s)−π1(a|s)), Qπ1(s, a) 〉 a
= lim δs→0
1
δs
[〈 (π2(a|s+ δs · ej) −π1(a|s+ δs · ej) ), Qπ1(s, a) 〉 a
− 〈 (π2(a|s) −π1(a|s) ), Qπ1(s, a) 〉 a ] . (33)
We can rearrange the finite difference terms to get〈 ∂ ∂s(j) (π2(a|s)−π1(a|s)), Qπ1(s, a) 〉 a
= lim δs→0
1
δs
[〈 (π2(a|s+ δs · ej) +π1(a|s) ), Qπ1(s, a) 〉 a
− 〈 (π2(a|s) +π1(a|s+ δs · ej) ), Qπ1(s, a) 〉 a ] . (34)
Equivalently, we can divide and multiply the inner products by a factor of 2, to make the inner product arguments resemble mixture distributions:〈 ∂
∂s(j) (π2(a|s)− π1(a|s)), Qπ1(s, a) 〉 a
= lim δs→0
2
δs [〈π2(a|s+ δs · ej) + π1(a|s) 2 , Qπ1(s, a) 〉 a
− 〈π2(a|s) + π1(a|s+ δs · ej)
2 , Qπ1(s, a) 〉 a ] . (35)
The RK duality can now be used to bound this difference as∣∣∣∣〈 ∂∂s(j) (π2(a|s)− π1(a|s)), Qπ1(s, a)〉a ∣∣∣∣ (36)
≤ lim δs→0
2
δs
[ W ( π2(a|s+ δs · ej) + π1(a|s)
2 , π2(a|s) + π1(a|s+ δs · ej) 2
) · Lip(Qπ1(s, a); a) ] ,
which can be simplified as∣∣∣∣〈 ∂∂s(j) (π2(a|s)− π1(a|s)), Qπ1(s, a)〉a ∣∣∣∣
≤ 2 · Lip(Qπ1(s, a); a) · ∂ ∂s′(j) W
( π2(a|s′) + π1(a|s)
2 , π2(a|s) + π1(a|s′) 2 )∣∣∣∣ s′=s . (37)
Combining Inequality (37) and Equation (32), we obtain the bound in the theorem.

heading: A.5 THE LIPSCHITZ CONTINUITY AND TRANSITION STABILITY ASSUMPTIONS
text: There are three key groups of assumptions made in the derivation of our policy improvement lower bound. First is the existence of Qπ-function Lipschitz constants. Second is the transition dynamics Lipschitz-continuity assumptions. Finally, we make an assumption about the stability of the transition dynamics. Next, we will discuss the meaning and the necessity of these assumptions in the same order.
A.5.1 ON THE EXISTENCE OF THE LIP(Qπ, a) CONSTANT
The Lip(Qπ, a) constant may be undefined when either the reward function or the transition dynamics are discontinuous. Examples of known environments with undefined Lip(Qπ, a) constants include those with grazing contacts which define a discontinuous transition dynamics. In practice, even for environments that do not satisfy Lipschitz continuity assumptions, there are mitigating factors; practical Qπ functions are reasonably narrow-bounded in a small trust region neighborhood, and since we use non-vanishing exploration scales and trust regions, a bounded interpolation slope can still model the Q-function variation effectively. We should also note that a slightly stronger version of this assumption is frequently used in the context of Lipschitz MDPs (Pirotta et al., 2015; Rachelson & Lagoudakis, 2010; Asadi et al., 2018). In practice, we have not found this to be a substantial limitation.

heading: A.5.2 THE TRANSITION DYNAMICS LIPSCHITZ CONTINUITY ASSUMPTION
text: Assumptions 5 and 6 of the main paper essentially represent the Lipschitz continuity assumptions of the transition dynamics with respect to actions and states, respectively. If the transition dynamics
and the policy are deterministic, then these assumptions are exactly equivalent to the Lipschitz continuity assumptions. Assumptions 5 and 6 only generalize the Lipschitz continuity assumptions in a distributional sense.
The necessity of these assumptions is a consequence of using metric measures for bounding errors. Traditional non-metric bounds force the use of full-support stochastic policies where all actions have non-zero probabilities (e.g., for the KL-divergence of two policies to be defined, TRPO needs to operate on full-support policies such as the Gaussian policies). In those analyses, since all policies share the same support, the next state distribution automatically becomes smooth and Lipschitz continuous with respect to the policy measure even if the transition dynamics was not originally smooth with respect to its input actions. However, metric measures are also defined for policies of non-overlapping support. To be able to provide closeness bounds for future state visitations of two similar policies with non-overlapping support, it becomes necessary to assume that close-enough actions or states must be yielding close-enough next states. In fact, this is a very common assumption in the framework of Lipschitz MDPs (See Section 2.2 of Rachelson & Lagoudakis (2010), Section 3 of Asadi et al. (2018), and Assumption 1 of Pirotta et al. (2015)).

heading: A.5.3 THE TRANSITION DYNAMICS STABILITY ASSUMPTION
text: Before moving to relax the γLµ < 1 assumption, we will make a few definitions and restate the previous lemmas and theorems under these definitions. We define Lµ1,µ2,π to be the infimum non-negative value that makes the equation W (P(µ1, π),P(µ2, π)) = Lµ1,µ2,πW (µ1, µ2) hold. Similarly, Lµ1,µ2,π is defined as the infimum non-negative value that makes the equation W (P(µ, π1),P(µ, π)) = Lµ,π1,π2W (π1(·|µ), π2(·|µ)) hold. For notation brevity, we will also denote LPt(µ,π1),Pt(µ,π2),π2 and LPt(µ,π1),π1,π2 by L̃ (t) µ,π1,π2 and L̂ (t) µ,π1,π2 , respectively.
Under these definitions, Lemma A.1 evolves into
W (P(µ1, π1),P(µ2, π2)) ≤ Lµ1,µ2,πW (µ1, µ2) + Lµ1,π1,π2W (π1, π2). (38) We can apply a time-point recursion to this lemma and have
W (P(Pt(µ, π1), π1),P(Pt(µ, π2), π2)) ≤ LPt(µ,π1),π1,π2W (π1, π2) + LPt(µ,π1),Pt(µ,π2),π2W (P
t(µ, π1),Pt(µ, π2)) (39) , which can be notationally simplified to
W (Pt(µ, π1),Pt(µ, π2)) ≤ L̂(t−1)µ,π1,π2W (π1, π2) + L̃ (t−1) µ,π1,π2W (P t−1(µ, π1),Pt−1(µ, π2)). (40)
These modifications lead Lemma A.2 to be updated accordingly into
W (Pt(µ, π1),Pt(µ, π2)) ≤ C(t)L;µ,π1,π2 ·W (π1, π2) (41) , where we have
C (t) L;µ,π1,π2 := t∑ k=1 L̂(t−k)µ,π1,π2 k−1∏ i=1 L̃(t−i)µ,π1,π2 . (42)
By a simple change of variables, we can have the equivalent definition of
C (t) L;µ,π1,π2 := t∑ k=1 L̂(k−1)µ,π1,π2 t−1∏ i=k+1 L̃(i)µ,π1,π2 . (43)
Now, we would replace the γLµ < 1 assumption with the following assumption.
The Transition Dynamics Stability Assumption: A transition dynamics P is called stable if and only if the induced {L̃(t)µ,π1,π2}t≥0 and {L̂ (t) µ,π1,π2}t≥0 sequences satisfy
CL := sup µ,π1,π2,t
C (t) L;µ,π1,π2
= sup µ,π1,π2,t t∑ k=1 L̂(k−1)µ,π1,π2 t−1∏ i=k+1 L̃(i)µ,π1,π2 <∞. (44)
To help understand which {L̃(t)µ,π1,π2}t≥0 and {L̂ (t) µ,π1,π2}t≥0 sequences can satisfy this assumption, we will provide some examples:
• The ∀t : L̃(t)µ,π1,π2 = c1 > 1 and ∀t : L̂ (t) µ,π1,π2 = c2 sequences violate the dynamics stability
assumption.
• The ∀t : L̃(t)µ,π1,π2 ≤ 1 and ∀t : L̂ (t) µ,π1,π2 = O( 1 t2 ) sequences satisfy the dynamics stability
assumption.
• supt L̃ (t) µ,π1,π2 < 1 guarantees the dynamics stability assumption.
• ∀t ≥ t0 : L̃(t)µ,π1,π2 < 1 guarantees the dynamics stability assumption no matter (1) how big t0 is (as long as it is finite), or (2) how big the members of the finite set {L̃(t)µ,π1,π2 |t < t0} are.
If the dynamics stability assumption holds with a constant CL, one can define a L̄µ constant such that CL = Lπ ∑∞ t=0(γL̄µ)
t. Then, we can replace all the Lµ instances in the rest of the proof with the corresponding L̄µ constant, and the results will remain the same without any change of format.
The L̃(t)µ,π1,π2 and L̂ (t) µ,π1,π2 constants can be thought as tighter versions of Lµ and Lπ, but with dependency on π1, π2, µ and the time-point of application. Having γLµ < 1 is a sufficient yet unnecessary condition for this dynamics stability assumption to hold. Vaguely speaking, Lµ is an expansion rate for the state distribution distance; it tells you how much a divergence in the state distribution will expand after a single application of the transition dynamics. Having effective expansion rates that are larger than one throughout an infinite horizon trajectory is a sign of the system instability; some change in the initial state’s distribution could cause the observations to diverge exponentially. While controlling unstable systems is an important and practical challenge, none of the existing reinforcement learning methods is capable of learning effective policies on such environments. Roughly speaking, having the dynamics stability assumption guarantees that the expansion rates cannot be consistently larger than one for infinite time steps.

heading: A.6 CHOICE OF C1 AND C2
text: Since the TDPO algorithm operates using the metric Wasserstein distance, thinking about how normalizing actions and rewards affect the corresponding optimization objective builds insight into how to set these coefficients properly. Say we use the same dynamics, only the new actions are scaled up by a factor of β and the rewards are scaled up by a factor of α:
anew = β · aold rnew = α · rold. (45)
If the policy function approximation class remained the same, the policy gradient would be scaled by a factor of αβ (i.e., ∂ηnew ∂anew = αβ · ∂ηold ∂aold
). Therefore, one can easily show that the corresponding new regularization coefficient and trust region sizes can be obtained by
Cnew = α
β2 · Cold (46)
and δnewmax = β · δoldmax. (47)
We used equal regularization coefficients (i.e., C1 = C2 = C), and the process to choose them can be summarized as follows: (1) Define C = 3600 · α · β−2, δmax = β/600 and σq = β/60 (where σq is the action disturbance parameter used for DeVine), (2) using prior knowledge or by trial and error determine appropriate action and reward normalization coefficients. The reward normalization coefficient α was sought to be approximately the average per-step discounted reward difference of a null policy and an optimal policy. We used a reward scaling value of α = 5 and an action scaling value of β = 5 for the non-locally rewarded pendulum and β = 1.5 for the long-horizon legged robot. Both environments had a per-step discounted reward of approximately −5 for a null policy and non-positive rewards, justifying the choice of α.

heading: A.7 PROOF OF THEOREM 4.1
text: We restate Theorem 4.1 below for reference and now prove it.
Theorem 4.1. Assume a finite horizon MDP with both deterministic transition dynamics P and initial distribution µ, with maximal horizon length of H . Define K = H · dim(A), a uniform ν, and q(s;σ) = π1(s) + σej in Algorithm 2 with ej being the jth basis element for A. If the (j, tk) pairs are sampled to exactly cover {1, . . . ,dim(A)} × {1, . . . ,H}, then we have
lim σ→0 ∇π2Aπ1(π2) ∣∣ π2=π1 = ∇π2ηπ2 ∣∣ π2=π1 . (48)
Proof. According to the advantage decomposition lemma, we have ∇π2ηπ2 ∣∣ π2=π1 = 1 1− γ Es∼ρπ1µ [∇π2A π1(s, π2)] ∣∣ π2=π1 . (49)
Due to the fact that the transition dynamics, policies π1 and π2, and initial state distribution are all deterministic, we can simplify Equation (49) to
∇π2ηπ2 ∣∣ π2=π1 = H−1∑ t=0 γt · ∇π2Aπ1(st, π2) ∣∣ π2=π1 , (50)
where st is the state after applying the policy π1 for t time-steps. We can use the chain rule to write ∇π2Aπ1(st, π2) ∣∣ π2=π1 = ∇π2Aπ1(st, at) ∣∣ at=π2(st) π2=π1
= dim(A)∑ j=1 ∇π2a (j) t · ∂ ∂a (j) t Aπ1(st, at) ∣∣ at=π2(st) π2=π1 . (51)
To recap, Equations (50), (50), and (51) can be summarized as
∇π2ηπ2 ∣∣ π2=π1 = H−1∑ t=0 γt dim(A)∑ j=1 ∇π2a (j) t · ∂ ∂a (j) t Aπ1(st, at) ∣∣ at=π2(st) π2=π1 . (52)
Under the assumption that the (j, t) pairs are sampled to exactly cover {1, . . . ,dim(A)}×{1, . . . ,H}, we can simplify the DeVine oracle to
Aπ1(π2) = 1
K H−1∑ t=0 dim(A)∑ j=1 [ dim(A) · γt ν(t) · (π2(st)− π1(st)) T (q(st; j, σ)− π1(st)) (q(st; j, σ)− π1(st))T (q(st; j, σ)− π1(st))
·Aπ1(st, q(st; j, σ)) ] . (53)
From the q definition, we have q(st; j, σ)− π1(st) = σej and (q(st; j, σ)− π1(st))T (q(st; j, σ)− π1(st)) = σ
2. Since ν is uniform (i.e., ν(t) = 1H ) and K = H · dim(A), we can take the policy gradient of Equation (53) and simplify it into
∇π2Aπ1(π2) ∣∣ π2=π1 = H−1∑ t=0 dim(A)∑ j=1 [ γt · ∇π2(π2(st)− π1(st))Tej · Aπ1(st, q(st; j, σ)) σ ] . (54)
Since, Aπ1(st, π1(st)) = 0, we can write
lim σ→0
Aπ1(st, q(st; j, σ))
σ = lim σ→0 Aπ1(st, π1(st) + σej)−Aπ1(st, π1(st)) σ
= ∂
∂a (j) t
Aπ1(st, at) ∣∣ at=π1(st) . (55)
Also, by the definition of the gradient, we can write
∇π2(π2(st)− π1(st))Tej = ∇π2a (j) t . (56)
Combining Equations (55) and (56), and applying them to Equation (54), yields
lim σ→0 ∇π2Aπ1(π2) ∣∣ π2=π1 = H−1∑ t=0 dim(A)∑ j=1 γt · ∇π2a (j) t · ∂ ∂a (j) t Aπ1(st, at) ∣∣ at=π2(st) π2=π1 . (57)
Finally, the theorem can be obtained by comparing Equations (57) and (52).

heading: A.8 QUADRATIC MODELING OF POLICY SENSITIVITY REGULARIZATION
text: First, we will build insight into the nature of the LWG(π1, π2; s) =W (π2(a|s), π1(a|s))×∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣ s′=s ∥∥∥∥ 2
(58)
term. It is fairly obvious that W (π2(a|s), π1(a|s)) ∣∣ π2=π1 = 0. (59) If π2 = π1, then the two distributions π2(a|s′)+π1(a|s) 2 and π2(a|s)+π1(a|s′)
2 will be the same no matter what s′ is. In other words,
π1 = π2 ⇒ ∀s′ : W ( π2(a|s′) + π1(a|s)
2 , π2(a|s) + π1(a|s′) 2
) = 0. (60)
This means that∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣
s′=s ∥∥∥∥ 2 ∣∣∣∣ π2=π1 = 0. (61)
The Taylor expansion of the squared Wasserestein distance can be written as W (π2(a|s), π1(a|s))2 ∣∣ θ2=θ1+δθ = 1 2 δθTH2δθ + h.o.t.. (62)
Considering (60) and similar to the previous point, one can write the following Taylor expansion∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣
s′=s ∥∥∥∥2 2 ∣∣∣∣ θ2=θ1+δθ = δθTH1δθ + h.o.t.. (63)
According to above, LWG is the geometric mean of two functions of quadratic order. Although this makes LWG of quadratic order (i.e., limδθ→0 LWG(αδθ)LWG(δθ) = α
2 holds for any constant α), this does not guarantee that LWG is twice continuously differentiable w.r.t. the policy parameters, and may not have a defined Hessian matrix (e.g., f(x1, x2) = |x1x2| is of quadratic order, yet is not twice differentiable). To avoid this issue, we compromise on the local model. Using the AM-GM inequality and for any arbitrary positive α, one can bound the LWG term into two quadratic terms:
LWG(π1, π2; s) ≤ 1
2
( 1
α2 W (π2(a|s), π1(a|s))2+
α2 ∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣ s′=s ∥∥∥∥2 2 ) . (64)
Therefore, by defining LG2(π1, π2; s) := ∥∥∥∥∇s′W(π2(a|s′) + π1(a|s)2 , π2(a|s) + π1(a|s′)2 )∣∣∣∣ s′=s ∥∥∥∥2 2 , (65) C ′1 := C1·α2 2 , and C ′ 2 := (C2 + C1 2α2 ) the new surrogate will have the twice-differentiable form
Lπ1(π2) = 1
1− γ · Es∼ρπ1µ [A π1(s, π2)]
− C ′1 · Es∼ρπ1µ [ LG2(π1, π2; s) ] − C ′2 · Es∼ρπ1µ [ W (π2(a|s), π1(a|s))2 ] .
(66)
C ′1 and C ′ 2 will be the corresponding regularization coefficients for the surrogate defined in (13). Due to the arbitrary α used in bounding, no constrain governs the C ′1 and C ′ 2 coefficients. Therefore, C ′ 1 and C ′2 can be chosen without constraining each other.
A.9 IMPLEMENTATION DETAILS FOR THE ENVIRONMENT WITH NON-LOCAL REWARDS
We used the stable-baselines implementation (Hill et al., 2018), which has the same structure as the original OpenAI baselines (Dhariwal et al., 2017) implementation. We used the “ppo1” variant since no hardware acceleration was necessary for automatic differentiation and MPI parallelization was practically efficient. TDPO, TRPO, and PPO used the same function approximation architecture with two hidden layers, 64 units in each layer, and the tanh activation. TRPO, PPO, DDPG, and TD3 used their default hyper-parameter settings. We used Xavier initialization (Glorot & Bengio, 2010) for TDPO, and multiplied the outputs of the network by a factor of 0.001 so that the initial actions were small and nearly zero. We also confirmed that TDPO does not decrease in performance when using an identical network initialization to PPO/TRPO. TD3’s baseline implementation was amended to support MPI parallelization just like TRPO, PPO, and DDPG. To produce the results for DDPG and TD3, we used hyperparameter optimization both with and without the tanh final activation function that is common for DDPG and TD3 (this causes the difference in initial payoff in the figures). However, under no conditions were DDPG and TD3 able to solve these environments effectively, suggesting that the deterministic search used by TDPO is operating in a qualitatively different way than the stochastic policy optimization used by DDPG and TD3. Note that we made a thorough attempt to compare DDPG and TD3 fairly, including trying different initializations, different final layer scalings/activations, different network architectures, and performing hyperparameter optimization. Mini-batch selection was unnecessary for TDPO since optimization for samples generated by DeVine was fully tractable. The confidence intervals in all figures were generated using 1000 samples of the statistics of interest.
For designing the environment, we used Dhariwal et al. (2017)’s pendulum dynamics and relaxed the torque thresholds to be as large as 40 N m. The environment also had the same episode length of 200 time-steps. We used the reward function described by the following equations:
R(st, at) = CR ·R(τ) · 1{t = 200} R(τ) = RFreq(τ) +ROffset(τ) +RAmp(τ)
RFreq(τ) = 0.1 · [ fmax∑ f=fmin Θ+std(f) 2 − 1 ]
ROffset(τ) = − ∣∣∣∣Θ(f = 0)200 − θTarget Offset ∣∣∣∣ = −∣∣∣∣( 1200 199∑ t=0 θt ) − θTarget Offset ∣∣∣∣ RAmp(τ) = hpiecewise ( ΘAC θTarget Amp. − 1 )
(67)
where
• θ is the pendulum angle signal in the time domain. • Θ is the magnitude of the Fourier transform of θ. • Θ+ is the same as Θ only for the positive frequency components. • ΘAC is the normalized oscillatory spectrum of Θ:
ΘAC =
√ Θ+TΘ+
200 . (68)
• hpiecewise is a piece-wise linear error penalization function:
hpiecewise(x) = −x · 1{x ≥ 0}+ 10−4x · 1{−x ≥ 0}. (69)
• Θ+std is the standardized positive amplitudes vector:
Θ+std = Θ+√
Θ+TΘ+ + 10−6 . (70)
• CR = 1.3×104 is a reward normalization coefficient, and was chosen to yield approximately the same payoff as a null policy would yield in the typical pendulum environment of Dhariwal et al. (2017).
• θTarget Offset is the target offset, θTarget Amp. is the target amplitude, and [fmin, fmax] is the target frequency range of the environment.
All methods used 48 parallel workers. The machines used Xeon E5-2690-v3 processors and 256 GB of memory. Each experiment was repeated 25 times for each method, and each run was given 6 hours or 500 million samples to finish.
A.10 IMPLEMENTATION DETAILS FOR THE ENVIRONMENT WITH LONG HORIZON AND RESONANT FREQUENCIES
For the robotic leg, we used exactly the same algorithms with the same parameters as described in Section A.9 above.
We used the reward function described by the following equations:
R = Rposture +Rvelocity +Rfoot offset +Rfoot height +Rground force +Rknee height +Ron-air torques (71)
with
Rposture = −1× [∣∣∣∣θknee + π2 ∣∣∣∣+ ∣∣∣∣θhip + π4 ∣∣∣∣]
Rvelocity = −0.08× [ |ωknee|+ |ωhip| ] Rfoot offset = −10× [ |xfoot| · 1{zknee < 0.2}
] Rground force = −1× [ |fz −mg| · 1{fz < mg} · 1touchdown
] Rfoot height = −1× [ |zfoot| · 1touchdown
] Rknee height = −15×
[∣∣zknee − ztargetknee ∣∣ · 1touchdown] Ron-air torques = −10−4 × [ (τ2knee + τ 2 hip) · (1− 1touchdown) ] (72)
where
• θknee and θhip are the knee and hip angles in radians, respectively.
• ωknee and ωhip are the knee and hip angular velocities in radians per second, respectively.
• xfoot and zfoot are the horizontal and vertical foot offsets in meters from the desired standing point on the ground, respectively.
• xknee and zknee are the horizontal and vertical knee offsets in meters from the desired standing point on the ground, respectively.
• fz is the vertical ground reaction force on the robot in Newtons.
• m is the robot weight in kilograms (i.e., m = 0.76 kg).
• g is the gravitational acceleration in meters per second squared.
• 1touchdown is the indicator function of whether the robot has ever touched the ground.
• ztargetknee is a target knee height of 0.1 m.
• τknee and τhip are the knee and hip torques in Newton meters, respectively.
All methods used 72 full trajectories between each policy update, and each run was given 16 hours of wall time, which corresponded to almost 500 million samples. This experiment was repeated 75 times for each method. The empirical mean of the discounted payoff values were reported without any performance or seed filtration. The same hardware as the non-local rewards experiments (i.e., Xeon E5-2690-v3 processors and 256 GB of memory) was used.

heading: A.11 GYM SUITE BENCHMARKS
text: While it is clear that our deterministic policy gradient performs well on the new control environments we consider, one may naturally wonder about its performance on existing RL control benchmarks. We ran our method on a suite of Gym environments and include four representative examples in Figure 3. Broadly speaking, our method (TDPO) performs slightly worse on average than others, but occasionally performs much better as seen in the Swimmer-v3 environment. We speculate that these gym environments are reasonably robust to injected noise, and this may mean that stochastic policy gradients can more rapidly and efficiently explore the policy space, or there may be other algorithmic enhancements that are needed for fully deterministic policy gradients in these cases. The experiments granted each method 72 parallel MPI workers for about 144 million steps (i.e., 2 million sequential steps), and the returns were averaged over 100 different seeds for each method. Since the computational cost of running both DDPG and TD3 were high, we only included TD3 since it was shown to outperform DDPG in earlier benchmarks.

heading: A.12 RUNNING TIME COMPARISON
text: Figure 4 depicts a comparison of each method’s running time per million steps. These plots show the combination of both the simulation (i.e., environment sampling) and the optimization (i.e., computing the policy gradient and running the conjugate gradient solver) time. It is clear that our method (TDPO) is generally faster than the other algorithms. This is mainly due to the computational efficiency of the DeVine gradient estimator, which summarizes two full trajectories in a single state-action-advantage tuple which can significantly reduce the optimization time. That being said, these relative comparisons could vary to a large extent (1) under different processor architectures, (2) with more (or less) efficient implementations, or (3) when running environments whose simulation time constitutes a significantly larger (or smaller) portion of the total running time.

heading: A.13 OTHER SWINGING PENDULUM VARIANTS
text: Multiple variants of the pendulum with non-local rewards were used, each with different frequency targets and the same reward structure. Table 3 summarizes the target characteristics of each variant.
The main variant was shown in the paper. Figures 5, 6, 7, 8, 9, 10, 11, and 12 show similar results for the second to ninth variants. To focus on our method’s ability to solve all these variants efficiently, we only show the performance of our method (TDPO) on all variants in Figure 13. Overall, we found TRPO, PPO, DDPG, and TD3 to occasionally find the correct offset. They either excited the natural or the maximum (not the desired) frequency of the pendulum, but they were not able to drive the desired frequency and amplitude. TDPO was able to achieve the desired oscillations (and thus high rewards) in all variants.

heading: A.14 NOTES ON HOW TO IMPLEMENT TDPO
text: In short, our method (TDPO) is structured in the same way TRPO was structured; both TDPO and TRPO use policy gradient estimation, and a conjugate-gradient solver utilizing a Hessian-vector product machinery. On the other hand, there are some algorithmic differences that distinguish TDPO from TRPO. First of all, TRPO uses line-search heuristics to adaptively find the update scale; no such heuristics are applied in TDPO. Second, TDPO uses the DeVine advantage estimator, which requires storing and reloading pseudo-random generator states. Finally, the Hessian-vector product machinery used in TDPO computes Wasserstein-vector products, which is slightly different from those used in TRPO. The hyper-parameter settings and notes on how to choose them were discussed
in Sections A.9, A.6, and A.10. We will describe how to implement TDPO, and focus on the subtle differences between TDPO and TRPO next.
As for the state-reset capability, our algorithm does not require access to a reset function to arbitrary states. Instead, we only require to be able to start from the prior trajectory’s initial state. Many environments, including the Gym environments, instantiate their own pseudo-random generators and only utilize that pseudo-random generator for all randomized operations. This facilitates a straightforward implementation of the DeVine oracle; in such environments, implementing an arbitrary state-reset functionality is unnecessary, and only reloading the pseudo-random generator to its configuration prior to the trajectory would suffice. In other words, the DeVine oracle can store the initial configuration of the pseudo-random generator before asking for a trajectory reset, and then start sampling. Once the main trajectory is finished, the pseudo-random generator can be reloaded, thus producing the same initial state upon a reset request. Other time-step states can then be recovered by applying the same preceeding action sequence.
To optimize the quadratic surrogate, the conjugate gradient solver was used. Implementing the conjugate gradient algorithm is fairly straightforward, and is already included in many common automatic differentiation libraries. The conjugate gradient solver is perfect for situations where (1) the Hessian matrix is larger than can efficiently be stored in the memory, and (2) the Hessian matrix includes many nearly identical eigenvalues. Both of these conditions apply well for TDPO, as well as for TRPO. Instead of requiring the full Hessian matrix to be stored, the conjugate gradient solver only requires a Hessian-vector product machinery v → Hv, which must be specifically implemented for TDPO. Our surrogate function can be viewed as
L(δθ) = gT δθ + C ′ 2
2 δθTHδθ
where the Hessian matrix can be defined as
H = H2 + C ′1 C ′2 H1,
H1 := ∇2θ′Es∼ρπkµ [ LG2(π′, πk; s) ] ,
H2 := ∇2θ′Es∼ρπkµ
[ W (π′(a|s), πk(a|s))2 ] . (73)
In order to construct a Hessian-vector product machinery v → Hv, one can design an automaticdifferentiation procedure that returns the Hessian-vector product. Many automatic-differentiation packages already include functionalities that can provide a Hessian-vector product machinery of a given scalar loss function without computing the Hessian matrix. This can be used to implement the Hessian-vector product machinery in a straightforward manner; one only needs to provide the scalar quadratic terms of our surrogate, and would obtain the Hessian-vector product machinery in return. On the other hand, this may not be the most computationally efficient approach, as our problem exhibits a more specific structure. Alternatively, one can implement a more elaborate and specifically designed Hessian-vector product machinery by following these three steps:
• Compute the Wasserstein-vector product v → H2v according to Algorithm 3. • Compute the Sensitivity-vector product v → H1v according to Algorithm 4. • Return the weighted sum of H1v and H2v as the final Hessian-vector product Hv.
One may also need to add a conjugate gradient damping to the conjugate gradient solver (i.e., return βv + Hv for some small β as opposed to returning Hv purely), which is also done in the TRPO method. This may be important when the number of policy parameters is much larger than the sample size. Setting β = 0 may yield poor numerical stability if H had small eigenvalues, and setting large β will cause the conjugate gradient optimizer to mimic the gradient descent optimizer by making updates in the same direction as the gradient. The optimal conjugate gradient damping may depend on the problem and other hyper-parameters such as the sample size. However, it can easily be picked to be a small value that ensures numerical stability.
Once the conjugate gradient solver returned the optimal update direction H−1g, it must be scaled down by a factor of C ′2 (i.e., δθ ∗ = H−1g/C ′2). If δθ ∗ satisfied the trust region criterion (i.e.,
1 2δθ ∗THδθ∗ ≤ δ2max), then one can make the parameter update (i.e., θnew = θold + δθ∗) and proceed to the next iteration. Otherwise, the proposed update δθ∗ must be scaled down further, namely by α, such that the trust region condition would be satisfied (i.e., 12 (αδθ
∗)TH(αδθ∗) = δ2max) before making the update θnew = θold + αδθ∗.
Algorithm 3 Wasserstein-Vector-Product Machinery Require: Current Policy π1 with parameters θ1. Require: The vector v with the same dimensions as θ1. Require: An observation s.
1: Compute the action for the observation s with |A| elements.
a|A|×1 :=  π (1)(s)
... π(|A|)(s)  . (74) This vector should be capable of propagating gradients back to the policy parameters when used in automatic differentiation software. 2: Define t to have be a constant vector with the same shape as a. It could be populated with any values such as all ones. 3: Define the scalar ã := aT t. 4: Using back-propagation, find the gradient
∇θã = |A|∑ i=1 ti∇θai = |A|∑ i=1 ti [ ∂ai ∂θ1 · · · ∂ai∂θ|Θ| ] . (75)
5: Compute the following dot-product:
〈∇θã, v〉 = ( |A|∑ i=1 ti · ∂ai ∂θ1 ) · v1 + · · ·+ ( |A|∑ i=1 ti · ∂ai ∂θ|Θ| ) · v|Θ|. (76)
6: Using automatic differentiation, take the gradient w.r.t. the t vector.
ãθ,v := ∇t〈∇θã, v〉 =  ∂a1 ∂θ1 · v1 + · · ·+ ∂a1∂θ|Θ| · v|Θ|
... ∂a|A| ∂θ1 · v1 + · · ·+ ∂a|A| ∂θ|Θ| · v|Θ|
 =  ∂a1 ∂θ1 · · · ∂a1∂θ|Θ| ...
... ∂a|A| ∂θ1 · · · ∂a |A| ∂θ|Θ|  v (77) 7: Compute the dot product 〈ãθ,v, ã〉. 8: Using back-propagation, take the gradient w.r.t. θ, and return it as the gain-vector-product.
∇θ〈ãθ,v, ã〉 =  ∂a1 ∂θ1 · · · ∂a1∂θ|Θ| ...
... ∂a|A| ∂θ1 · · · ∂a |A| ∂θ|Θ|
 T  ∂a1 ∂θ1 · · · ∂a1∂θ|Θ| ...
... ∂a|A| ∂θ1 · · · ∂a |A| ∂θ|Θ|
 v (78)
Algorithm 4 Sensitivity-Vector-Product Machinery Require: Current Policy π1 with parameters θ1. Require: The vector v with the same dimensions as θ1. Require: An observation s.
1: Compute the action to observation Jacobian matrix
J|A|×|S| :=  ∂π(1)(s) ∂s1 · · · ∂π(1)(s) ∂s(|S|) ... ...
∂π(|A|)(s) ∂s(1) · · · ∂π (|A|)(s) ∂s|S|  . (79) This can either be done using finite-differences in the observation using
∂π(i)(s)
∂s(j) ' π (i)(s+ ds · ej)− π(i)(s) ds
(80)
(which may be a bit numerically inaccurate), or using automatic differentiation. In any case, this matrix should be a parameter tensor capable of propagating gradients back to the parameters when used in automatic differentiation software. 2: Define J̃ to be the vectorized (i.e., reshaped into a column) J matrix, with |AS| = |A| × |S| rows and one column. 3: Define t to have be a constant vector with the same shape as J̃ . It could be populated with any values such as all ones. 4: Define the scalar Jt := J̃T t. 5: Using back-propagation, find the gradient
∇θJt = |A|∑ i=1 |S|∑ j=1 ti,j∇θJi,j = |A|∑ i=1 |S|∑ j=1 ti,j [ ∂Ji,j ∂θ1 · · · ∂Ji,j∂θ|Θ| ] . (81)
6: Compute the following dot-product.
〈∇θJt, v〉 = ( |A|∑ i=1 |S|∑ j=1 ti,j · ∂Ji,j ∂θ1 )× v1 + · · ·+ ( |A|∑ i=1 |S|∑ j=1 ti,j · ∂Ji,j ∂θ|Θ| )× v|Θ| (82)
7: Using automatic differentiation, take the gradient w.r.t. the t vector.
(∇θJ)v := ∇t〈∇θJt, v〉 =  ∂J1,1 ∂θ1 · v1 + · · ·+ ∂J1,1∂θ|Θ| · v|Θ|
... ∂J|A|,|S| ∂θ1 · v1 + · · ·+ ∂J|A|,|S| ∂θ|Θ| · v|Θ|

=  ∂J(1,1) ∂θ1 · · · ∂J (1,1) ∂θ|Θ| ...
... ∂J(|A|,|S|)
∂θ1 · · · ∂J
(|A|,|S|)
∂θ|Θ|  v (83) 8: Reshape (∇θJ)v into a column vector and name it J̃θ,v . 9: Compute the dot product 〈J̃θ,v, J̃〉.
10: Using back-propagation, take the gradient w.r.t. θ, and return it as the gain-vector-product.
∇θ〈J̃θ,v, J̃〉 =  ∂J(1,1) ∂θ1 · · · ∂J (1,1) ∂θ|Θ| ...
... ∂J(|A|,|S|)
∂θ1 · · · ∂J
(|A|,|S|)
∂θ|Θ|
 T  ∂J(1,1) ∂θ1 · · · ∂J (1,1) ∂θ|Θ| ...
... ∂J(|A|,|S|)
∂θ1 · · · ∂J
(|A|,|S|)
∂θ|Θ|
 v (84)

Review summaries:
The paper proposes a deterministic policy gradient method using the Wasserstein distance to quantify the difference of deterministic policies. It also introduces a new deterministic estimator for the policy gradient with theoretical justification. The paper uses a similar path as TRPO to first introduce a surrogate objective then replace the complex term by a tractable term for practical performance. The method has been shown to work well under long-horizon and non-local reward environments whereas existing methods struggle to solve.

This paper introduces a policy gradient method based on deterministic policies and deterministic gradient estimates. By assuming such a deterministic setting, the authors show that the proposed technique can estimate gradients on long-horizon tasks without the need to inject noise into the system for exploration. The authors propose using a Wasserstein-based model to perform regularization in this particular setting involving deterministic policies. They also establish the conditions under which it is possible to guarantee monotonic policy improvement.

This paper proposes a new deterministic policy gradient method (TDPO). The main idea and its derivations are based on the use of a deterministic Vine (DeVine) gradient estimator and the Wasserstein metric. The paper shows that a closed-form computation of Wasserstein distance can be derived without any approximation on the deterministic policies. A monotonic policy improvement guarantee is provided. Deterministic Vine (DeVine) is proposed as an optimization method for the proposed surrogate loss function, which is also inspired by the TRPO's surrogate loss. The proposed method is evaluated on two new robotic control domains that are proposed by the authors. The comparisons are against PPO, TRPO, DDPG, and TD3. Additional experiment results on Open AI Gym environments are provided in Appendix.