Year: 2021

Abstract:
Adversarial imitation learning techniques are based on modeling statistical divergences using agent and expert demonstration data. However, unbiased minimization of these divergences is not usually guaranteed due to the geometry of the underlying space. Furthermore, when the size of demonstrations is not sufficient, estimated reward functions from the discriminative signals become uncertain and fail to give informative feedback. Instead of formulating a global cost at once, we consider reward functions as an iterative sequence in a proximal method. In this paper, we show that rewards dervied by mirror descent ensures minimization of a Bregman divergence in terms of a rigorous regret bound of O(1/T ) for a particular condition of step sizes {ηt}t=1. The resulting mirror descent adversarial inverse reinforcement learning (MD-AIRL) algorithm gradually advances a parameterized reward function in an associated reward space, and the sequence of such functions provides optimization targets for the policy space. We empirically validate our method in discrete and continuous benchmarks and show that MD-AIRL outperforms previous methods in various settings.

heading: 1 INTRODUCTION
text: Inverse reinforcement learning (IRL) is an algorithm of learning ground-truth rewards from demonstrations of an expert that acts optimally with respect to an unknown reward function. Traditional IRL approaches (Ng & Russell, 2000; Abbeel & Ng, 2004; Ziebart et al., 2008) solve the imitation problem based on iterative algorithms, alternating between the reward estimation process and a reinforcement learning (RL) algorithm. Recent adversarial imitation learning (AIL) studies (Ho & Ermon, 2016; Fu et al., 2017) focus on formulating a statistical divergence minimization with a fine-tuned representation of the target expert probability distribution (Ghasemipour et al., 2020).
One method to approximate the distance between probability distributions is using Bregman divergences, a family of metric-like functions induced by a strongly convex function. Jeon et al. (2020) generalized the AIL framework to an optimization problem with respect to a Bregman divergence between policy distributions where the method shares on the essence of regularized Markov decision processes (Geist et al., 2019). Through the lens of differential geometries, this implies that a policy distribution and a reward function can be associated with geometric constraints specified by a convex regularizer (Shima, 2007). Apparently, the limitation of regularized IRL studies naturally comes from a geometric property in which minimizing a Bregman divergence does not guarantee unbiased progression due to the constraints of the underlying space (Butnariu & Resmerita, 2006).
The success of machine learning has been the key to learning such divergences with high precision. When the demonstration size is sufficiently large, the discriminator approximation of AIL can model the entire expert policy distribution internally using nonlinear approximators such as neural networks (Finn et al., 2016; Fu et al., 2017). However, there are challenging tasks that require long action sequences to be solved. When states can only be rarely visited by imitation learning agents, the uncertainty of estimated expert policy often substantially affects the discrimination quality. Therefore, the development of a robust IRL method that is tolerant to unreliable discriminative signals induced by imperfect demonstrations is necessary.
The reasonable question is how to learn with cost estimations that are not precisely determined throughout the learning process. There have been considerable achievements in dealing with temporal costs in the online learning domain (Fiat & Woeginger, 1998; Hazan, 2019). An online learning
algorithm predicts a sequence of parameters {wt}Tt=1 ⊂ W that is incurred by a cost functionFt(·) at each iteration. The most ordinary approach in online learning is stochastic gradient descent (SGD): wt+1 = wt − ηt∇Ft(wt). SGD is a desirable algorithm when the setW is the Euclidean geometry since it can provide unbiased progression for minimizing Ft (Raskutti & Mukherjee, 2015). However, policies appear in manifolds of probability distributions; a gradient may not be the direction of the steepest descent in this case due to geometric constraints (Amari, 1998; 2016).
In this paper, we propose an iterative imitation learning framework that interprets the RL-IRL scheme similar to proximal optimization methods (Boyd et al., 2004; Amari, 2016). We identify two issues in the AIL studies that have characteristics of unconstrained updates: (1) a divergence does not guarantee informative global directions to match the expert policy due to the geometry, and (2) a representation of divergence often cannot be accurately obtained due to insufficient data.
Our method is motivated by a template of optimization algorithms specified by a convex function, called mirror descent (MD; Nemirovsky & Yudin 1983). For parameters sequences, cost functions, and step sizes, an MD update for a strongly convex function Ω is derived as
∇Ω(wt+1) = ∇Ω(wt)− ηt∇Ft(wt), (1) where ∇Ω is a bijective transformation that links the primal space and the dual space of gradients. Different from these standard MD formulations, our methodology draws a sequence of functions on a space formulated by a regularized reward operator which is an alternative to the dual space. Consequently, the reward functions are projected optimization targets for the space of policies.
Our Contributions. We propose a novel IRL algorithm which facilitates the agent in robustly imitating the expert. Our work is complementary to previous regularized IRL studies; we introduce NEW R1, R2 a geometric perspective for optimizing rewards and derive solutions with theoretical guarantees. • Instead of a monolithic estimation process of a global solution, we derive a sequence of reward
functions that provides local optimization targets for the space of policies (Section 3). • We prove that rewards derived by an MD algorithm guarantee convergent divergence minimization
performance along with a rigorous regret bound (Section 4). • We propose mirror descent adversarial inverse reinforcement learning (MD-AIRL), a novel IRL
algorithm that can be easily implemented on top of the standard AIL framework (Section 5). • We validate the outperforming performance of MD-AIRL on benchmarks with large discrete ac-
tion spaces and continuous action spaces (Section 6).

heading: 2 PRELIMINARIES
text: Notation. For finite setsX and Y , we define Y X as a set of functions fromX to Y . ∆X (∆YX ) is a set of (conditional) probabilities over X (conditioned on Y ). For a function f ∈ RS×A and a policy π ∈ ∆SA, fs and πs denote shorthand notation of f(s, ·) and π(·|s), respectively. A Markov decision process (MDP) is defined as a tuple (S,A, P, r, γ) with the state space S, the action space A, the transition kernel P ∈ ∆S×AS , the reward function r ∈ RS×A and the discount factor γ ∈ [0, 1). A Bregman divergence with respect to a convex function Ω for a state s is defined as DΩ(πs∥π̂s) = Ω(πs)− Ω(π̂s)− ⟨∇Ω(π̂s), πs − π̂s⟩A. We write Ω(π) = [Ω(πs)]s∈S and similarly for∇Ω. Regularized RL & IRL. We consider the RL-IRL framework in regularized MDPs (Geist et al., 2019), where the policy is optimized along with a causal convex regularizer. The objective is to find π which maximizes the expected discounted sum of rewards with a strongly convex regularizer Ω:
maximize π∈∆SA
JΩ(π, r) := Eπ [∑∞ i=0 γi { r(si, ai)− Ω ( π(·|si) )}] , (2)
where the subscript π on the expectation denotes samples generated from the MDP and policy π.
Consider the convex conjugate of qsE = ∇Ω(πsE): Ω∗(qsE) = maxπs∈∆A⟨πs, qsE⟩A − Ω(πs). Differentiating both sides with respect to qsE , the gradient of conjugate∇Ω∗ maps qE to the unique optimal policy πE . Let Π be a bounded, open set for ∆SA ⊂ Π and R : Π → RS×A be a reward operator which maps a policy to a reward function. The output r̂ = R(πE) is a representation of the expert’s behavior; but finding such operator is an ill-posed problem because every function r̂ that makes the state-action value function qE is a valid solution of regularized IRL (Geist et al., 2019). Recently, Jeon et al. (2020) proposed a tractable solution, which we refer to the function as the regularized reward function. We rewrite the previous notation by defining the regularized reward operator ΨΩ.
Definition 1 (Regularized reward operators). Define the regularized reward operator ΨΩ : π 7→ ψπ ψπ(s, a) := Ω ′(s, a;π)− 〈 πs,∇Ω(πs) 〉 A +Ω(π s) ∀s ∈ S, a ∈ A, (3)
for Ω′(s, ·;π) := ∇Ω(πs) = [∇pΩ(p)]p=π(·|s).
By using the operator in discriminative networks, Jeon et al. (2020) proposed regularized adversarial IRL (RAIRL) as a generalization of AIL in terms of minimizing a specific Bregman divergence.
Mirror Descent. Let w be a parameter on a set W and Ft : W → R be a convex cost function from a class of functions F . Replacing the L2 proximity term of proximal gradient descent with the Bregman divergence (Gutman & Peña, 2018), the proximal form of the MD update is given as
wMDt+1 = argmin w∈W
〈 ∇Ft(wMDt ), w − wMDt 〉 W + αtDΩ ( w ∥∥wMDt ), (4)
where αt := 1/ηt denotes an inverse of the step size ηt. MD was developed as a generalization of SGD, where the local geometry is specified by a Bregman divergence (Gunasekar et al., 2020).

