DPO: Differential reinforcement learning with application to optimal configuration search
Chandrajit Bajaj Minh Nguyen
















Minh, Bajaj


Department of Computer Science 

The University of Texas at Austin

Department of Mathematics

The University of Texas at Austin

Oden Institute for Computational Engineering and Sciences

The University of Texas at Austin


DPO: Differential reinforcement learning with application to optimal configuration search
    Chandrajit Bajaj 1, 3Minh Nguyen 2, 3
    
=========================================================================================






Reinforcement learning (RL) with continuous state and action spaces remains one of the most challenging problems within the field. Most current learning methods focus on integral identities such as value functions to derive an optimal strategy for the learning agent. In this paper, we instead study the dual form of the original RL formulation to propose the first differential RL framework that can handle settings with limited training samples and short-length episodes. Our approach introduces Differential Policy Optimization (DPO), a pointwise and stage-wise iteration method that optimizes policies encoded by local-movement operators. We prove a pointwise convergence estimate for DPO and provide a regret bound comparable with current theoretical works. Such pointwise estimate ensures that the learned policy matches the optimal path uniformly across different steps. We then apply DPO to a class of practical RL problems which search for optimal configurations with Lagrangian rewards. DPO is easy to implement, scalable, and shows competitive results on benchmarking experiments against several popular RL methods.



§ INTRODUCTION

Reinforcement learning (RL) is a powerful machine learning paradigm with applications in many areas of engineering and biological sciences <cit.>. In RL, an agent interacts with the outside environment to learn the best action maximizing a certain cumulative reward. Such unsupervised and feedback-based learning approach resembles human thinking, giving promising solutions to challenging tasks without proper labeled data. Nevertheless, RL is not widely applied and popularized as a large-scale industrial solution like other fields such as computer vision or natural language processing. In fact, current RL algorithms are more complex and sample-expensive than their deep learning counterparts. Thus, they offer lower-quality solutions under limited data settings and restricted theoretical sample-efficiency guarantees except in special cases. 


In this paper, we aim to make a very first step towards practical RL algorithms, particularly in the case of continuous state and action spaces, for limited data scenarios. Mathematically speaking, an RL agent seeks an optimal policy, the function mapping from the agent's current state to its optimal action. RL algorithms rely on the cumulative reward objective to assess a particular policy. In applications involving physical systems, the path of the learned actions and states ideally should resemble the optimal path across every time steps to capture the underlying physics. Unfortunately, under the extreme situations of limited training samples, this is met with several challenges. First, since the value function is a sum quantity, the learned policy can put more bias on a particular term and less focus on others to temporarily boost the sum value. Such policy might also search for the best state/action in terms of rewards without concerning about the optimal path as a whole, leading to a sub-optimal path that don't conform to physics laws. Second, convergence guarantees <cit.> for RL algorithms mainly look at the value function's improvement, making it difficult to assess the learned path's quality across all individual time steps. These weaknesses motivate us to consider an alternative formulation with exclusive focus on the dynamics induced by the optimal policy and individual data points along the learned trajectories (see <ref>). As a direct solution to this new formulation, referred to as differential RL, we develop differential policy optimization (DPO), a simple yet effective policy optimization algorithm. We then prove a theoretical pointwise convergence estimate for DPO that allows policy's assessment on the whole path. This estimate is also crucial to our derivation of the theoretical sample-efficiency guarantees in the next step. 


To theoretically evaluate an RL algorithm's sample efficiency, researchers often use a quantity called regret, which is the gap between the cumulative rewards of the trained policy and the optimal policy over all training samples and episodes. While several research works <cit.> have addressed the discrete settings with their optimal regret bounds, the case of continuous state and action spaces remains challenging. Under a generic assumption, RL regret bounds can be untractable. Thus, a mild assumption of Lipschitz MDP <cit.> was applied, and the corresponding regret has a lower bound of Ω(K^d+1/d+2) <cit.>, where K is the number of training episodes and d is the dimension of the combined state-action space. Within our framework, we also prove a regret bound for DPO of similar order O(K^2d+3/2d+4) (see <ref>). Nonetheless, these bounds come with some performance deterioration due to their dependence on state-action dimension. Several works impose stronger assumptions to obtain more desirable bounds. Maran et al. <cit.>, for instance, considers cases with smooth transition dynamics and rewards up to order ν with a regret bound O(K^3d/2+ν+1/d+2(ν+1)), which becomes the desirable O(K^1/2) when ν→∞. Vakili and Olkhovskaya<cit.> study transition and reward functions in reproducing kernel Hilbert space (RKHS) induced by the Matern covariance function with parameter m. This assumption yields the the regret bound O(K^d+m+1/d+2m+1) that again becomes O(K^1/2) when m →∞. Our work takes a different approach to achieve more helpful regret bounds. Instead of considering more specialized problems, we propose a narrower hypothesis class for policy approximation. In particular, we assume that black box environment can provide a hypothesis class where functions differ from the true policy either a d-weakly convex or concave function and linearly bounded (see <ref> and <ref>). DPO algorithm then can attain a regret bound of order O(K^5/6) (see <ref>) that is independent of the state-action dimension.


Contributions. Our main contributions can be summarized as:


  * We propose a novel learning framework with a practical policy optimization algorithm called differential policy optimization (DPO) (see <ref>) to address RL problem in differential form, where the dynamics is evaluated directly without expected cumulative rewards. 

  * We provide a rigorous proof for DPO algorithm's pointwise convergence estimates together with its regret bounds analysis.

  * We demonstrate DPO's effectiveness in limited data settings through application to complex physics-based and engineering problems with Lagrangian energy rewards.


Organization. In Section 2, we introduce a novel framework of differential reinforcement learning and its learning algorithm called DPO. Here, we also demonstrate the duality between our framework and original RL formulation. In Section 3, we outline lemmas to prepare for the convergence proof and regret analysis for the DPO algorithm in Section 4. In Section 5, we apply DPO algorithm to three RL tasks with Lagrangian energy cost, and perform benchmarkings against TRPO<cit.>, PPO <cit.> and SAC <cit.>.



§ DIFFERENTIAL REINFORCEMENT LEARNING
 


 §.§ Problem statement