heading: 3 ITERATIVE RL-IRL AS A PROXIMAL OPTIMIZATION METHOD
text: We interpret the RL-IRL framework as a variant of proximal optimization methods and consider two sequences {πt}∞t=1 and {ψt}∞t=1 that denote the learning policies and reward functions, respectively. Associated reward function. An updated point in MD can be uniquely projected to the desired space (∆SA in our case) using a Bregman projection operator PΩ that locates iterative points to the feasible region, i.e. πt+1 = PΩ(π̃t+1) := argminπ∈∆SA [DΩ(π∥π̃t+1)]s∈S for π̃t+1 ∈ Π. To avoid these computations, we preemptively constrain a reward operator to satisfy the projection invariance with respect to regularized RL process of the associated state-action value function q̂: PΩ(∇Ω∗(q̂)) = ∇Ω∗(q̂). According to Lemma 1 of Jeon et al. (2020), a regularized reward function ψπ can replace q̂, since the induced Bregman divergence allows the learning in a greedy manner. FIX R3 The projection invariance of the operator ΨΩ can be shown by a bijective relation with ∇Ω∗ in the space of ∆SA by the following lemma. Lemma 1 (Natural isomorphism). Let ψ ∈ ΨΩ(∆SA) for ΨΩ(Π) = {ψ | ψ(s, a) = ψπ(s, a), ∀s∈ S, a∈A, π∈Π }. Then,∇Ω∗(ψ) is unique with respect to ψ, and for every π = ∇Ω∗(ψ), π ∈ ∆SA.
Policy Space
∆SA
Π
πt
πt+1
Reward Space
ψt
ψt+1
ΨΩ
∇Ω∗
FIX R3
As a result, we consider an updated reward function ψt+1 ∈ ΨΩ(∆SA) as a projected target of MD, FIX R3 which can be associated by parameterization of ∆SA. For instance, a bijective mapping of ψt+1 can NEW R3 represent a softmax policy for a discrete space, or a Gaussian policy for a continuous space, using its parameters. The subsequent RL process at t-th step becomes finding the next iteration πt+1 by NEW R3
FIX R3maximize π∈∆SA
JΩ(π, ψt+1) = Eπ [∑∞
i=0 γiDΩ
( πsi ∥∥πsit+1)], πt+1 = ∇Ω∗(ψt+1). (5)
The equation shows that an arbitrary regularized RL algorithm with the regularizer Ω forms a cumulative discounted sum of Bregman divergences, thus the mapping to πt+1 can be achieved.
Online imitation learning. The necessity of online learning setup comes from observing that our RL-IRL processes do not retain ψE = ΨΩ(πE), a representation of πE , during training. Instead, consider a random process {π̄E,t}∞t=1 where the estimation π̄E,t resides in a closed, convex neighborhood of πE , estimated by a separate estimation algorithm such as maximum likelihood methods. Substituting ψE to ψπ̄E,t in JΩ(π, ψE) turns the RL objective into an online learning problem:
FIX R3minimize π∈∆SA
Ft(π) = −JΩ(π, ψπ̄E,t) = Eπ [∑∞
i=0 γiDΩ
( πsi ∥∥π̄siE,t)]. (6)
For a deeper understanding of our setup, we consider an example in Figure 2. Assume that policies of the learning agent and the expert follow multivariate Gaussian distributions initially at π1 = N ([0, 0]T, I) and πE = N ([5, 3]T,ΣE) for |ΣE| < 1. The policy πt is trained by a cost function DΩ(·∥π̄E,t) with the proximity constraint 1ηDΩ(·∥πt) where the Gaussian policy π̄E,t is also fitted using stochastic samples from πE at each iteration, starting from π̄E,1 = π1. We first observe that the step size constant η affects the training speed in the early phase. Notably, the performance of certain cases exceeds the baselines of DΩ(π̄E,t∥πE) by choosing the step size effectively low as η < 1. This suggests that there is a clear advantage of the online imitation learning setup for unreliable π̄E,t.
MD update rules. Plugging each divergence of the cumulative cost Ft to the template of Equation (4), the optimization process for the subsequent ψt+1 of the MD formulation is derived as1
ψt+1 := ΨΩ(πt+1), π s t+1 = argmin
πs∈∆A 〈 ∇DΩ(πst ∥π̄sE,t)︸ ︷︷ ︸ ∇Ω(πst )−∇Ω(π̄sE,t) , πs − πst 〉 A + αtDΩ(π s∥πst )
= argmin πs∈∆A
DΩ(π s∥π̄sE,t)−DΩ(πs∥πst ) + αtDΩ(πs∥πst )
= argmin πs∈∆A ηtDΩ(π s∥π̄sE,t)︸ ︷︷ ︸
estimated expert
+(1− ηt)DΩ(πs∥πst )︸ ︷︷ ︸ learning agent , ηt = 1/αt, ∀s ∈ S, (7)
where the gradient of DΩ is taken with respect to its first argument πst . The objective of reward NEW R3 learning is analogous to finding an interpolation at each iteration where the point is controlled by the step size ηt. Figure 3 shows that the region of πt (defined by a norm) gradually decreases
1Bregman divergences are generally intractable to be computed, unless the policy is a specific parametric model (e.g. exponential families). See the works of Nielsen & Nock (2011); Jeon et al. (2020) and Appendix B.
when ηt > ηt+1. Note that solving the optimization of Ft requires interaction between πt and the dynamics of the given environment; thus, the RL process in Equation (5) plays the essential role of sequential learning by the value measures of the reward function.

heading: 4 ANALYSES ON STEP SIZES, CONVERGENCE, AND REGRETS
text: Online IRL. We define an online cost ft(πt, τt) which involves a policy and a trajectory as inputs. The trajectory τt = {si}∞i=0 is available at the t-th step by executing the policy πt in the environment. An online cost function for reward ψ̄E,t = ΨΩ(π̄E,t) can be expressed as
FIX R3ft(πt, τt) = ∑∞ i=0 γi {(∑ a πt(a|si)ψ̄E,t(si, a) ) − Ω ( πt(·|si) )} . (8)
We refer the objective of the learning is to find a unique fixed point π∗ ∈ Π that minimizes E[ft(π, τ)], where the expectation is taken over the entire steps (i.e. limt→∞ Eτ1,··· ,τt [f(πt, τt)]). NEW R3 Taking (stepwise) gradient for each π(·|s), π∗ is found by E[∇Ω(π∗(·|s)) − ∇Ω(π̄E,t(·|s))] = 0 when t → ∞, hence ∇Ω(π∗) = limt→∞ E[∇Ω(π̄E,t(·|s))]. The assumption of πE = π∗ allows the particular situation when the estimation algorithm of π̄E,t is actually convergent with t → ∞. This also allows general situations where the estimated expert policy is not stationary; the algorithm finds the unique point according to the expectation of gradient by scheduling the step size ηt.
αt ηt
1− ηt
t
α1
η1
Under the satisfaction of both cases, a sequence that is divergent in its series (e.g. harmonic series∑∞ n=1 1 n ) is suitable as in Figure 4. In our arguments, a policy conditioned by a state is in a Banach space called the Lp space (RA, ∥·∥p), where ∥·∥p denotes a p-norm, where we assume 1 < p ≤ 2 for the convergence in the dual Lq space (1/p + 1/q = 1). The proofs are in Appendix A. Theorem 1 (Stepsize considerations). Let Ω be strongly convex, ∇Ω be Lipschitz continuous, and the associated Bregman divergences are bounded. Assume infπ∈∆SA E[ft(π, τ)] > 0. Then limT→∞ Eτ1,···,τT [∑∞ i=0DΩ(π si ∗ ∥πsiT ) ] = 0 if and only if Eq. (9) is satisfied.
(a) If lim t→∞
ηt = 0, then T ∈ N, n < T , and c > 0 exist s.t. Eτ1,···,τT [ fT (πT , τT ) ] ≥ cT−n .
(b) If the step size takes the form ηt = 4t+1 , then Eτ1,···,τT [∑∞ i=0 γ iDΩ ( πsi∗ ∥∥πsiT )] = O(1/T ).
Theorem 2 addresses the convergence of the algorithm in a specific case when πE can be achieved. Theorem 2 (Convergence in optimal cases). Let Ω be strongly convex,∇Ω be Lipschitz continuous, and the associated Bregman divergences are bounded. Assume π1 ̸= πE and infπ∈∆SA E[ft(π, τ)] = 0. Then, limt→∞ Eτ1,···,τt [ ft(πt, τt) ] = 0 if and only if ∑∞ t=1 ηt = ∞. If ηt ≡ η1, then there exist c1, c2 ∈ (0, 1) such that c1 ·A1 ≤ AT ≤ c2 ·A1, where {At}∞t=1 denote a sequence of At = FIX R3 Eτ1,···,τt [ DΩ(π∗∥πt) ] .
Proposition 1 provides a sufficient condition for the almost certain convergence of an MD algorithm by imposing the stronger condition of step size in Equation (10). Proposition 1 (Convergence in general cases). Assume infπ∈∆SA E[ft(π, τ)] > 0. If the step size sequence satisfies Eq. (10), then we have limt→∞ ∑∞ i=0 γ iDΩ(π si ∗ ∥π si t ) = 0 almost surely.
Regrets. Define the regret at t-th iteration for a sequence of cost functions {ft}t∈N as Rt = 1 t ∑t j=1fj(πj , τj)− infπ∈∆SA { 1 t ∑t j=1fj(π, τj) } . (11) In the optimal case of infπ∈∆SA E[ft(π, τt)] = 0, Rt is bounded by O(1/T ) since ft inherits the property of Bregman divergence so that the infimum is achieved by 0 at πE . By Proposition 1, the updates regarding the policy converge when the step size sequence abides by Equation (10). Therefore, even for the general case of infπ∈∆SA E[ft(π, τ)] > 0, RT is bounded to O(1/T ).
Algorithm 1 Mirror Descent Adversarial Inverse Reinforcement Learning. 1: Input: an expert trajectory dataset (τ∗t )Tt=1, a regularized reward function ψϕ ∈ ΨΩ(∆SA), an
agent policy πθ, an estimate policy πν , a neural network b : S → R, λ, α1, αT ∈ R+. 2: for t← 1 to T do 3: αt ← α1 + αT−α1T−1 (t− 1) 4: Collect rollout trajectories τt by using the agent policy πθ. 5: Optimize b via binary logistic regression for Db to classify τ∗t and τt. 6: Optimize πν via binary logistic regression for Dθ,ν to classify τ∗t and τt. 7: Optimize ψϕ with the MD objective in Equation (7) using αt, τ∗t , τt, πθ, and πν . 8: Train πθ via a regularized actor-critic method to maximize ψλϕ(s, a) with regularizer λΩ(·). 9: Output: πθ, ψλϕ .

heading: 5 MIRROR DESCENT ADVERSARIAL INVERSE REINFORCEMENT LEARNING
text: In this section, parameters θ, ϕ, and ν are presented representing agent policy, reward, and expert policy functions respectively parameterized with neural networks. On top of the standard AIL algorithm, we propose the MD-AIRL method, learning with a dual discriminator architecture motivated by previous studies regarding multiple discriminators (Nguyen et al., 2017; Chongxuan et al., 2017).
AIL methods have dealt divergences between joint densities of states and actions (Ghasemipour et al., 2020). To clearly distinguish matching overall state densities and imitating specific behavior, we aim to disentangle these concepts and propose two structured discriminators:
Dθ,ν(s, a) = σ ( log { πν(a|s)/πθ(a|s) } + b(s) ) and Db(s) = σ(b(s)),
where b : S → R is a neural network regarding states and σ(·) denotes the sigmoid function. The discriminators Db and Dθ,ν are trained with binary cross-entropy losses regarding trajectories:
maximize ν EπE [logDθ,ν(s, a)] + Eπt [log(1−Dθ,ν(s, a))], (12)
maximize b EπE [logDb(s)] + Eπt [log(1−Db(s))], (13)
Note that b and θ are not trained for learning Dθ,ν . Let ρπ ∈ ∆S denote normalized state visitation distribution of π, which is defined as ρπ(s) := (1 − γ)Eπ [∑∞ i=0 γ iI{si = s} ]
where I{·} is an indicator function. Since x 7→ α log σ(x) + β log(1 − σ(x)) attains its maximum at σ(x) = αα+β (Goodfellow et al., 2014), optimality is achieved when πν = πE and b(s) = log ρπE (s)/ρπ(s).
Let ψϕ denote the regularized reward function where the parameter ϕ is trained using a step size ηt, the agent policy πθ, and the estimation of expert policy πν :
minimize ϕ
Es∼τ̄t [ ηtDΩ(πϕ(·|s)∥πν(·|s)) + (1− ηt)DΩ(πϕ(·|s)∥πθ(·|s)) ] , (14)
where πϕ denotes the transformed policy ∇Ω(ψϕ) and τ̄t denotes mini-batches of states using the both agent and expert trajectories. The algorithm adjusts the term ηt = 1/αt with a harmonic progression by linearly increasing αt by the range of [α1, αT ], which is derived from our analyses.
For a hyperparameter λ > 0, the MD-AIRL reward function is defined by a linear combination:
ψλϕ(s, a) = λψϕ(s, a) + b(s).
By using arbitrary regularized RL which uses λΩ(·) as the regularization function, the reward learning of ψλϕ regarding the agent πθ is decomposed into the following two terms:
Eπθ [ ψλϕ(s, a)− λΩ(πθ(·|s)) ] = λEπθ [ ψϕ(s, a)− Ω(πθ(·|s)) ] −DKL(ρπθ∥ρπE )
= −λEπθ [ DΩ(πθ(·|s)∥πϕ(·|s)) ] −DKL(ρπθ∥ρπE ),
Minimizing the first term of Eπθ [ DΩ(π s θ∥πsϕ) ] represents the online learning of MD formulation. However, ψϕ that is trained on τ̄t cannot cover the entire reachable states since the state visitation is heavily misaligned in challenging sequential decision problems. Therefore, we propose the second auxiliary term of generalized KL divergence DKL(ρπθ∥ρπE ), playing an additional role of facilitating the supports of state visitation densities to be matched properly. Algorithm 1 summarizes the entire procedure of MD-AIRL. See Appendix C for detailed implementation.

heading: 6 EXPERIMENTAL RESULTS
text: For the RL algorithm, we implemented a RAC (Yang et al., 2019) method that is a generalization of SAC (Haarnoja et al., 2018) in terms of the regularization choice. We considered the class of separable regularizers Ω(p) = −Ea∼p[φ(p(a))]: (1) Shannon entropy (φ(x) = log(x)), (2) Tsallis entropy Tq (φ(x; q) = 1q−1 (x
q−1 − 1), q = 2 by default), (3) exp regularizer (φ(x) = e− ex), (4) cos regularizer (φ(x) = cos(π2x)), and (5) sin regularizer (φ(x) = 1 − sin π 2x). We evaluated our approach on three topics (bandits, multiple goals, and MuJoCo environments). The main comparative method was RAIRL with a density-based model (RAIRL-DBM) since this model shares the identical level of expressiveness as our method when the parameterization of ∆SA is specified.

heading: 6.1 LARGE-SCALE MULTI-ARMED BANDITS
text: We first considered multi-armed bandit problems, where the cardinality of action spaces is varied. Learning the optimal distribution of πE becomes challenging as |A| increases, because the frequency of each sample becomes sparse due to the curse of dimensionality (Bellman et al., 1957). For each experiment, a stateless expert distribution πE was generated by the parameters of softmax distribution πE(i) = exp(zi)/ ∑ j exp(zj) where the logits zi were randomly initialized to a uniform distribution. We set the action size to |A| = 102, 103, 104 and restricted the sample size to 16. Figure 5 shows that the Bregman divergence was large for MD-AIRL at the early training phase, because we chose the initial step size η1 to be greater than 1 (α1 = 0.5). MD-AIRL exceeded the discriminative performance of RAIRL after certain steps, while the progression of RAIRL mostly stopped at local minima. Table 1 shows that MD-AIRL achieved overall lower Bregman divergence on average when three different cardinalities and five regularizers were considered. MD-AIRL outperformed RAIRL in four cases by choosing effectively low step size at the ηT to be less than 1 (αT = 2). These results match properties of MD algorithms and our convergence analyses. Therefore, we argue that a constrained update rule with appropriate step sizes is necessary for robust reward acquisition and imitation for the situations when the total number of data samples is limited.

heading: 6.2 CONTINUOUS MULTI-GOAL ENVIRONMENT
text: We then considered an environment with a two-dimensional continuous state space. In this envioronment, an agent is a point mass initialized at the origin, and the four goals are located in the four cardinal directions. To draw meaningful reward surface, we considered multivariate Gaussian dis-
tribution policies with full covariance matrices. We parameterized a covariance matrix of Gaussian policy using lower triangular matrix which is an outcome of LDL decomposition (see Appendix C).
Figure 6 (a) shows trajectories generated by the trained agent. Figure 6 (b) shows that MD-AIRL achieved higher entropy for reaching the goals. Figure 6 (c) shows reward surfaces with regularizers, which was calculated by ψϕ(s, a) + φ(πθ(a|s)) for each point of a ∈ A and s = (5,−1). During the training, the MD reward was similar to the estimated ground truth using adversarial training. However, the surface of MD-AIRL became flatter than the ground-truth estimation when πt was sufficiently close to the expert behavior. As a result, we claim that drastic changes in the target distribution, which are one of the typical characteristics of adversarial frameworks, are prevented. We argue that these characteristics mitigate overfitting caused by unreliable discriminatve signals.

heading: 6.3 CONTINUOUS CONTROL: MUJOCO
text: We validated MD-AIRL on MuJoCo continuous control tasks (Brockman et al., 2016). We assumed diagonal Gaussian policies for both learner’s policy π and expert policy πE . Instead of the tanh-squashed policy (Haarnoja et al., 2018), we used the hyperbolized environment assumption of RAIRL, which means that tanh is regarded as a part of the environment.
For each tasks, we considered two different numbers of episodes collected by an expert policy. In Figure 7, the performance of MD-AIRL, RAIRL, and behavior cloning (bc; Pomerleau, 1991) al-
gorithms are shown with the expert and random agent performance. In terms of sample efficiency, MD-AIRL outperformed RAIRL on Hopper-v3, Walker-v3, and Ant-v3. The performance gaps between MD-AIRL and RAIRL were more prominent in 4 episodes of expert demonstrations. The training curves in Figure 8 indicate that MD-AIRL showed lower variance than RAIRL especially after early phase of training. It can be concluded that MD-AIRL inherits the sample efficiency of AIL algorithms in challenging RL tasks benchmarks. Additionally, the algorithm is highly stable with respect to limited sample sizes, which is in alignment with our theoretical analyses.