Differential reinforcement learning's main goal is to solve the problem 𝒟 which aims to find the optimal dynamics G: Ω→Ω underlying the optimal trajectory:

    x_0 = X ∼ρ_0, x_1 = G(X), x_2 = G(x_1) = G^(2)(X), ⋯, x_H-1 = G^(H-1)(X)

Here Ω is a compact domain in ℝ^d, and H is the number of steps in an episode. G^(k)(X) is the k-times composition function G(G(⋯ G(x) ⋯)). Moreover, ρ_0 is the distribution of the starting point x_0. In this work, we assume ρ_0 is a continuous function. By interacting with an environment ℬ (see <ref>), we learn a policy G_θ parameterized by θ that approximates G.


An environment ℬ with respect to an adversarial distribution ρ_0 and a score function g is a black-box system that allows you to check the quality of a policy G_θ. More concretely, ℬ inputs a policy G_θ and outputs the trajectories (G_θ^(k)(x))_k=0^H-1, and their associated scores (g(G_θ^(k)(x)))_k=0^H-1. Recall that H is the number of steps in an episode.


The function g allows one to assess the performance of the current policy G_θ and then perform necessary updates to get to the the desired dynamics G. In this paper, we assume that there is a first-order relation between the score function g and the dynamics function G:

    G = A · dg,  where  A  is a known constant matrix

This relation, which generalizes the assumption under a Hamiltonian system, will be explained in <ref>.



 §.§ Differential policy optimization (DPO) algorithm

Differential policy optimization (DPO) (see <ref>) is an “time-expansion (Dijkstra-like)” algorithm that iteratively uses approriate on-policy samples for policy update to ensure an increase in policy quality over each iteration. This algorithm has a similar idea as Trust Region Proximal Policy Optimization. However, thanks to the differential approach that focuses on pointwise estimates, DPO becomes much simpler and easier to implement/optimize in practice compared to its RL counterpart.





 §.§ Original RL formulation

We first recall the usual RL mathematical formulation, which has a general idea similar to <ref> but with different details. A finite-horizon discounted Markov decision process (MDP) is defined by the 6-tuple (𝒮, 𝒜, , r, ρ_0, γ), where 𝒮 and 𝒜 be the set of states and of actions, : 𝒮×𝒜×𝒮→ℝ is the transition probability distribution, r: 𝒮×𝒜→ℝ is the reward function, ρ_0: 𝒮→ℝ is the distribution of initial state and γ∈ (0, 1) is the discount factor. At time t, the agent choose an action a_t ∈𝒜 given its current state s_t ∈𝒮 based on π(a_t|s_t):

    s_0 ∼ρ_0(s_0), a_t ∼π(a_t|s_t), s_t+1∼(s_t+1|s_t, a_t)

The goal is to maximize the expected γ-discounted cumulative reward 𝒥 = _π[∑_t=0^H-1γ^t r(s_t, a_t)]. For a particular policy π, let V_π denote the following value function obtained by executing π:

    V_π(s_t) = _a, s_1, ⋯[∑_t=0^H-1γ^t r(s_t, a_t) | s_0 = s]

The value function is then defined as V(s): = _π V_π(s) and is used across several RL algorithms. 



 §.§ Duality between Standard RL and Differential RL

To illustrate the duality between two approaches, we first consider a more specialized problem with Lagrangian reward (kinetic minus potential):

    π^* = max_π[∑_k=0^H-1γ^-k(a_t_k^2/2 - ℱ(s_t_k)) ] 
       subject to the linear dynamics: s_k+1 = s_k + Δ a_k,  for step size Δ.

Here ℱ is a potential energy function unknown to the learning agent on the state space 𝒮. The action a_t_k = π(s_t_k)  be the output action from a policy π at discrete time t_k = kΔ. T = t_H = HΔ is the terminal time, and γ is the magnification factor on the energy. To see the duality, we first look at the continuous-version of <ref>, which is the optimal control problem of finding the path that optimizes:

    J    = ∫_0^T β^-t(a^2/2 - ℱ(q)) dt 
       subject to ṡ = a

where γ^-k = β^-t_k = β^-Δ k so that we have the relation γ = β^Δ. With this continuous-time format, the duality can be expressed via the following <ref>


The optimal trajectory/dynamics of <ref> in the generalized coordinates has the form:

    ṡ   = β^t p 
    ṗ    = -β^-t∇ℱ(s)



First of all, the associated Hamiltonian HF(p, s, a) for <ref> is pa - β^-t(a^2/2 - ℱ(q)). By Pontryagin maximum principle <cit.>, HF_a = p - β^-ta = 0 so that the optimal action a^* = β^tp. Hence, the reduced Hamiltonian hf(s, p):= HF(s, p, a^*(s, p)) has the form:

    hf(s, p) = β^tp^2 - β^t p^2/2 + β^-tℱ(s) = β^t p^2/2 + β^-tℱ(s)

As a result, other conditions from Pontryagin maximum principle imply that:

    ṡ   = hf_p = β^t p 
    ṗ   = -hf_s = -β^-t∇ℱ(s)



The above dynamics <ref> is called optimal control flow (OCF). Agents learning this type of dynamics strive to arrive at a stable position with low potential energy ℱ. By considering x = (s, p), our RL problem can be approximated by solving the abstract problem 𝒟 in <ref> with the score function g(x) = g(s, p) = 0.5 p^2 + ℱ(s) and a constant sympletic matrix A. Similar to this specific example, duality between the traditional RL and our differential formulation can be obtained via the reduced Hamiltonian function and is explicitly expressed through the Legendre transform <cit.>. This type of duality is a generalized version of the duality between Lagrangian mechanics and Hamiltonian mechanics in classical physics. The latter, which is similar to our formulation, imposes a physical condition regarding the Hamiltonian's derivative that leads to a more regularized path:

    ṡ = hf_p  and ṗ = -hf_s