heading: 7 RELATED WORKS
text: Statistical manifolds. The Hessian of strongly convex function forms a metric tenser of a manifold called Hessian geometries (Shima, 2007). Bregman divergences are similar to these metrics, providing useful metric-like properties (Butnariu & Resmerita, 2006). Using the Bregman divergence allows to solve many optimization problems, generalizing traditional approaches such as least squres (Boyd et al., 2004; Hiriart-Urruty & Lemaréchal, 2004). Probability distributions can be considered as points in a geometric space. A representative statistical manifold is the information geometry induced by the Fisher-Rao metric (Amari, 2016; Bauer et al., 2016; Nielsen, 2020).
Regularized IRL. Energy-based policies (i.e. Boltzmann distributions) have been appeared in early IRL researches (Ramachandran & Amir, 2007; Neu & Szepesvári, 2007; Babeş-Vroman et al., 2011). Notably, MaxEnt IRL (Ziebart et al., 2008; 2010) is a representative IRL algorithm based on information theory. Other statistical entropies have also been applied to the imitation learning problem, such as the Tsallis entropy derived from Tsallis statistics (Lee et al., 2018). Compared to RAIRL which also can use various convex regularizers (Jeon et al., 2020), our work allows more NEW R1, R2 realistic situations where the expert policy cannot be precisely attained due to insufficient data, thanks to theoretical foundations originated from optimization studies.
Mirror descent. MD is closely related to algorithms regarding non-Euclidean geometries with discretization of steps such as natural gradients (Amari, 1998; Raskutti & Mukherjee, 2015; Gunasekar et al., 2020). On the primal space, the infinitesimal limit of MD step corresponds to a Riemannian gradient flow (Do Carmo, 2016; Gunasekar et al., 2020). The online MD algorithms possess rigorous regret bounds (Srebro et al., 2011; Lei & Zhou, 2020); thus they can be highly efficient in terms of the number of evaluations until convergence.

heading: 8 DISCUSSION AND CONCLUSIONS
text: In this paper, we presented a novel IRL framework. We provided an mirror descent solution on reward functions and corresponding theoretical arguments. We proposed MD-AIRL, a practical FIX R3 adversarial IRL framework that can solve challenging imitation learning tasks. We verified that the proposed method has clear advantages over previous AIL methods in terms of robustness. As the reward hypothesis is grounded in obtainining a robust representation of the expert policy using IRL, we argued that current IRL studies lack robustness and theoretical guarantees for practical situations.
Considering RL and its inverse problem with geometric perspectives is vital for achieving desired goals in realistic situations. Although our work covers various online imitation learning methods with MD, it does not include some other cases when the proximity term is of other statistical divergence families such as f-divergence (Amari, 2016). Additionally, the boundedness and continuity assumptions on Bregman divergences in our analyses are usually justified, but outliers exist such as NEW R3 KL divergences. More sophisticated analyses on these general cases remain as future works.

heading: A PROOFS
text: We define the space of policy as an object ∆SA, which is a vector space formed by a collection of |S| elements of unit (|A| − 1)-simplexes: ∆A = { x1e1 + · · · + x|A|e|A| | ∑|A| i=1 θi = 1 and θi ≥
0 for i ∈ A }
. We assume that each space of simplex is a subset of a specific Banach space called Lp space (RA, ∥·∥), where ∥·∥ is a p-norm on A. The dual space of Lp space for 1 < p <∞ is Lq space (RA, ∥·∥∗), where ∥·∥∗ is defined as a q-norm (1/p + 1/q = 1). We assume 1 < p ≤ 2 for the convergence property in the dual Lq space. We start with following definitions.
Definition 2 (Lipschitz constants). Given two metric spaces (X, dX) and (Y, dY ) where dX denotes the metric on set X and dY is the metric on set Y , a function f : X → Y is called Lipschitz continuous if there exists a real constant k ≥ 0 such that, for all x1 and x2 in X ,
dY (f(x1), f(x2)) ≤ k · dX(x1, x2). (15)

heading: In particular, a function f is called Lipschitz continuous if there exists a constant k ≥ 0 such that,
text: ∥f(x1)− f(x2)∥∗ ≤ k · ∥x1 − x2∥, ∀x1, x2 (16)
where norms ∥·∥ and ∥·∥∗ are endowed with each space X and Y respectively. For the smallest L that substitutes k, L is called the Lipschitz constant and f is called the L-Lipschitz continuous function.
Definition 3 (Discrete-time martingale). If a stochastic process {Zt}t≥1 satisfies E[|Zn|] ≤ ∞ and
1) E[Zn+1|X1, . . . , Xn] ≤ Zn, 2) E[Zn+1|X1, . . . , Xn] = Zn, 3) E[Zn+1|X1, . . . , Xn] ≥ Zn,
then, the stochastic process {Zt}t≥1 is called a 1) submartingale, 2) martingale, and 3) supermartingale, respectively, with respect to a filtration {Xt}t≥1.
The following analyses and proofs follow the results appeared in previous literatures for general aspects (Nemirovsky & Yudin, 1983; Gunasekar et al., 2020; Srebro et al., 2011; Lei & Zhou, 2020; Beck & Teboulle, 2003; Raskutti & Mukherjee, 2015). Our analyses extend existing results to imitation learning, and they are also highly general to cover various online methods for sequential decision problems.

heading: A.1 PROOF OF LEMMA 1
text: Proof of Lemma 1. The conjugate operator of ψsπ satisfies the following identity (Lemma 1 of Jeon et al. (2020))
Ω∗(ψsπ) = max π̃s∈RA ⟨π̃s, ψsπ⟩A − Ω(π̃s)
= max π̃s∈RA
⟨π̃s,∇Ω(πs)⟩A − 〈 πs,∇Ω(πs) 〉 A +Ω(π s)− Ω(π̃s)
= min π̃s∈RA
Ω(π̃s)− Ω(πs)− ⟨∇Ω(πs), π̃s − πs⟩A
= min π̃s∈RA
DΩ(π̃ s∥πs)
for every state s ∈ S. By the property of Bregman divergence, and the convexity of DΩ(π̃s∥πs) with respect to π̃s (Geist et al., 2019; Acharyya et al., 2013), the optimal condition is obtained by the unique maximizing argument π̃(·|s) = π(·|s). By taking gradient to both sides with respect to ψsπ we yield π = ∇Ω∗(ψπ). If there is another π̃ ̸= π that makes ψπ̃ = ψπ , this contradicts the property of unique maximizing arguments for conjugates. Therefore, ψπ is uniquely defined for each π and ∇Ω∗(ψ) ∈ ∆SA.
We model the reward approximator regarding actions as ψ ∈ ΨΩ(∆SA) = {ψ | ψ(s, a) = ψπ(s, a),∀s ∈ S, a ∈ A, π ∈ ∆SA}, where ψπ is defined ψπ(s, a) := Ω′(s, a;π) − Ea′∼π(·|s) [ Ω′(s, a′;π) ] + Ω(π(·|s)). Note that ΨΩ(π) is different from ∇Ω(π), but it is shifted
by the amount of −Ea′∼π(·|s) [ Ω′(s, a′;π) ] + Ω(π) for every element of ψπ(s, ·). Therefore, the operators ΨΩ and ∇Ω are link functions that form natural isomorphisms of functions.