In the differential form <ref>, we couple state s and momentum p, which contains information from both state and action, into a single vector x = (s, p) with dimension d = d_𝒮 + d_𝒜 (sum of state and action spaces' dimensions). The above Hamiltonian condition is a special case of <ref>.


We also note an interesting property of the optimal solution for <ref>. From  <ref>, for a = β^t p, ȧ = β^tṗ + (logβ) β^t p = -∇_s ℱ(s) + logβ a. Hence, <ref> is equivalent to:

    ṡ   = a 
    ȧ   = (-logβ) a + ∇_s ℱ_s

or

    [ ṡ; ȧ ] + [         0; (-logβ) a ] - [  0  I; -I  0 ][ ∇_s hf; ∇_a hf ] = 0

This is similar to the continuous-version of the Nesterov accelerated method, which has the form of a damped Hamiltonian flow: 

    [ ẋ_̇ṫ; ṗ_̇ṫ ]+[       0; α_t p_t ] - [  0  I; -I  0 ][ ∇_x HF(x_t,p_t); ∇_p HF(x_t,p_t) ] = 0

where HF is the associated Hamiltonian flow. 



§ BASIC NOTATIONS AND SUPPORTING LEMMAS

In this section, we state all supporting lemmas to prove the main estimate for <ref> (DPO) in <ref>.




<ref> defines the number of training samples needed to allow derivative approximation transfer. This definition is used to derive the number of samples needed for <ref> (DPO).

Recall that ρ_0 is the a continuous density for starting states. Suppose we are given a function g: Ω→Ω, a hypothesis space ℋ consists of function h ∈ℋ that approximates g, two positive constants ϵ and δ. We define the function N(g, ℋ, ϵ, δ) to be the number of samples needed so that if we approximate g by h ∈ℋ via N(g, ℋ, ϵ, δ) training samples, then with probability of at least 1 - δ, we have the following estimate bound on two function gradients:

    ∇ g(X) - ∇ h(X) < ϵ

In other words, we want N(g, ℋ, ϵ, δ) to be large enough so that the original approximation can transfer to the derivative approximation above. If no such N(g, ℋ, ϵ, δ) exists, let N(g, ℋ, ϵ, δ) = ∞


<ref> and <ref> below estimates N(g, ℋ, ϵ, δ) (defined in <ref>) in terms of the required threshold error ϵ. 
[]lemmaDerivativeApproxGeneralHypothesis
Suppose that the hypothesis space ℋ for approximating the target function g: Ω⊂ℝ^d →ℝ^d consists of neural network appproximators with bounded weights and biases. Further assume that function g and neural network appproximators h ∈ℋ are twice continuously differentiable with bounded first and second derivative by some constant C. One way this assumption can be satisfied is to choose activation functions that are twice continuously differentiable. Then N(g, ℋ, ϵ, δ) is upper-bounded by O(ϵ^-(2d+4)), where we only ignore quadratic factors of δ, polynomial terms of d, and other logarithmic terms.



For p ∈ℕ, a function h is called a p-weakly convex function if for any x ∈Ω, there exists a sufficiently small neighborhood U of x so that:

    h(y) ≥ h(x) + ∇ h(x)(y-x) - C y-x^p

Note that any convex function is p-weakly convex for any p ∈ℕ.



A function h is linearly bounded, if for some constant C > 0, and for any x ∈Ω, there exists a sufficiently small neighborhood U of x so that:

    |h(y) - h(x)| ≤ C∇ h(x)y-x
 

Any Lipschitz function on compact domain that has non-zero derivative is linearly bounded. One such example is the function e^αx^2 on domains that don't contain 0.

[]lemmaDerivativeApproxWeakConvexHypothesis
Suppose that, possibly with knowledge from outside environment or from certain policy experts, the hypothesis space ℋ is reduced to a smaller family of functions of the form g + h, where h is a linearly bounded and p-weakly convex function for some p ≥ d. Then N(f, ℋ, ϵ, δ) can be upper-bounded by the factor O(ϵ^-6) that is independent of the dimension d. 



There are many hypothesis spaces consisting of infinitely many elements that satisfy <ref>. One such hypothesis space is: 

    ℋ = h(x):=∑_i=1^D a_i e^b_i x^2,   0 ≤ a_i ≤ A, 0 < b ≤ b_i ≤ B

for constant A, B, b > 0 and for domain Ω doesn't contain 0. <ref> also holds for concave functions and their weak versions.


Finally, we state the supporting <ref> for the Main <ref> which proves the pointwise estimates for DPO algorithm.
[]lemmaBoundSeq
Given L and ϵ > 0, define two sequences α_j_j ≥ 0 and ϵ_j_j ≥ 0  recursively as follows:

    α_0 = 0     and α_j = Lα_j-1 + ϵ
    ϵ_1 = ϵ    and ϵ_k+1 = Lα_k + ϵ + Lϵ_k

Then for each k, we get:

    ϵ_k ≤kL^kϵ/L-1




§ THEORETICAL ANALYSIS FOR DIFFERENTIAL POLICY OPTIMIZATION

This section shows the convergence of differential policy optimization (DPO, <ref>) based on generalization pointwise estimates. We then use this result to derive suitable regret bounds for DPO.



 §.§ Main theorem


Suppose we're given a threshold error ϵ, a probability threshold δ, and number of steps per episode H. Assume that N_k_k = 1^H-1 is the sequence of numbers of samples used at each stage <ref> (DPO) so that:

    N_1    = N(g, ℋ_1, ϵ, δ),  and for  k ∈2, H-1:
    
    N_k    = maxN(g_θ_k-1, ℋ_k, ϵ, δ_k-1/(k-1)), N(g, ℋ_k, ϵ, δ_k-1/(k-1))

We further assume that there exists a Lipschitz constant L > 0 such that both the true dynamics G and the policy neural network approximator G_θ_k at step k with regularized parameters have their Lipschitz constant at most L for each k ∈1, H. Then, for a general starting point X, with probability at least 1 - δ, the following generalization bound for the trained policy G_θ_k holds for all k ∈1, 2, ⋯, H-1:

    _XG_θ_k^(j)(X) - G^(j)(X) < jL^jϵ/L-1 for all 
      1 ≤ j ≤ k

Note that when N_k →∞, the errors approach 0 uniformly for all j given a finite terminal time T.



Define δ_k = δ/3^m-k = 3 δ_k-1. Let α_k and ϵ_k be two sequences associated with Lipschitz constant L and threshold error ϵ as in <ref>. We prove the generalization bound statement by induction on the stage number k that for probability of at least 1 - δ_k,

    _XG_θ_k^(j)(X) - G^(j)(X) < ϵ_j  for all 
      1 ≤ j ≤ k

By <ref>, proving this statement also proves <ref>.


The bound for the base case k = 1 is due to the definition of N_1 = N(g, ℋ_1, ϵ, δ) that allows the approximation of g by g_θ_1 transfers to their derivatives G and G_θ_1 with error threshold ϵ and probability threshold δ. Assume that the induction hypothesis is true for k. We prove that for a starting (random variable) point X, the following error holds with a probability of at least 1 - δ_k+1 = 1 - 3δ_k:

    _XG_θ_k+1^(j)(X) - G^(j)(X) < ϵ_j  for all 
      j ≤ k+1


First, from induction hypothesis, with probability of at least 1 - δ_k:

    _XG_θ_k^(j)(X) - G^(j)(X) < ϵ_j  for all 
      j ≤ k


In stage k+1, all previous stages' samples up to stage k-1 is used for G_θ_k+1. As a result, we can invoke the induction hypothesis on k to yield the same error estimate for G_θ_k+1 on the first k sample points with probability 1 -δ_k:

    _XG_θ_k+1^(j)(X) - G^(j)(X) < ϵ_j  for all 
      j ≤ k


Recall from <ref> that g_θ_k+1 is trained to approximate g_θ_k to ensure that the updated policy doesn't deviate too much from current policy. For j ∈1, ⋯, k-1, g_θ_k+1∈ℋ_k+1 approximates g_θ_k on N_k+1 samples of the form G_θ_k^(j)(X^i) for i ∈1, ⋯, N_k+1. Since N_k+1≥ N(g_θ_k, ℋ_k+1, ϵ, δ_k/k) allows derivative approximation transfer, for probability of at least 1 - δ_k/k, _XG_θ_k+1(G_θ_k^(j)(X)) - G_θ_k(G_θ_k^(j)(X)) < ϵ. Hence, under a probability subspace Γ with probability of at least (1 - (k-1)δ_k/k), we have:

    _XG_θ_k+1(G_θ_k^(j)(X)) - G_θ_k(G_θ_k^(j)(X)) < ϵ for all  1 ≤ j < k


We prove by induction on j that under this probability subspace Γ, we have:

    _XG_θ_k+1^(j)(X) - G_θ_k^(j)(X) < α_j  for all 
      1 ≤ j ≤ k


In fact, for the induction step, one get:

    _XG_θ_k+1^(j)(X) - G_θ_k^(j)(X)≤_XG_θ_k+1(G_θ_k+1^(j-1)(X)) - G_θ_k+1(G_θ_k^(j-1)(X))
          + _XG_θ_k+1(G_θ_k^(j-1)(X)) - G_θ_k(G_θ_k^(j-1)(X))
       ≤ L_XG_θ_k+1^(j-1)(X) - G_θ_k^(j-1)(X) + ϵ≤ Lα_j-1 + ϵ = α_j


Finally, we look at the approximation of g by g_θ_k+1∈ℋ_k+1 on the specific sample points G_θ_k^(k)(X^i)_i=1^N_k+1. Definition of N_k+1≥ N(g, ℋ_k+1, ϵ, δ_k/k) again allow derivative approximation transfer so that with probability at least 1-δ_k/k:

    _XG_θ_k+1(G_θ_k^(k)(X)) - G(G_θ_k^(k)(X)) < ϵ


Now consider the probability subspace 𝒮 under which 3 inequalities <ref>, <ref>, and <ref> hold. The subspace 𝒮 has the probability measure of at least 1 - (δ_k + (k-1)δ_k/k + δ_k/k) = 1 - 2δ_k = 1 - (δ_k+1-δ_k), and under 𝒮, we have:

    _XG_θ_k+1^(k+1)(X) - G^(k+1)(X)≤_XG_θ_k+1(G_θ_k+1^(k)(X)) - G_θ_k+1(G_θ_k^(k)(X))
         +_XG_θ_k+1(G_θ_k^(k)(X)) - G(G_θ_k^(k)(X)) + _XG(G_θ_k^(k)(X)) - G(G^(k)(X))
       ≤ L_XG_θ_k+1^(k)(X) - G_θ_k^(k)(X)
     + ϵ + L_XG_θ_k^(k)(X) - G^(k)(X)
       ≤ Lα_k + ϵ + Lϵ_k = ϵ_k+1


Merging this inequality with probability subspace where the inequality in <ref> holds leads to the estimate on the final step for stage k+1 for the induction step.


Together with the general pointwise estimates for DPO algorithm in <ref>, <ref> and <ref> allow us to explicitly state the number of training episodes required for two scenarios considered in this work: one work with general neural network approximators and the other with more restricted (weakly convex and linear bounded) difference functions:

 
In <ref> (DPO), suppose we are given fixed step size and fixed number of steps per episode H. Further assume that for all k ∈1, H-1, ℋ_k is the same everywhere and is the hypothesis space ℋ consisting of neural network approximators with bounded weights and biases. Then with the sequence of numbers of training episodes N_k = O(ϵ^-(2d+4)), the pointwise estimates <ref> hold. 


 
Again, in <ref> (DPO), suppose we are given fixed step size and fixed number of steps per episode H. Suppose ℋ_k is a hypothesis subspace consisting of h ∈ℋ_k so that h - g and h - g_θ_k-1 are both p-weakly convex and linearly bounded. Then with the sequence of numbers of training episodes N_k = O(ϵ^-6), we obtain the pointwise estimates <ref>. 



In order to justify the Lipschitz constant L assumption, we note that the dynamics function G and its approximator often have the discretized dynamics form: G(x) = x + Δ F(X), where Δ is the step size, and F is a bounded function. In this case, even when the step size is infinitely small, <ref> (DPO)'s training converges with reasonable numbers of training episodes.


The Lipschitz constant L of G in this case is bounded by 1 + CΔ for some constant C > 0. By scaling, WLOG, assume that for finite-time horizon problem, the terminal time T = 1 so that the number of steps is m=1/Δ. Hence, for N_k = O(1/Δ^2p), ϵ=O(Δ^p) for p > 2, the error bounds in <ref> are upper-bounded by:

    G_θ_k^(k)(X) - G^(k)(X)≤kL^kϵ/L-1≤1/Δ(1+C/N_k)^N_kO(Δ^p) 1/Δ≤ e^C O(Δ^p-2) → 0

for k ≤ H-1, as Δ→ 0.




 §.§ Regret bound analysis

In this section, we define explicitly the regret mentioned in <ref> and then derive two regret bounds for DPO algorithm (<ref>). Suppose K episodes are used during the training process and suppose  a policy π^k is applied at the beginning of the k-th episode with the starting state s^k (sampled from adversary distribution ρ_0) for k ∈1, K. If we focus on the total number of training episodes used, and assuming that number of steps H is fixed. Then, the quantity Regret is defined as the following function of the number of episodes K:

    Regret(K) = ∑_k=1^K (V(s^k) - V_π^k(s^k))

We now derive an upper bound on Regret(K) via <ref>:
 
Suppose that number of steps per episode H is fixed and relatively small. If in <ref> (DPO), the number of training samples N_k have the scale of O(ϵ^-μ), regret bound for is upper-bounded by O(K^(μ-1)/μ). As a result, for the general case in <ref>, we obtain a regret bound of O(K^(2d+3)/(2d+4)). For the special cases with restricted hypothesis space in <ref> the regret bound is approximately O(K^5/6). 


For a fixed H, <ref> gives us the estimate _X[G_θ_k^(j)(X) - G(X)] ≈ O(ϵ) for the state-action pair X at step j between (N_1 + ⋯ + N_k-1 + 1)^th and (N_1 + ⋯ + N_k)^th episodes during stage k. Assuming a mild Lipschitz condition on reward function, the gap between the optimal reward and the reward obtained from the learned policy at step j during these episodes is also approximately O(ϵ). Summing these uniform bounds over all j and over all episodes gives: 

    Regret(K) ≤ H(N_1 + ⋯ + N_H-1) Cϵ = CHKϵ

Since N_k = O(ϵ^-μ), K = HO(ϵ^-μ). With a fixed H, ϵ = O(K^-1/μ). Hence, <ref> leads to Regret(K) ≤ K O(K^-1/μ) = O(K^(μ-1)/μ) as desired.



§ EXPERIMENTS ON LAGRANGIAN REWARDS PROBLEMS

We consider 3 specific problems under the form of <ref>: materials deformation, topological materials deformation, and molecular dynamics.



 §.§ Materials deformation

In materials engineering <cit.>, the process of deforming a raw figure of a particular material such as wood, metal, or plastic into a desired shape is essential. As each deformation is costly in terms of energy usage, we demand a minimum number of deformations. In abstract terms, given a particular shape Γ_0 in the space of shapes 𝒮(shape), one wants to deform Γ_0 as fast as possible to get to a desired shape Γ^*, where the desired property is characterized by minimizing a functional cost 𝒞: 𝒮(shape) →ℝ. Assuming that the deformation/hamming machine can only incrementally deform the current object's outer part, we proceed to formulate this problem in terms of <ref> as follows: the state s for a 2D object Γ is the vector of discrete points on its boundary. The object Γ = Γ(s) is then constructed from s by the 1D cubic spline approximation so that ℱ(s) = 𝒞(Γ(s)). The action a is defined as the hamming machine's modification at the boundary's discrete points. For experiments, we choose the functional cost to be the dimensional-less quantity: 𝒞(Γ)=Peri(Γ)√(Area(Γ)), where Peri and Area are perimeter and area of Γ. The distribution ρ_0 gives a random polygon at the start of each episode.



 §.§ Topological materials deformation

Rather than restricting the deformation to the boundary only, one allows the hamming machine to modify any point inside and even change the topological structure of the object. For this new setting, the state s of Γ is now a discretization vector of Γ's level-set function on a fixed grid. The original Γ can then be constructed via a bicubic spline approximation <cit.>. The action a is the modification of the level-set function at grid points. The starting distribution ρ_0 is an uniform distribution on grid points.





 §.§ Molecular dynamics

Molecule's configuration is the 3D spatial arrangement between its atoms and is essential to the molecule's functioning. For a given molecule, we're interested in capturing its movement from a random configuration to a stable configuration, where the potential energy is minimized. With respect to the formulation <ref>, the state s of each molecular configuration is a vector of its dihedral angles denoted by (ϕ_j, ψ_j)_j, where j are the angle indices. Given these angles, we can recover the atom coordinates (x_i)_i=1^N = X((ϕ_j, ψ_j)_j) as a function of dihedral angles. Such coordinates allow us to calculate the functional ℱ(s) = ℱ((ϕ_j, ψ_j)_j) = U(X((ϕ_j, ψ_j)_j)) = U((x_i)_i=1^N) where U is the energy function on the atom coordinates given by PyRosetta Python package <cit.>. The action a is the modification of the molecule's dihedral angles. The choice of the dihedral angles is due to the rigidity property of molecule. For this experiment, we consider the octa-alanine with the primary sequence consisting of eights amino acids A. The starting distribution ρ_0 is purposely chosen to be an uniform distribution over a very small compact intervals around zero to see how the agent performs with limited exploration from the start.



 §.§ Benchmark results

Models:  We run 4 models including <ref> (DPO), Trust Region Policy Optimization (TRPO) <cit.>, Proximal policy optimization (PPO) <cit.>, and Soft Actor-Critic (SAC) <cit.> on 3 problems with RL formulation in <ref> as outlined above. We additionally run 3 benchmark models TRPO, PPO, SAC on these 3 problems with the straightforward negative energy reward r(s, a) = -ℱ(s). We call the 3 trained models on this type of straightforward rewards: Straight TRPO, Straight PPO, and Straight SAC, bringing in a total of 6 benchmark models. Implementations of benchmark models come from Stable-Baselines3 package <cit.>.


Parameters:  For each case, we use 200 test episodes with normalized time horizon [0, 1] (T = 1). The number of episodes, the number of time steps, and the magnification factor γ are outlined in <ref>. Here γ is specific to each problem and is independent of the learning methods. All training and testing have been conducted on Nvidia A100 GPU.



Metrics:  The main metrics are the evaluation costs, which are the values of ℱ evaluated along the test episodes. The goal is to arrive at the minimum ℱ-values in a fast and smooth manner.


Results:  From <ref> and <ref>, DPO (<ref>) performs the best across 3 problems. <ref> showing evaluation costs with uncertainty at each time step further demonstrates that DPO yields more desirable optimal paths.




Related model-based reinforcement learning algorithms.  Several model-based RL algorithms <cit.> are not applicable here because:

    
  * Policy search with back-propagation through time requires more information than what the system ℬ could provide. For instance, SVG <cit.> and iLQR <cit.> require explicit forms of ℱ's derivatives, while PILCO <cit.> needs the convolution of ℱ with a Gaussian kernel. For simple problems such as Acrobot or Ant-walk, such information can be hand-calculated. But this is impossible for the 3 complex problems considered in this section. For the first two, ℱ is a combination of a geometric functional together with cubic spline approximation operator. For the third problem, the energy functional contains molecular statistics unknown constants.
    
    
  * Shooting algorithms <cit.> require the agent to perform replanning at every time step. They need to execute the action sequence starting from an arbitrary time t, while <ref> only allows t = 0.




§ CONCLUSION

In this paper, we develop a new approach called differential reinforcement learning focusing on the quality of the individual data points along the learned path. The duality between our problem and the original RL formulation is demonstrated as a generalization of the classical duality between Lagrangian and Hamiltonian mechanics. We propose a simple yet effective algorithm called differential policy optimization (DPO, see <ref>) as a concrete solution to our approach. Theoretically, DPO not only possesses pointwise convergence guarantees but also achieves regret bounds comparable to current literature. Experimentally, DPO demonstrates promising performance on 3 physics-based problems with complex reward systems where the learning opportunities are limited under a small number of trajectories and number of time steps per episode. In the future, we aim to take the next few steps in promoting practical RL algorithms for continuous state and action spaces. These include: improving solution algorithms for differential RL both theoretically and empirically, and drawing more connections between our new approach to current state-of-the-art RL methods.


Author acknowledgement: This research was supported in part by a grant from the NIH DK129979, in part from the Peter O’Donnell Foundation, the Michael J. Fox Foundation, Jim Holland-Backcountry Foundation, and in part from a grant from the Army Research Office accomplished under Cooperative Agreement Number W911NF-19-2-0333.


Author Contributions: This research work (Differential RL and DPO) is developed by Minh Nguyen and improved by Chandrajit Bajaj and Minh Nguyen. Alphabetical order is currently used to order author names. The implementation is done by Minh Nguyen and is available at https://github.com/mpnguyen2/dpohttps://github.com/mpnguyen2/dpo.







splncs04





§ BASIC ALGORITHMIC LEARNING THEORY UNDER I.I.D SETTING

We first introduce basic learning theory result on independent identically distributed (i.i.d) samples. Proofs for all lemmas in this section can be found at <cit.> and can also be found in any modern learning theory literature.


Notations: Throughout this section, 𝒳 be the feature space, 𝒴 be the label space, 𝒵 = 𝒳×𝒴. Let ℋ be a hypothesis space consisting of hypothesis h: 𝒳→𝒴. Let l: ℋ×𝒵→ℝ_+ be a loss function on labeled sample z = (x, y) for x ∈𝒳 and y ∈𝒴. Our loss function will have the following particular form: l(h, (x, y)) = ϕ(h(x), y) for a given associated function ϕ: 𝒴×𝒴→𝐑_+. We use capital letters to represent random variables. We also use the following notations to define the sets: ℋ∘Z_1, ⋯, Z_n: = (h(Z_1), ⋯ h(Z_n) ∈ℝ^n, h ∈ℋ
and ℒ∘Z_1, ⋯, Z_n: = l(h, Z_1), ⋯ l(h, Z_n) ∈ℝ^n, h ∈ℋ. Also define e(h) to be the average loss over new test data e(h):=_Z[l(h, Z)], and E(h) be the empirical loss over n ∈ℤ_+ i.i.d samples Z_1, ⋯, Z_n: E(h):= 1/n∑_i=1^n l(h, Z_i).



The Rademacher complexity of a set 𝒯 is defined as:

    Rad(𝒯) = [sup_t ∈𝒯1/n∑_i=1^n B_i t_i]

for Bernoulli random variables B_i ∈-1, 1.



Suppose ϕ(., y) is γ-Lipschitz for any y ∈𝒴 for some γ > 0. Then:

    [sup_h ∈ℋe(h) - E(h)]    ≤ 2 [Rad(ℒ∘Z_1, ⋯, Z_n)] 
       ≤ 2γ[Rad(ℋ∘Z_1, ⋯, Z_n)]




For a hypothesis space ℋ, a loss function l bounded in the interval [0, c], and for n i.i.d labeled samples Z_1, ⋯, Z_n, with probability of at least 1 - δ, the following bound holds:

    sup_h ∈ℋ (e(h) - E(h)) < 4 Rad(ℒ∘Z_1, ⋯, Z_n) + c√(2log(1/δ)/n)




For hypothesis space ℋ consisting of (regularized) neural network approximators with weights and biases bounded by a constant, there exists a certain constant C_1, C_2 > 0 so that for a set of n random variables Z_1, ⋯, Z_n: 

    Rad(ℋ∘Z_1, ⋯, Z_n) ≤1/√(n)(C_1 + C_2√(log d))



Throughout this paper, we assume that the optimization error can be reduced to nearly 0, so that E(h) ≈ 0 if the hypothesis space ℋ contains function to be learned. From <ref>, the average estimation error generally scales with c√(2log(1/δ)n)


§ PROOFS OF LEMMAS IN <REF>

First of all, note that the continuous density ρ_0 for starting states on Ω is bounded below by a constant.  
*

Suppose we're given ϵ > 0. Take ϵ_1 > 0 so that 16C ϵ_1 < ϵ/(2d), and ϵ_2 = 0.5. Now take n ∈ℕ large enough so that C_1√(log(1/(δ/(3d))n)< ϵ/(2d) for an appropriate constant C_1. Then take m ∈ℕ large enough so that (1-c_1ϵ_2ϵ_1^d)^m < δ/(3n), where c_1 is an appropriate geometric constants (see paragraphs below). Let M = m + n.


Train a neural network function h ∈ℋ to approximate the target function g on N samples, where N is given by:

    N = log((6M)/δ)/(Cϵ_1^2δ/(6M))^2 or √(log(1/(δ/6M))/N) = C ϵ_1^2 δ/6M = N(g, ℋ, ϵ, δ)


Before going to the main proof, we dissect N to get its asymptotic rate in terms of ϵ and d. First of all ϵ_1 = O(ϵ/2d). Next, n ≈log(1/δ)/(ϵ/(2d))^2, and m ≈log(3n/δ) ϵ_1^-d. Hence by ignoring log terms, polynomial terms in d, and quadratic factor of δ, M = m + n is approximately O(ϵ^-d).  As a result, N ≈ϵ^-(2d+4). 


Choose a set S = X_1, ⋯, X_m, Y_1, ⋯, Y_n consisting of M = m + n random samples from distribution ρ_0 which are independent with g: m samples X_1, ⋯, X_m and n sample Y_1, ⋯, Y_n.


We apply <ref> to ℋ with i.i.d labeled samples Z = (X, g(X)) and loss function l with the associated ϕ(y, ŷ) = |y - ŷ|, which is 1-Lipschitz for a fixed ŷ. In this case, from <ref>, for probability of at least 1 - δ/(6M), _U[|h(U)-g(U)|] < Cϵ_1^2δ/(6M) for the random variable U ∈ S. As a result, by Markov inequality, with probability at least 1 - δ/(6M) -δ/(6M) = 1 - δ/(3M), |h(U)-g(U)| < C ϵ_1^2 for each U ∈X_1, ⋯, X_m, Y_1, ⋯, Y_n. Hence there exists a probability subspace Γ with probability at least 1 - δ/(3M) * M = 1-δ/3 so that |h(U)-g(U)| < C ϵ_1^2 for all U ∈ S.


The probability that a particular sample (random variable) X_i is in the hyper-cone with an conic angle difference of ϵ_2 surrounding the direction (∇ h(Y) - ∇ g(Y)) of the hyper-spherical ϵ_1-circular neighborhood of Y is at least c_1ϵ_2ϵ_1^d. Here a ϵ_1-circular neighborhood here include points with radius sizes between ϵ_1/2 and ϵ_1. The probability that no m samples is in this cone is at most (1 - c_1ϵ_2ϵ_1^d)^m < δ/(3n). As a result, on Γ except a subspace with probability less than δ/(3n), there exists k (depends on both Y and X_1, ⋯, X_m) so that:

    ϵ_1/2    < X_k - Y < ϵ_1 
    
    (∇ h(Y) - ∇ g(Y)) · (X_k - Y)    > (1-ϵ_2) ∇ h(Y) - ∇ g(Y)X_k - Y


Second-order Taylor expansion for f and g at each Y yields:

    h(X_k)    = h(Y) + ∇ h(Y) (X_k-Y) + 1/2X_k-Y^2 U_h(X_k, Y) 
    
    g(X_k)    = g(Y) + ∇ g(Y) (X_k-Y) + 1/2X_k-Y^2 U_g(X_k, Y)

where U_h and U_g are the second derivative terms of h and g in respectively, and are bounded by C. As a result:

    (1-ϵ_2) ∇ h(Y) - ∇ g(Y)(ϵ_1/2) 
       < (1-ϵ_2) ∇ h(Y) - ∇ g(Y)X_k - Y
       < (∇ h(Y) - ∇ g(Y)) · (X_k - Y) 
       < |h(Y) - g(Y)| + |h(X_k) - g(X_k)| + 2C X_k-Y^2 
       < 2Cϵ_1^2 + 2C X_k-Y^2 < 4C ϵ_1^2

Thus, ∇ h(Y) - ∇ g(Y) < 8C ϵ_1/(1-ϵ_2) = 16Cϵ_1. Therefore, on the subspace Γ_0 with probability of at least 1 - δ/3 - n*(δ/(3n)) = 1 -  (2δ)/3, ∇ h(Y) - ∇ g(Y) < 16Cϵ_1 for all samples Y ∈Y_1, ⋯, Y_n.


In order to prove that _X∇ h(X) - ∇ g(X) < ϵ and thus finishing the proof for N(g, ℋ, ϵ, δ) = O(ϵ^-(2d+4)), we only need to prove bounds on individual components of ∇ h(X) - ∇ g(X): 

    _X(∇ h(X) - ∇ g(X))_s < ϵ/d

where x_s is the s^th component of a vector x ∈ℝ^d. To this end, if we have a bound of order O(1/√(n)) for the Rademacher complexity of the hypothesis space ℋ_s = ∂_s h, h ∈ℋ for s ∈1, d, then we can invoke the corollary <ref> on ℋ_s to get:

    _X(∇ h(X) - ∇ g(X))_s   < 1/n∑_i = 1^n (∇ h(Y_i) - ∇ g(Y_i))_s + ϵ/(2d) 
       < ϵ/(2d) + ϵ/(2d) = ϵ/d

on Γ_0 except a set of probability of at most δ/(3d). Then we finish the proof of <ref> by summing all inequalities over d components.


For the final argument to get a Rademacher complexity bound on ℋ_s, we use the chain rule to express each element of ℋ_s as h_1⋯ h_R, where R is the fixed number of layers in ℋ's neural network architecture. Here each h_i can be expressed as the composition of Lipschitz (activation) functions and linear functions alternatively with bounded weight and bias. Those elements h_i then form another neural network hypothesis space. By invoking <ref>, we obtain a bound of order O(1/√(n)) on individual h_i's. To connect these h_i's, we express the product h_1⋯ h_R as:

    h_1⋯ h_R    = ∏_i=1^R (h_i + D - D) = ∑_W ⊆ [R] (-D)^R-|W|∏_i ∈ W (h_i + D)) 
       = ∑_W ⊆ [R] (-D)^R-|W|exp(∑_i ∈ Wlog(h_i + D)))

Here D is a constant large enough to make the log function well-defined and to make Lipschitz constant of the log function bounded above by another constant. Now we use the simple bounds on Rad(𝒯+𝒯') and Rad(f∘𝒯) for some sets 𝒯 and 𝒯' and some Lipschitz function f to derive the Rademacher complexity bound of the same order O(1/√(n)).


Another simpler way to achieve a similar result is to upper-bound ∇ h(X) - ∇ g(X) by ∇ h(Y_i) - ∇ g(Y_i) + 2C X-Y_i and use the probability subspace in which one of the Y_i is close enough to X. However, we need a similar argument for <ref>, so we moved forward with the current approach.


*

Suppose we're given ϵ > 0. Take ϵ_1 = ϵ^α with α = 1/d, and set ϵ_2 = 1/2. Now choose m ∈ℕ large enough so that (1-c_1ϵ_2ϵ_1^d)^m < ϵ^2, where c_1 is an appropriate geometric constant. From here, we can see that m ≈ O(ϵ_1^-d) = O(ϵ).

Choose N ≈ O(ϵ^-6) ∈ℕ:

    N = O(log(d/δ)/ϵ^6)  so that √(log(1/(δ/d))/N)≈ O(ϵ^3)


Train a neural network function g to approximate f on N samples Y_1, ⋯, Y_N so that g(Y_i) ≈ f(Y_i) for i ∈1, N. We prove that this N is large enough to allow derivative approximation transfer.


Choose β_k = k/d and δ_k = k/d for k ∈0, d. We prove by induction on k ∈0, d that there exists a subspace of probability at least 1 - δ_k so that ∇ f(Y) - ∇ g(Y) < C_0ϵ^β_k for all Y ∈Y_1, ⋯, Y_n and for some constant C_0. 


For k = 0, the bound is trivial. For the induction step from k to k+1, we first consider the loss function l_k of the form l_k(h, (x, y)) = l_k(h, (x, g(x))) = ϕ_k(h(x), g(x)), where ϕ_k(y, ŷ) = clip(|y - ŷ|, 0, Cϵ^α + β_k)^d. Here clip is to denote a clip function, and, in this case, is obviously a Lipschitz function with Lipschistz constant 1. First, the loss function l_k is bounded by (Cϵ^α + β_k)^d. Now note the following simple inequality:

    |a^d - b^d| = |a-b||∑_k=0^d-1 a^k b^d-1-k| < |a-b| d (Cϵ^α + β_k)^d-1

for a, b < Cϵ^α + β_k. The inequality shows that ϕ_k has the Lipschitz constant bounded by d(Cϵ^α + β_k)^d-1. We are now ready to go the main step of the induction step.


Condition on Y, we repeat the same argument in <ref>'s proof to show that except for a subspace with probability less than ϵ^2, there exists j ∈1, m (depends on both Y and X_1, ⋯, X_m) so that:

    ϵ_1/2    < X_j - Y < ϵ_1 = ϵ^α
    
    (∇ h(Y) - ∇ g(Y)) · (X_j - Y)    > (1-ϵ_2) ∇ h(Y) - ∇ g(Y)X_j - Y

Under this subspace, because h-g is linearly bounded, 

    |h(X_j)-g(X_j)| ≤ C ∇ h(Y) - ∇ g(Y)Y-X_j < C ϵ^α + β_k

As a result, |h(X_j)-g(X_j)|^d = ϕ_k(h(X_j), g(X_j)). Next, since h - g is p-weakly convex and hence d-weakly convex, we have:

    (1-ϵ_2) ∇ h(Y) - ∇ g(Y)(ϵ_1/2) 
       < (1-ϵ_2) ∇ h(Y) - ∇ g(Y)X_j - Y
       ≤ (∇ h(Y) - ∇ g(Y)) · (X_j - Y) 
       ≤ |h(X_j)-g(X_j)| + |h(Y)-g(Y)| + C_1 X_j-Y^d 
       < |h(Y) - g(Y)| + |h(X_j) - g(X_j)| + 2C_1 ϵ^α d
       = |h(X_j) - g(X_j)| + 2C_1 ϵ^2
       = ϕ_k(h(X_j), g(X_j))^1/d + 2C_1 ϵ^2
       ≤(∑_i=1^m ϕ_k(h(X_i), g(X_i)))^1/d + 2C_1 ϵ^2


By taking expectation over X_1, ⋯, X_m, we obtain

    (1-ϵ_2) ∇ h(Y) - ∇ g(Y)(ϵ_1/2)
       < (1-ϵ^2)((_X_1, ⋯, X_m[∑_i=1^m ϕ_k(h(X_i), g(X_i))])^1/d + 2C_1 ϵ^2) + C_2 ϵ^2
       < (m _X[ϕ_k(h(X), g(X))])^1/d + C_3 ϵ^2
       = (m _X[l_k(h, (X, g(X)))])^1/d + C_3 ϵ^2

for appropriate constants C_2, C_3 > 0.


By <ref> and <ref> on the Lipschitz loss function l_k bounded by (Cϵ^α + β_k)^d with Lipschitz constant d(Cϵ^α + β_k)^d-1, with probabilty of at least 1-δ/d, we can continue the sequence of upper-bounds:

    (1-ϵ_2) ∇ h(Y) - ∇ g(Y)(ϵ_1/2)
       ≤ C_4 (ϵ^-1 (ϵ^α + β_k)^(d-1)ϵ^3)^1/d + C_3 ϵ^2 
       ≤ C_5 ϵ^α + β_k + 1/d = C_5 ϵ^α + β_k+1

for appropriate constants C_4, C_5 > 0.


Hence, we finish the induction step because except on a subspace with probability at most 1 - δ_k - δ/d = δ_k+1, the inequality for induction hypothesis at (k+1) holds. For β_d = 1, we obtain ∇ h(Y) - ∇ g(Y) < ϵ for all Y ∈Y_1, ⋯, Y_N. By repeating the argument in <ref>'s proof, we obtain the expected bound 

    _X [∇ h(Y) - ∇ g(Y)] < C_6 ϵ

for an appropriate constant C_6 > 0. By rescaling by a constant, we showed that N(g, ℋ, ϵ, δ) ≈ O(ϵ^-6)


*

First, α_j = (L^j-1 + ⋯ + 1) ϵ. Hence,

    ϵ_k/L^k   ≤L(L^k-2+⋯+1) + 1/L^kϵ + ϵ_k-1/L^k-1
       = L^k-1/L^k(L-1)ϵ +  ϵ_k-1/L^k-1
    < ϵ/L-1 + ϵ_k-1/L^k-1

Hence, by simple induction, ϵ_kL^k < (k-1)ϵL-1 + ϵ_1L < kϵL-1. As a result, ϵ_k < kL^kϵL-1 as desired.