heading: A.2 PROOF OF THEOREM 1
text: We consider the unique fixed point of π∗ as the solution of infπ∈∆SA limt→∞ E[ft(π, τt)] where the expectation indicates that we take consideration for all τt and ft. By equating derivatives to zero, we write ∇Ω(π∗) = limt→∞ E[∇Ω(π̄E,t)]). We provide more general results than previous imitation learning, as some of following statements include the case of infπ∈∆SA E[ft(π, τ)] > 0, which means the estimates {π̄E,t}∞t=1 do not converge to the fixed point of π∗, so limt→∞∥π∗− π̄E,t∥ > 0. Thus, MD-AIRL and other online MD algorithms allow realistic settings such as scarcity of data or imperfect demonstrations.
We first introduce a key relationship regarding cumulative gradients in our online MD setting. Lemma 2. Let {πt}∞t=1, {π̄E,t}∞t=1, and {ηt}∞t=1 be policy, estimate, and step size sequences, respectively. The subsequent policy πt+1 in Equation (5) is obtained by an RL algorithm using the derivation of ψt+1 in Equation (7), resulting to the following equation:
πt+1(·|s) = argmin πs∈∆A ηtDΩ(π s∥π̄sE,t) + (1− ηt)DΩ(πs∥πst ), ∀s ∈ S. (17)
We have for t ∈ N, ηt ( ∇Ω(πt)−∇Ω(π̄E,t) ) = ∇Ω(πt)−∇Ω(πt+1) (18)
where we write∇Ω(π) = [ ∇Ω(π(·|s)) ] s∈S
Proof of Lemma 2. Since the optimization is convex with respect to each πs, we equate the derivatives to 0 at πt+1 as
ηt(∇Ω(πst+1)−∇Ω(π̄sE,t)) + (1− ηt) ( ∇Ω(πst+1)−∇Ω(πst ) ) = 0, ∀s ∈ S.
Finally, we derive Equation (18) as ηt(∇Ω(πst+1)−∇Ω(π̄sE,t)) + (1− ηt) ( ∇Ω(πst+1)−∇Ω(πst ) ) = 0
⇐⇒ ∇Ω(πst+1)− ηt∇Ω(π̄sE,t)− (1− ηt)∇Ω(πst ) = 0 ⇐⇒ ∇Ω(πst )−∇Ω(πst+1) = ηt ( ∇Ω(πst )−∇Ω(π̄sE,t) ) .
Therefore, the proof is complete.
Lemma 2 indicates that the distances between dual maps are equivalent to ηt∥∇Ω(π̄sE,t)−∇Ω(πst )∥∗ Therefore, if the step size decreases as limt→∞ ηt = 0, limt→∞∥∇Ω(πt)−∇Ω(πt+1)∥∗ = 0; thus, our argument in Section 3 is reasonable when Ω is strongly smooth.
Next, we reintroduce the three-point identity as follows. Lemma 3 (Three-point identity). Let πa,πb, and πc be any policies with a given state. We have the following identity:〈
∇Ω(πa)−∇Ω(πb), πc − πb 〉 A = DΩ(πc∥πb)−DΩ(πc∥πa) +DΩ(πb∥πa).
Proof of Lemma 3. This can be derived using the definition of divergence as follows. DΩ(πc∥πb)−DΩ(πc∥πa) +DΩ(πb∥πa) = Ω(πc)− Ω(πb)− 〈 ∇Ω(πb), πc − πb 〉 A
− Ω(πc) + Ω(πa) + 〈 ∇Ω(πa), πc − πa 〉 A
+Ω(πb)− Ω(πa)− 〈 ∇Ω(πa), πb − πa 〉 A
= 〈 ∇Ω(πa)−∇Ω(πb), πc − πb 〉 A.
Therefore, the proof is complete.
We now introduce basic identities regarding a Bregman divergence Lemmas 4 and 5 that are used to address progress of the algorithm. Lemma 4. Let πa,πb, and πc be any policies with a given state. The following identity holds.
DΩ(πc∥πb)−DΩ(πc∥πa) = DΩ(πa∥πb) + 〈 ∇Ω(πa)−∇Ω(πb), πc − πa 〉 A. (19)
Proof of Lemma 4. By Lemma 3, we have
DΩ(πc∥πb)−DΩ(πc∥πa) = −DΩ(πb∥πa) + ⟨∇Ω(πa)−∇Ω(πb), πc − πb⟩A. Utilizing the identity of two Bregman divergences
DΩ(π, π̃) +DΩ(π̃, π) = 〈 ∇Ω(π)−∇Ω(π̃), π − π̃ 〉 A, (20)
we separate πc − πb into πc − πa and πa − πb and write the rest of derivation as followings. DΩ(πc∥πb)−DΩ(πc∥πa)
= −DΩ(πb∥πa) + ⟨∇Ω(πa)−∇Ω(πb), πa − πb⟩A︸ ︷︷ ︸ Eq. (20) +⟨∇Ω(πa)−∇Ω(πb), πc − πa⟩A
= DΩ(πa∥πb) + 〈 ∇Ω(πa)−∇Ω(πb), πc − πa 〉 A
Therefore, we obtain the desired identity.
Lemma 5. Let πa,πb, and πc be any policies with a given state. The following identity holds. DΩ(πb∥πa)−DΩ(πc∥πa) = − 〈 ∇Ω(πc)−∇Ω(πa), πc − πb 〉 A +DΩ(πb∥πc). (21)
Proof of Lemma 5. By Lemma 3, we have
DΩ(πb∥πa)−DΩ(πc∥πa) = −DΩ(πc∥πb) + ⟨∇Ω(πa)−∇Ω(πb), πc − πb⟩A. We separate∇Ω(πa)−∇Ω(πb) into∇Ω(πa)−∇Ω(πc) and∇Ω(πc)−∇Ω(πb) and write the rest of derivation as followings.
DΩ(πb∥πa)−DΩ(πc∥πa) = −DΩ(πc∥πb) + ⟨∇Ω(πc)−∇Ω(πb), πc − πb⟩A︸ ︷︷ ︸
Eq. (20)
+⟨∇Ω(πa)−∇Ω(πb), πc − πb⟩A
= DΩ(πb∥πc) + 〈 ∇Ω(πa)−∇Ω(πc), πc − πb 〉 A
Therefore, we obtain the desired identity.
We show the key arguments to prove Theorem 1 in the following lemma. Lemma 6. Assume infπ∈∆SA E[ft(π, τ)] > 0. Assume also that Ω is ω-strongly convex and ∇Ω is L-Lipschitz continuous for ω ≥ 0 and L ≥ 0. If limt→∞ Eτ1,···,τt [∑∞ i=0 γ iDΩ ( πsit ∥∥πsiE )] = 0 for πE , then {ηt}∞t=1 satisfies Equation (9). Furthermore, if Ω is strongly smooth, then Theorem 1 (a) holds with some constants n ∈ N and c > 0.
Proof of Lemma 6. First, we show the condition of limt→∞ ηt = 0. Assuming all states are reachable (Fu et al. 2017, Definition B.1), the condition limt→∞ Eτ1,···,τt [∑∞ i=0 γ iDΩ ( πsit ∥∥πsiE )] = 0, NEW R3
implies limt→∞ Eτ1,···,τt [ ∥πt − πE∥ ] = 0, where ∥·∥ is a matrix norm induced by the p-norm on A. Then, our aim is to show that the gradient of strong convex function for πt converges to∇Ω(πE) as
lim t→∞ Eτ1,···,τt [∥∥∇Ω(πt)−∇Ω(πE)∥∥∗] = 0. (22)
To prove this point, we use the continuity of ∇Ω at πE; for any ε > 0, there exists some 0 < δ ≤ 1 such that ∥∇Ω(π)−∇Ω(πE)∥∗ < ε whenever ∥π − πE∥ < δ. When ∥π − πE∥ ≥ δ, we apply the L-Lipschitz continuity assumption to find∥∥∇Ω(π)−∇Ω(πE)∥∥∗ ≤ L∥π − πE∥, (23) where ∥·∥∗ is a matrix norm induced by the q-norm. Combining Equation (22) and Equation (23), we know that
Eτ1,···,τt [ ∥∇Ω(πt)−∇Ω(πE)∥∗ ] ≤ ε+ L · Eτ1,···,τt [ ∥πt − πE∥ ] . (24)
However, limt→∞ Eτ1,···,τt [ ∥πE − πt∥ ] = 0 ensures the existence of some n ∈ N such that for
t > n, there holds Eτ1,···,τt [ ∥πE − πt∥ ] < ε/L. Applying this inequality to Equation (24), we have
Eτ1,···,τt [ ∥∇Ω(πt)−∇Ω(πE)∥∗ ] < 2ε for some t > n.
For temporal estimation, let us denote the infimum of the expectation ℓ = infπ∈∆SA E[∇Ω(πt) − ∇Ω(π̄E,t)] > 0. From Lemma 2, we have ηt(∇Ω(πt)−∇Ω(π̄E,t)) = ∇Ω(πt)−∇Ω(πt+1). Taking the expectations on both sides of this equation yields
ηtℓ ≤ ηtEτ1,···,τt+1 [ ∥∇Ω(πt+1)−∇Ω(π̄E,t)∥∗ ] = Eτ1,···,τt+1 [ ∥∇Ω(πt)−∇Ω(πt+1)∥∗ ] .
Hence the convergence of the point∇Ω(πt) is confirmed by taking the limit: limt→∞ ηt = 0. Next, we show ∑∞ t=1 ηt = ∞. By the ω-strong convexity by the L-Lipschitz continuity of Ω, we can find inequalities as〈 ∇Ω(π)−∇Ω(π̃), π − π̃ 〉 A ≤ L∥π − π̃∥ 2 ≤ 2L ω DΩ(π∥π̃). (25)
We note that ∥πt+1 − π̄E,t∥ ≤ ∥πt − π̄E,t∥ so that there is a constant ε that satisfies E[∥πt+1 − π̄E,t+1∥] ≥ E[∥πt+1− π̄E,t∥]+ε. Therefore, taking expectations in Equation (21) (setting πa = π̄E,t, πb = πt+1, and πc = πt), for limt→∞ ηt = 0 and the strongly convex Ω, we can find
Eτ1,···,τt+1 [ DΩ(πt+1∥π̄E,t+1) ] ≥ Eτ1,···,τt+1 [ DΩ(πt+1∥π̄E,t) ] + ε′
≥ (1− aηt)Eτ1,···,τt [ DΩ(πt∥π̄E,t) ] + Eτ1,···,τt+1 [ DΩ(πt+1∥πt) ] + ε′ (Eq. (21))
≥ (1− aηt)Eτ1,···,τt [ DΩ(πt∥π̄E,t) ] + ε′′ (26)
for some t and 0 < ε′ < ε′′. The positive constant a = 2L/ω is derived by the inequalities in Equation (25). We omit the conditioned states (superscripts) for simplicity of the derivation.
Since limt→∞ ηt = 0, we can also find a constant n ∈ N such that ηt ≤ (3a)−1 for t ≥ n. Applying the inequality 1− x > exp(−2x) for x ∈ (0, 1/3], we derive another inequality
Eτ1,···,τt+1 [ DΩ ( πt+1 ∥∥π̄E,t+1)] ≥ exp(−2aηt)Eτ1,···,τt[DΩ(πt∥∥π̄E,t)], ∀t ≥ n (27) Applying this for t = T − 1, . . . , n yields
Eτ1,···,τT [ DΩ ( πT ∥∥π̄E,T)] ≥ [ T∏
t=n+1
exp ( −2aηt )] Eτ1,···,τn [ DΩ ( πn ∥∥π̄E,n)]
= exp ( −2a ·
T∑ t=n+1 ηt ) Eτ1,···,τn [ DΩ ( πn ∥∥π̄E,n)].
(28)
Using Equation (28), we conclude Eτ1,···,τn [ DΩ(πn+1∥π̄E,n+1) ] > 0. Otherwise, we have
Eτ1,···,τn [ DΩ(πn∥π̄E,n) ] = Eτ1,···,τn+1 [ DΩ(πn+1∥π̄E,n+1) ] = 0
according to Equation (27), which leads to Eτ1,···,τn [ ∥πn − π̄E,n∥2 ] = Eτ1,···,τn+1 [ ∥πn+1 −
π̄E,n+1∥2 ]
= 0. This implies πn = π̄E,n = πn+1almost surely, leading to E[ft(πt, τt)] = 0. This is a contradiction to the previous assumption infπ∈∆SA E[ft(π, τt)] > 0, thus Eτ1,···,τn+1 [ DΩ(πn+1∥π̄E,n+1) ] > 0. Let us suppose the ideal case that the estimation process
learns the exact πE in t→∞. In order to satisfy the limit limT→∞ Eτ1,···,τT [ DΩ(πT∥π̄E,T ) ] = 0 we
see from Equation (28) that ∑∞
t=1 ηt =∞.
Now, we show that Theorem 1 (a) holds. Since Ω is ω-strongly convex, so Ω∗ is (ω−1)-strongly smooth with respect to ∥·∥∗. Additionally, the L-Lipschitz continuity of ∇Ω implies L-strong smoothness of Ω; thus, naturally, Ω∗ is L-strongly convex.
Combining these, for ∀t ≥ n, the condition ηt ≤ (3a)−1 induces Eτ1,···,τt+1 [ DΩ(πt+1∥π̄E,t+1) ] ≥ (1− aηt)Eτ1,···,τt [ D(πt∥π̄E,t) ] + (2L)−1Eτ1,···,τt+1 [ ∥∇Ω(πt)−∇Ω(πt+1)∥2∗ ] ,
and by Lemma 2, we get Eτ1,···,τt+1 [ DΩ(πt+1∥π̄E,t+1) ] ≥ (1− aηt)Eτ1,···,τt [ DΩ(πt∥π̄E,t) ] + (2L)−1η2tEτ1,···,τt+1 [ ∥∇Ω(πt)−∇Ω(π̄E,t)∥2∗ ] .
Using the Cauchy-Schwarz inequality, we obtain a lower bound of the last term as Eτ1,···,τt [ ∥∇Ω(πt)−∇Ω(π̄E,t)∥2∗ ] ≥ { Eτ1,···,τt [ ∥∇Ω(πt)−∇Ω(π̄E,t)∥∗ ]}2 ≥ ℓ2. Thus, we obtain the final inequality as
Eτ1,···,τt+1 [ DΩ(πt+1∥π̄E,t+1) ] ≥ (1− aηt)Eτ1,···,τt [ DΩ(πt∥π̄E,t) ] + (2L)−1(ηtℓ) 2, ∀t ≥ n.
Applying this inequality from t = T ≥ n+ 1 to t = n+ 1, we obtain
Eτ1,···,τT+1 [ DΩ(πT+1∥π̄E,T+1) ] ≥ Eτ1,···,τn [ DΩ(πn∥π̄E,n) ] T∏ t=n+1 (1− aηt)
+ (2L)−1ℓ2 T∑
t=n+1
η2t T∏ k=t+1 (1− aηk)
≥ (2L)−1ℓ2 T∑
t=n+1
η2t T∏ k=t+1 (1− aηk).
By the Cauchy-Schwarz inequality and our bound 0 < 1− aηk ≤ 1 for k ≥ n, we have T∑
t=n+1
ηt T∏ k=t+1 (1− aηk) ≤
{ T∑
t=n+1
η2t T∏ k=t+1 (1− aηk)
}1/2 (T − n)1/2.
Hence T∑
t=n+1
η2t T∏ k=t+1 (1− aηk) ≥ 1 a2(T − n)
( T∑
t=n+1
aηt T∏ k=t+1 (1− aηk)
)2
= 1
a2(T − n)
( T∑
t=n+1
(1− (1− aηt)) T∏
k=t+1
(1− aηk)
)2
= 1
a2(T − n)
( T∑
t=n+1
[ T∏
k=t+1
(1− aηk)− T∏
k=t
(1− aηk)
])2
≥ 1 a2(T − n)
( T∑
t=n+1
1− T∏
k=t
(1− aηk)
)2
≥ 1 a2(T − n) (1− (1− aηn+1))2 = η2n+1 T − n
Therefore, we obtain the lower bound of
Eτ1,···,τT [ DΩ(πT+1∥π̄E,T+1) ] ≥ η2n+1(2L) −1ℓ2
T − n .
Since the Bregman divergence is bounded for all states, the sequence {γiDΩ(πsit ∥π si t )} will converge as i→∞. Applying the monotone convergence theorem, we can interchange expectation and summation, which yields
Eτ1,···,τT [ ∞∑ i=0 γiDΩ ( πT (·|si) ∥∥π̄E,T (·|si))] = ∞∑ i=0 Eτ1,···,τT [ γiDΩ(πT (·|si)∥π̄E,T (·|si)) ] =
∞∑ i=0 γiEτ1,···,τT [ DΩ(πT (·|si)∥π̄E,T (·|si)) ] ≥ η2n+1(2L− 2Lγ)−1ℓ2
T − n , ∀T ≥ n.
This verifies Theorem 1 (a) with the constant c = η2n+1(2L− 2Lγ)−1ℓ2.
Last, we show convergence to a unique fixed point of π∗ using the particular form of ηt in Equation (9). Lemma 7. If {ηt}∞t=1 satisfies Eq. (9), limt→∞ Eτ1,···,τt [∑∞ i=0 γ iDΩ(π s ∗∥πst ) ] = 0. Furthermore,
if the step size takes the form ηt = 4t+1 , then Eτ1,···,τT [∑∞ i=0 γ iDΩ ( πsi∗ ∥∥πsiT )] = O(1/T ).
Proof of Lemma 7. According to Lemma 4 and the Bregman divergence of the conjugate DΩ∗ , the one step progress regarding π̄E,t can be written as DΩ(π∗∥πt+1)−DΩ(π∗∥πt) = 〈 ∇Ω(πt)− Ω(πt+1), π∗ − πt 〉 A +DΩ(πt∥πt+1) = ηt⟨∇Ω(πt)−∇Ω(π̄E,t), π∗ − πt⟩+DΩ∗(∇Ω(πt+1)∥∇Ω(πt)), where we omit a given state. As ω-strong convexity of Ω implies the (ω−1)-strong smoothness of Ω∗, we have
DΩ∗(∇Ω(πt)∥∇Ω(πt+1)) ≤ 1
2ω ∥∇Ω(πt)−∇Ω(πt+1)∥2∗ = η2t 2ω ∥∇Ω(πt)−∇Ω(π̄E,t)∥2∗ (29)
We bound ∥∇Ω(πt)−∇Ω(π̄E,t)∥2∗ by 2∥∇Ω(πt)−∇Ω(π∗)∥2∗ +2∥∇Ω(π∗)−∇Ω(π̄E,t)∥2∗ , following the work of Lei & Zhou (2020). Since ∇Ω is cocoercive with 1L bye the Lipschitz continuity, we obtain ∥∇Ω(πt)−∇Ω(π∗)∥2∗ ≤ L⟨∇Ω(π∗)−∇Ω(πt), π∗ − πt⟩ thus
DΩ(π∗∥πt+1)−DΩ(π∗∥πt) ≤ ηt⟨∇Ω(π∗)−∇Ω(π̄E,t), π∗ − πt⟩
− (1− ηtL ω )ηt⟨∇Ω(π∗)−∇Ω(πt), π∗ − πt⟩+ η2t ω
[ ∥∇Ω(π∗)−∇Ω(π̄E,t)∥2∗ .
(30)
By taking expectation, it follows that there exists n ∈ N such that ηt ≤ ω2L for t ≥ n holds Eτt [ DΩ(π∗∥πt+1) ] ≤ DΩ ( π∗∥πt ) − ηt
2 DΩ(π∗∥πt) + zη2t , (31)
where z is the constant z = 1ωE[∥∇Ω(π∗) − Ω(π̄E,t)∥ 2 ∗ ]. Let {At}∞t=1 denote a sequence of At = Eτ1,···,τt [ DΩ(π∗∥πt) ] . Then we have
At+1 ≤ ( 1− ηt
2
) At + zη 2 t , ∀t ≥ n. (32)
For a constant h > 0, we claim that At1 < h for some t1 > n ′. Assume that this is not true, and we find some t2 ≥ t1 such that At > h, ∀t ≥ t2. Since limt→∞ ηt = 0, there are some t > t3 > t2 that ηt ≤ h4b . However, Equation (32) tells us that for t ≥ t3,
At+1 ≤ ( 1− ηt
2
) At + zη 2 t ≤ At3 − h
4 t∑ k=t′γ ηk → −∞ (as t→∞).
This is a contradiction, which verifies At < h for t > n′. Since limt→∞ ηt = 0, we can find some ηt that makes At to be monotone decreasing. Then, we can conclude that the nonnegative sequence {At}∞t=1 converges by iteratively applying the upper bounds. We now prove Theorem 1 (b) under the condition and the choice ηt = 4t+1 of the step size sequence. The estimate becomes
At+1 ≤ ( 1− 2
t+ 1
) At + 16z
(t+ 1)2 , ∀t ≥ n.
It follows the recurrence relation as t(t+ 1)At+1 ≤ (t− 1)tAt + 16z, ∀t ≥ n.
Applying this relation iteratively, we obtain the general form of inequality. (T − 1)TAT ≤ (n− 1)nAn + 16z(T − n), ∀T ≥ n,
therefore we get the inequality as follows: Eτ1,···,τT [ DΩ(π∗∥πT ) ] ≤
(n− 1)nEτ1,···,τn [ DΩ(πE∥πn) ] (T − 1)T + 16z T , T ≥ n.
By applying the monotone convergence theorem similar to Lemma 6, we conclude that Eτ1,···,τT [ ∑∞ i γ iDΩ(π∗(·|si)∥πT (·|si))] = O(1/T ).

heading: A.3 PROOF OF THEOREM 2
text: Necessity. We rewrite the inequality in Equation (26) as Eτ1,···,τt+1 [ DΩ ( πt+1 ∥∥π̄E,t+1)] ≥ (1− 2Lω−1ηt)Eτ1,···,τt[DΩ(πt∥∥π̄E,t)]. (33) Since we assume that ηt converge to 0 from previous arguments, onsider the step size sequence FIX R3 0 < ηt ≤ ω(2+κ)L for κ > 0 and t ≥ n where n ∈ N. Denote a constant ã = 2+κ 2 log 2+κ κ and apply the elementary inequality
1− x ≥ exp(−ãx), ∀0 < x ≤ 2 2 + κ
From Equation (33), it can be obtained that Eτ1,···,τt+1 [ DΩ ( πt+1 ∥∥π̄E,t+1)] ≥ exp(−2ãLω−1ηt)Eτ1,···,τt[DΩ(πt∥π̄E,t)]. Applying this inequality iteratively for t = n, . . . , T − 1, then gives
Eτn,···,τT [ DΩ(πT∥π̄E,T ) ] ≥ T−1∏ t=n exp ( −2ãLω−1ηt ) DΩ(πn∥π̄E,n)
= exp { −2ãLω−1
T−1∑ t=n ηt
} DΩ(πn∥π̄E,n).
From the assumption πE ̸= πn, we have DΩ(πn∥π̄E,n) > 0. The convergence limt→∞ Eτ1,···,τt [ DΩ(πt∥π̄E,t) ] = 0 then implies ∑∞ t=1 ηt =∞.
Sufficiency. Here we use the estimate (32) derived in the proof of Lemma 7. However, in the optimal case, z = 1ωE[ft(π∗, τ)] = 0, so (32) takes the form (we can choose n = 1 by Equation (31)
At+1 ≤ ηt 2 At, ∀t ∈ N. (34)
This implies that for any 0 < h, there must exist some integer t1 ∈ N such that At ≤ γ for t ≥ t1. Otherwise, At > γ for every t ≥ t2 with t2 ≥ t1, which leads to a contradiction:
At+1 ≤ At2 − h
2 t∑ k=t̃γ ηk → −∞ (as t→∞).
Equation (34) also tells us that the sequence {At}∞t=1 of nonnegative numbers is monotone decreasing. Hence At ≤ h for every t ≥ t1, which proves the convergence of At
lim t→∞
Eτ1,···,τt [ DΩ(π∗∥πt) ] = lim
t→∞ At = 0.
We now prove the second point in Theorem 2 which is under the special choice of constant step size sequence ηt ≡ η1. It follows from Equation (33) that AT ≥ (1− 2Lω−1η1)T−1A1. Hence, Eq (34) translates to
At+1 ≤ (1− η1/2)At, from which we find AT ≤ (1 − η1/2)T−1A1 by iteration starting from t = 1. This verifies the theorem with c1 = (1− 2Lη1ω ) and c2 = (1− η1 2 ).

heading: A.4 PROOF OF PROPOSITION 1
text: The proof is based on Doob’s forward convergence theorem.
Theorem 3 (Doob’s forward convergence theorem). Let {Xt}t∈N be a sequence of nonnegative random variables and let {Ft}t∈N be a filtration with Ft ⊂ Ft+1 for every t ∈ N. Assume that E [ Xt+1|Ft ] ≤ Xt almost surely for every t ∈ N. Then the sequence {Xt} converges to a nonnegative random variable X∞ almost surely.
We follow the proof of Lemma 7 and apply Eq. (30). Since ⟨π∗ − πt,∇Ω(π∗) − ∇Ω(πt)⟩ ≥ 0, Equation (30) implies
Eτt [ DΩ(π∗∥πt+1) ] ≤ DΩ(π∗∥πt) + η2t ω Eτt [ ∥∇Ω(π∗)−∇Ω(π̄E,t)∥2∗ ] , ∀t ≥ n. (35)
The condition ∑∞
t=1 η 2 t <∞ enables us to define a stochastic process {Xt} by
Xt = DΩ(π∗∥πt+1) + 1 ω E [ ∥∇Ω(π∗)−∇Ω(π̄E,t)∥2∗ ] ∞∑ i=t+1 η2i .
By (35), we know that Eτt [Xt+1] ≤ Xt for t ≥ n. Additionally, Xt ≥ 0. Therefore, the stochastic process {Xt}t−n+1≥1 is a submartingale (equivalently, {−Xt}t−n+1≥1 is a supermartingale). Then by Lemma 3, the sequence {Xt}t≥1 converges to a nonnegative random variable X∞ almost surely. According to Fatou’s Lemma and the convergence limt→∞ E [∑∞ i=0 γ iDΩ(π∗∥πt) ] = 0 for any states a proved by Lemma 7, we obtain NEW
E[X] = E [ lim t→∞ ∞∑ i=0 γiDΩ ( π∗(·|si) ∥∥πt(·|si))] ≤ (1− γ)−1 lim inf t→∞ E [ DΩ(π∗∥πt) ] = 0.
It follows that the sequence of discounted sum {∑∞
i=0 γ iDΩ(π∗(·|si)∥πt(·|si)) } t∈N converges to 0
almost surely.

heading: B TSALLIS ENTROPY AND ASSOCIATED BREGMAN DIVERGENCE AMONG FULL COVARIANCE MULTIVARIATE GAUSSIAN DISTRIBUTIONS
text: In this section, we briefly reintroduce the tractable method to derive Bregman divergences and regularized reward functions proposed by Nielsen & Nock (2011) and Jeon et al. (2020). Then, we identify specific paremeterization to model Gaussian distributions with full covariance matrices.
The standard form of the exponential family is represented as exp {〈 θ, t(x) 〉 − F (θ) + k(x) } . (36)
The parameterization is as follows:
θ = [ Σ−1µ − 12Σ −1 ] = [ θ1 θ2 ] ,
t(x) = [ x xxT ] ,
F (θ) = −1 4 θT1 θ −1 2 θ1 + 1 2 ln|−πθ−12 | = 1 2 µTΣ−1µ+ 1 2 ln(2π)d|Σ|,
k(x) = 0,
where we can analytically recover the multivariate Gaussian distribution (Nielsen & Nock, 2011) exp {〈 θ, t(x) 〉 − F (θ) + k(x) } (37)
= exp { µTΣ−1x− 1 2 tr ( Σ−1xxT ) − 1 2 µTΣ−1µ+ 1 2 ln(2π)d|Σ| } (38)
= 1
(2π)d/2|Σ|1/2 exp
{ µTΣ−1x− 1
2 xTΣ−1x− 1 2 µTΣ−1µ
} (39)
= 1
(2π)d/2|Σ|1/2 exp
{ 1
2 (x− µ)TΣ−1(x− µ)
} . (40)
For two distributions π and π̂ with k(x) = 0, Nielsen & Nock (2011) proposed the function I(·):
I(π, π̂;α, β) = ∫ π(x)απ̂(x)β dx = exp { F ( αθ + βθ̂ ) − αF (θ)− βF (θ̂) }
where the detailed derivation is as follows:∫ π(x)απ̂(x)β dx
= ∫ exp { α 〈 θ, t(x) 〉 − αF (θ) + β 〈 θ̂, t(x) 〉 − βF (θ̂) } dx
= ∫ exp {〈 αθ + βθ̂, t(x) 〉 − F ( αθ + βθ̂ )} exp { F ( αθ + βθ̂ ) − αF (θ)− βF (θ̂) } dx = exp { F ( αθ + βθ̂ ) − αF (θ)− βF (θ̂) }∫ exp {〈 αθ + βθ̂, t(x) 〉 − F ( αθ + βθ̂ )} dx
= exp { F ( αθ + βθ̂ ) − αF (θ)− βF (θ̂) } .

heading: B.1 TSALLIS ENTROPY OF FULL COVARIANCE GAUSSIAN DISTRIBUTIONS
text: For φ(x; q) = 1q−1 (x q−1 − 1), the Tsallis entropy can be written as Tq(π) := −Ex∼πφ(x; q) = ∫ π(x) 1− π(x)q−1
q − 1 dx
= 1−
∫ π(x)q dx
q − 1 =
1 q − 1 ( 1− I(π, π; q, 0) ) = 1− exp ( F (qθ)− qF (θ) ) q − 1 .
If π is a multivariate Gaussian distribution, we have
F (qθ) = q
2 µTΣ−1µ+
1 2 ln(2π)d|Σ| − 1 2 ln qd.
Since the covariance matrix is positive-definite, Cholesky decomposition can be applied, which separates the matrix to lower- and upper-triangle matrices. Likewise, we can apply LDL decomposition. Let us factorize the covariance Σ = Ldiag{σ21 , . . . , σ2d}LT where L denotes a unit lower triangular matrix produced by LDL decomposition. Then we have
F (qθ)− qF (θ) = (1− q) { d
2 ln 2π +
1 2 ln|Σ| − d ln q 2(1− q) } = (1− q) { d
2 ln 2π +
1 2 ln d∏ i=1 σ2i − d ln q 2(1− q)
}
= (1− q) d∑
i=1
{ ln 2π
2 + lnσi −
ln q
2(1− q)
} .
B.2 TRACTABLE FORM OF ψπ
For separable Ω, ψπ is written as (Jeon et al., 2020) ψπ(s, a) = −f ′(s, a) + Ea∼π[f ′(π(a|s))− φ(a|s)] where φ(x) = kq−1 (1− x q−1) and accordingly f(x) = xφ(x). For the gradient of f(·), we have
f ′(x) = k
q − 1 (1− qxq−1)
= k
q − 1 (q − qxq−1 − (q − 1))
= qk
q − 1 (1− xq−1)− k
= qφ(x)− k.
Taking the expectation yields Tsallis entropy as follows. Ex∼π [ −f ′(x;π) + φ(x) ] = Ex∼π [ k − qφ(x) + φ(x) ] = (1− q)T kq (π) + k.
For a multivariate Gaussian distribution π, the tractable form of Ex∼π [ −f ′(x) + φ(x) ] can be derived by using that of Tsallis entropy T kq (π) of π. Thus ψπ can be rewritten as
ψπ(s, a) = qφ(s) + (q − 1)T kq (π)
In the special case of q = 1 and k = 1, we have ψπ(s, a) = log π(a|s).

heading: B.3 BREGMAN DIVERGENCE WITH TSALLIS ENTROPY REGULARIZATION
text: We consider the following form of the Bregman divergence:∫ π(x) { f ′ ( π̂(x) ) − ω ( π(x) )} dx− ∫ π̂(x) { f ′ ( π̂(x) ) − ω ( π̂(x) )} dx
For ω(x) = kq−1 (1− x q−1), f ′(x) = kq−1 (1− qx q−1) = qω(x)− k, and k = 1, the above form is equal to ∫
π(x)
[ 1− qπ̂(x)q−1
q − 1
] dx− Tq(π)− (q − 1)Tq(π̂) + 1
= 1 q − 1 − q q − 1
∫ π(x)π̂(x)q−1 dx− Tq(π)− (q − 1)Tq(π̂) + 1
= q q − 1 − q q − 1
∫ π(x)π̂(x)q−1 dx− Tq(π)− (q − 1)Tq(π̂).
Let us define two multivariate Gaussians as follows:
π(x) = N (x;µ,Σ), µ = [µ1, · · · , µd]T,Σ = Ldiag(σ21 , · · · , σ2d)LT, π̂(x) = N (x; µ̂, Σ̂), µ̂ = [µ̂1, · · · , µ̂d]T, Σ̂ = L̂diag(σ̂21 , · · · , σ2d)L̂T,
where L and L̂ denote unit lower triangular matrix. We have∫ π(x)π̂(x)q−1 dx = I(π, π̂; 1, q − 1) = exp { F (θ′)− F (θ)− (q − 1)F (θ̂) } ,
where
θ = [ Σ−1µ − 12Σ −1 ] θ̂ = [ Σ̂−1µ
− 12 Σ̂ −1 ] θ′ = θ + (q − 1)θ̂ =
[ Σ−1µ+ (q − 1)Σ̂−1µ − 12 (Σ −1 + (q − 1)Σ̂−1) ] = [ θ′1 θ′2 ] and
F (θ) = 1
2 µTΣ−1µ+
1 2 ln(2π)d|Σ| = 1 2 (µ)TΣ−1µ+ d∑ i=1 ln 2π 2 + lnσi,
F (θ̂) = 1
2 µ̂TΣ̂−1µ̂+
1 2 ln(2π)d|Σ̂| = 1 2 (µ̂)TΣ̂−1µ̂+ d∑ i=1 ln 2π 2 + ln σ̂i,
F ( θ + (q − 1)θ̂ ) = −1
4 (θ′1) T(θ′2) −1(θ′1) +
1 2 ln|−π(θ′2)−1|
We can replace some difficult computations using LDL decomposition such as Σ−1 = L−1 diag(1/σ)(L−1)T and |Σ| = ∑d i=1 lnσi.
C IMPLEMENTATION DETAILS
Unnormalized rewards of the IRL algorithm often mislead the agent to take unnecessary awareness of termination in finite-horizon MDPs (Kostrikov et al., 2018). To solve the issue IRL algorithms need to remove the difference between regarding steps depending on MDP’s time. Doob’s optimal stopping theorem formally states that the expected value of a martingale at a stopping time is equal to its initial expectation. Assume a martingale makes the entire procedure as a fair game on average, which means nothing can be gained by stopping the play.
Theorem 4 (Doob’s optimal stopping theorem). Let a process {Xt}∞t=1 be a martingale and τ be a stopping time with respect to filtration {Ft}t≥1. Assume that one of the conditions holds: (a) τ is almost surely bounded, i.e., there exists a constant c ∈ N such that τ ≤ c. (b) τ has finite expectation and the conditional expectations of the absolute value of the martingale
increments almost surely bounded, more precisely, E[τ ] <∞ and there exists a constant c such that E [ |Xt+1 −Xt| ∣∣Ft] ≤ c almost surely on the event {τ > t} for all t ≥ 0. (c) There exists a constant c such that |Xmin{t,τ}| ≤ c almost surely for all t ≥ 0. Then Xτ is an
almost surely well-defined random variable and E[Xτ ] = E[X0]. Then X∞ is integrable and E[X∞] = E[X0]
Doob’s optimal stopping theorem states one of the necessary conditions of IRL reward of normalizing the reward measures and making them as a martingale even for finite-horizon benchmarks.

heading: C.1 NETWORK ARCHITECTURES
text: For all networks, we use networks with 2-layer MLP with 100 hidden units, respectively. We consider the reward model with two separate neural networks (ψϕ, b(·)) for the proposed reward function for λ > 0:
rϕ(s, a) = ψ λ ϕ(s, a) = λψϕ(s, a) + b(s),
Motivated by RAIRL-DBM, we consider the reward models in Figure 9. The model outputs reward for proximal updates trained by mirror descent and state-only discriminator network. Discriminating state visitation by b(·) is required because the reward function needs to consider every state (especially the state that cannot be visited by πE) until DKL(ρπ∥ρπE ) ≈ 0. Figure 9 (a) shows logits of the softmax distribution involved when calculating rewards when the action space is discrete. For continuous control (Figure 9 (b)), the architecture is similar, where the mean and covariance are used to compute a reward for a certain action.

heading: C.2 MULTI-GOAL ENVIRONMENT
text: Let the 2D coordinate denote the position of a point mass on the environment. In the multi-goal environment, the agent is initially located according to the normal distribution N (0, (0.1)2I). The four goals are located at (6, 0), (−6, 0), (0, 6), and (0,−6), where the agent can move a maximum of
1 unit per timestep for each coordinate. The ground-truth reward is given by the difference between successive values of a Gaussian mixture depicted as Figure 10.
We used the full covariance Gaussian distribution in this experiment (as well as the conceptual experiment in Figure 2. Note that the covariance matrix is positive-definite and symmetric. To achieve numerically stable computation, we applied LDL decomposition to covariance matrix involving unit lower- and upper-triangle matrices, and a diagonal matrix. As a result, the policy network outputs a vector (µ(s)T,σ(s)T, l(s)T)T where the additional part l(s) denotes d(d−1)2 entries of unit lower triangular matrix. Denote L(s) as a unit lower triangular matrix from l(s). For example, the covariance matrix can be reconstructed by
Σ(s) = L(s)[diag(σ(s))]L(s)T.
Additionally, the action samples can be calculated by
a = µ(s) +L(s)(σ(s) · z) z ∼ N (0, I)
Computing inverses, determinants and multiplications with unit triangular matrices and diagonal matrices can be efficiently performed by numerical libraries. Therefore, we can fully model the Bregman divergence and reward using neural networks as provided in Appendix B.

heading: C.3 MUJOCO EXPERIMENTS
text: Instead of directly using squashed policies proposed in SAC (Haarnoja et al., 2018), we assume the application of tanh as a part of the environment (known as hyperbolized environments of RAIRL (Jeon et al., 2020)). Specifically, after an action a is sampled from the policies, we pass tanh(a/1.01) ∗ 1.01 to the environment. We additionally clip the hyperbolized actions to 1, if the environment is not tolerant to the excessive values of action. Therefore, we can consider the standard diagonal Gaussian policy π(·|s) = N (µ(s),Σ(s)) where µ(s) denotes means [µ1(s), µ1(s), . . . , µd(s)]T and σ(s) denotes standard variance [σ1(s), σ2(s), . . . , σd(s)]
T. We soft-clip the standard deviation as σi(s) ∈ [ln 0.01, ln 2] for the stability using tanh. Last, we update the moving mean of intermediate values of regularized reward r(s, a)+λφ(π(a|s)) and update the RL algorithm with mean-subtracted rewards. In addition to the analyses of (Fu et al., 2017; Jeon et al., 2020) regarding reward shaping and normalization, meanzero rewards for training agents have the additional property of preventing termination awareness, as stated by the optimal stopping theorem (Theorem 4).

heading: C.4 HYPERPARAMETERS
text: Tables 2, 3, and 4 show the hyperparameters of the conducted experiments.

Review summaries:
This paper presents a novel mirror-descent adversarial inverse reinforcement learning (MD-AIRL) algorithm. MD-AIRL considers the reward function as an iterative sequence in a proximal method. MD-AIRL has been introduced with dense theoretical analysis and validated with diverse experiments covering both discrete and continuous action spaces. 

This paper proposes an alternative optimization technique (Mirror Descent) for adversarial inverse reinforcement learning, claiming it will resolve some of the issues with previous methods.   -The authors prove the convergence of the proposed optimization method and provide a regret bound.  -They experimentally show the benefits of the proposed methods in discrete and continuous action spaces.   

The paper proposes an algorithm MD-AIRL, where each iteration involves policy and reward updates, and an expert policy estimation. The algorithm is motivated by MD, especially the reward update is modeled as a projected update of an MD update. The sufficiency and necessities of the choices of stepsizes to guarantee convergence are theoretically justified. A practical algorithm with numerical experiments are also provided.