Year: 2019

Abstract:
We consider the Inverse Reinforcement Learning problem in Contextual Markov Decision Processes. In this setting, the reward, which is unknown to the agent, is a function of a static parameter referred to as the context. There is also an “expert” who knows this mapping and acts according to the optimal policy for each context. The goal of the agent is to learn the expert’s mapping by observing demonstrations. We define an optimization problem for finding this mapping and show that when it is linear, the problem is convex. We present and analyze the sample complexity of three algorithms for solving this problem: the mirrored descent algorithm, evolution strategies, and the ellipsoid method. We also extend the first two methods to work with general reward functions, e.g., deep neural networks, but without the theoretical guarantees. Finally, we compare the different techniques empirically in driving simulation and a medical treatment regime.

heading: 1 INTRODUCTION
text: We study sequential decision-making in a Contextual Markov Decision Process (CMDP, Hallak et al. (2015)), where the reward, while unknown to the agent, depends on a static parameter referred to as the context. For a concrete example, consider the dynamic treatment regime (Chakraborty & Murphy, 2014). Here, there is a patient and a clinician which acts to improve the patient’s health. The context is composed of static information of the patient (such as age and weight); the state is composed of the patient’s dynamic measurements (such as heart rate and blood pressure); and the clinician’s actions are a set of intervention categories (e.g., infusion). The reward is different for each patient (context), and there is a mapping from the context to the reward.
Recent trends in personalized medicine motivate this model – instead of treating the ”average patient”, patients are separated into different groups for which the medical decisions are tailored (Fig. 1b). For example, in Wesselink et al. (2018), the authors study organ injury, which may occur when a specific measurement (mean arterial pressure) decreases below a certain threshold. They found that this threshold varies across different patient groups (context). In other examples, clinicians set treatment goals for the patients, i.e., they take actions to make the patient measurements reach some pre-determined values. For instance, in acute respiratory distress syndrome (ARDS), clinicians argue that these treatment goals should depend on the static patient information (the context) (Berngard et al., 2016).
There are serious issues when trying to manually define a reward signal in real-world tasks. When treating patients with sepsis, for example, the only available signal is the mortality of the patient at the end of the treatment (Komorowski et al., 2018). This signal is sparse, and it is unclear how to manually tweak the reward to maximize the patient’s health condition (Leike et al., 2017; Raghu et al., 2017; Lee et al., 2019).
To address these issues, we propose the Contextual Inverse Reinforcement Learning (COIRL) framework. Similarly to Inverse Reinforcement Learning (Ng & Russell, 2000, IRL), we focus on trying to infer the mapping from contexts to rewards by observing experts. The main challenge in our problem is that for each context there is a different reward, hence, a different optimal policy for each context. Therefore, Apprenticeship Learning algorithms (Abbeel & Ng, 2004; Syed & Schapire, 2008) that try to mimic the expert cannot be used and, instead, we focus on directly learning the mapping.
In particular, our main contributions are:
1. We formulate COIRL with a linear mapping as a convex optimization problem. 2. We propose and analyze the sample complexity of three algorithms for COIRL: the mirrored
descent alg. (MDA), evolution strategies (ES), and the ellipsoid method. 3. For nonlinear mappings, we implement a deep learning version for MDA and ES (without
theoretical guarantees). 4. We compare these methods empirically on two frameworks: an autonomous driving simulator
(Abbeel & Ng, 2004) and a dynamic treatment regime (Komorowski et al., 2018).

heading: 2 PRELIMINARIES
text: Contextual MDPs: A Markov Decision Process (Puterman, 1994, MDP) is defined by the tuple (S,A, P, ξ, R, γ) where S is a finite state space, A a finite action space, P : S × S ×A→ [0, 1] the transition kernel, ξ the initial state distribution, R : S → R the reward function and γ ∈ [0, 1) is the discount factor. A Contextual MDP (Hallak et al., 2015, CMDP) is an extension of an MDP, and is defined by (C,S,A,M, γ) where C is the context space, andM is a mapping from contexts c ∈ C to MDPs: M(c) = (S,A, P,Rc, ξ, γ). In addition, each state is associated with a feature vector φ : S → [0, 1]k. Note that P and ξ are not context dependent. We consider a setting in which the reward for context c is a linear combination of the state features: R∗c(s) = f
∗(c)Tφ(s). The goal is to approximate f∗(c) using a function fW (c), with parameters W . This notation allows us to present our algorithms for any function approximator fW (c), and in particular, a deep neural network (DNN). For the theoretical analysis, we will further assume a linear setting, where f∗(c) = cTW ∗, fW (c) = cTW and that W ∗ is in some convex setW. We assume that c ∈ C = ∆d−1, the standard d− 1 dimensional simplex. This assumption makes the contexts bounded (which we use in our proofs), and it also allows a straight-forward expansion to a model in which the transitions are also a linear mapping of the context (Modi et al., 2018). One way of viewing this model is that each row in the mapping W ∗ is a base rewards coefficient vector, and the reward for a specific context is a convex combination of these base rewards.
We consider deterministic policies π : S → A which dictate the agent’s behaviour at each state. The value of a policy π for reward coefficients vector r is: V πr = Eξ,P,π[ ∑∞ t=0 γ tR(st)] = r Tµ(π)
where µ(π) := Eξ,P,π[ ∑∞ t=0 γ
tφ(st)] ∈ Rk is called the feature expectations of π. For the optimal policy with respect to (w.r.t.) a reward coefficients vector r, we denote the value by V ∗r . For any context c, π∗c denotes the optimal policy w.r.t. reward R ∗ c(s) = f
∗(c)Tφ(s) and π̂c(W ) denotes the optimal policy w.r.t. reward R̂c(s) = fW (c)Tφ(s).
Inverse Reinforcement Learning in CMDPs: In standard IRL, the goal is to learn a reward which best explains the behavior of an observed expert. The model describing this scenario is the MDP\R -
an MDP without a reward function (also commonly called a controlled Markov chain). Similarly, we denote a CMDP without a mapping of context to reward by CMDP\M. The goal in Contextual IRL is to approximate the mapping f∗(c) by observing an expert. The expert knows f∗(c), and for each context c, can provide a demonstration from π∗c .
Contextual dynamics: Learning a transition kernel and an initial state distribution that is parametrized by the context is an orthogonal problem to COIRL. Therefore, we focus only on a contextual reward which simplifies our analysis. Existing methods, such as in Modi et al. (2018), can be used to learn the mappings for the transition kernel and initial distribution in a contextual model. In conjunction with the simulation lemma (Kearns & Singh, 2002), these methods can extend our results to the more general CMDP setting.

heading: 3 OPTIMIZATION METHODS FOR COIRL
text: In this section, we propose and analyze optimization algorithms for minimizing the following loss function; Lemma 1 below justifies its use for COIRL.
Loss(W ) = Ec max π
[ fW (c) · ( µ(π)− µ(π∗c ) )] = Ec [ fW (c) · ( µ(π̂c(W ))− µ(π∗c ) )] . (1)
Lemma 1. Loss(W ) satisfies the following properties: (1) ∀W, Loss(W ) ≥ 0, and Loss(W ∗) = 0. (2) If Loss(W ) = 0 then ∀c ∈ C, the expert policy π∗c is the optimal policy w.r.t. reward cTW.
To evaluate the loss, the optimal policy π̂c(W ) and its features expectations µ(π̂c(W )) must be computed for all contexts. For a specific context, finding π̂c(W ) can be solved with standard RL methods such as Value or Policy Iteration. Computing µ(π̂c(W )) is equivalent to policy evaluation (solving linear equations).
The challenge is that Eq. (1) is is not differentiable in W . We tackle this problem using two methods for computing descent directions that do not involve differentiation: (i) computing subgradients and (ii) randomly perturbing the loss function. In addition, as the loss is defined in expectation over the contexts, computing it requires to calculate the optimal policy for all contexts. We deal with this issue at the end of Section 3.1. In the special case that fW (c) is a linear function, Eq. (1) is convex. The following Lemma characterizes Eq. (1) in this case.
Lemma 2. Let Llin(W ) = Ec [ cTW · ( µ(π̂c(W ))−µ(π∗c ) )] . We have that: (1) Llin(W ) is a convex
function. (2) g(W ) = Ec [ c ( µ(π̂c(W ))− µ(π∗c ) )] is a sub gradient of Llin(W ). (3) Llin is a Lipschitz continuous function, with Lipschitz constant L = 21−γ w.r.t. ‖·‖∞ and L = 2 √ dk 1−γ w.r.t. ‖·‖2.
A technical proof (by definition) is provided in the supplementary material. Note that g(W ) ∈ Rd×k; we will sometimes refer to it as a matrix and sometimes as a flattened vector, no confusion will arise. Remark 1. The Lipschitz of LLin(W ) is related to the simulation lemma (Kearns & Singh, 2002); a small change in the reward results in a small change in the optimal value. Remark 2. As g(W ) is a subgradient of Loss(W ), it can be used to back-propagate DNNs. Clearly, we cannot guarantee convexity (hence no theoretical guarantees), but we can design Loss(W ) to be Lipschitz continuous in W using the methods presented in Cisse et al. (2017); Arjovsky et al. (2017). Remark 3. The subgradient g(W ) is given in expectation over contexts, and in expectation over trajectories (feature expectations). We will later see how to replace it with an unbiased estimate, which can be computed by observing a single expert trajectory for a single context.

heading: 3.1 MIRRORED DESCENT FOR COIRL
text: Lemma 2 identifies LLin(W ) as a convex function and provides a method to compute its subgradients. A standard method for minimizing a convex function over a convex set is the subgradient projection algorithm (Bertsekas, 1997): wt+1 = ProjW{wt − αtg(wt)}, where f(wt) is a convex function, g(wt) is a subgradient of f(wt), and αt the learning rate.W is a convex set, and specifically, we consider the `2 ball (Abbeel & Ng, 2004) and the simplex (Syed & Schapire, 2008)1. We focus on
1Scaling of the reward by a constant does not affect the resulting policy, thus, these sets are not restricting.
Algorithm 1 MDA for COIRL input: a convex set W , T number of iterations initialize w1 ∈ W for t = 1, . . . , T do
Observe c, µ(π∗c ) Compute π̂c(W ), µ(π̂c(W )) Compute gt according to Lemma 2 if PSGD then
αt = (1− γ) √ 1 2dkt wt+1 = wt − αtgt if ‖wt+1‖ > 1 then
wt+1 = wt+1/ ‖wt+1‖2 else if Exponential weights then
αt = (1− γ) √ log(dk) 2t for i = 1, . . . , dk do wt+1(i) = wt(i) exp (−αtgt(i))
wt+1 = wt+1/ ∑ i wt+1(i)
return 1t ∑T t=1 wt
Algorithm 2 ES for COIRL
input: step sizes {αt}Tt=1s, noise STD σ, number of evaluations m and smoothing parameter ν > 0 initialize: W ∈ Rk for t = 1, . . . , T do
Observe c, µ(π∗c ) for j = 1, ...,m do
uj ∼ N k(0, σ2) Lossj(W ) = Loss ( W +
uj ||uj ||ν ) dLoss(W ) = ∑m j=1 Lossj(W ) uj ||uj ||ν If Loss(W − αtb dLoss(W )) < Loss(W ) then W = W − αtmσdLoss(W )
return W
a generalization of the subgradient projection algorithm that is called the mirror descent algorithm (Nemirovsky & Yudin, 1983, MDA): wt+1 = arg minw∈W { w · ∇f (wt) + 1αtDψ(w,wt) } , where Dψ(w,wt) is a Bregman distance2, associated with a strongly convex function ψ. The following theorem characterizes the convergence rate of MDA.
Theorem 1 (Convergence rate of MDA). Let ψ be a σ-strongly convex function onW w.r.t. ‖·‖, and let D2 = supw1,w2∈W Dψ(w1, w2). Let f be convex and L-Lipschitz continuous w.r.t. ‖·‖. Then,
MDA with αt = DL √ 2σ t satisfies: f ( 1 T ∑T s=1 xs ) − f(x∗) ≤ DL √ 2 σT .
We refer the reader to Beck & Teboulle (2003) and Bubeck (2015) for the proof. Next, we provide two MDA instances (see, for example Beck & Teboulle (2003) for derivation) and analyze them for COIRL.
Projected subgradient descent (PSGD): LetW be an `2 ball with radius 1. Fix || · ||2, and ψ(w) = 1 2 ||w|| 2 2. ψ is strongly convex w.r.t. || · ||2 with σ = 1. The associated Bregman divergence is given by Dψ(w1, w2) = 0.5||w1 − w2||22. Thus, mirror descent is equivalent to PSGD. D2 = maxw1,w2∈W Dψ(w1, w2) ≤ 1, and according to Lemma 2, L = 2 √ dk 1−γ . Thus, we have that after T
iterations Llin ( 1 T ∑T t=1 wt ) − Llin(w∗) ≤ O ( √ dk (1−γ) √ T ) . Exponential Weights (EW): Let W be the standard dk − 1 dimensional simplex. Let ψ(w) =∑ i w(i) log(w(i)). ψ is strongly convex w.r.t. || · ||1 with σ = 1. We get that the associated Bregman
divergence is given by Dψ(w1, w2) = ∑ i w1(i) log( w1(i) w2(i)
), also known as the Kullback-Leibler divergence. In addition, D2 = maxx,y∈W Dψ(w1, w2) ≤ log(dk) and according to Lemma 2, L = 21−γ . Furthermore, the projection onto the simplex w.r.t. to this distance amounts to a simple renormalization w ← w/||w||1. Thus, we get that MDA is equivalent to the exponential weights algorithm and Llin ( 1 T ∑T t=1 wt ) − Llin(w∗) ≤ O (√log(dk) (1−γ) √ T ) .
Practical MDA: One of the “miracles” of MDA is its robustness to noise. If we replace gt with an unbiased estimate g̃t, such that Eg̃t = gt and E ‖g̃t‖ ≤ L, we obtain the same convergence results as in Lemma 2 (Robbins & Monro, 1951) (see, for example, (Bubeck, 2015, Theorem 6.1)). Such an unbiased estimate can be obtained in the following manner: (i) sample a context ct, (ii) compute
2We refer the reader to Appendix C for definitions of the Bregman distance, the dual norm, etc.
µ(π∗ cTt wt ), (iii) observe a single expert demonstration τEi = {si0, a0, si1, a1, . . . , }, where ai is chosen by the expert policy π∗
cTt w ∗ (iv) let µ̂i = ∑ t∈[0,...,|τEi |−1]
γtφ(sit) be the accumulated discounted features across the trajectory such that Eµ̂i = µ(π∗c ).
The challenge is, that for µ̂i to be an unbiased estimate of µ(π∗cTt w∗), τ E i needs to be of infinite length. There are two ways in which we can tackle this issue. We can either (1) execute the expert trajectory online, and terminate it at each time step with probability 1− γ (as in (Kakade & Langford, 2002)), or (2) execute a trajectory of length H = 11−γ log(1/ H). The issue with the first approach is that since the trajectory length is unbounded, the estimate µ̂i cannot be shown to concentrate to µ(π∗c ) via Hoeffding type inequalities. Nevertheless, it is possible to obtain a concentration inequality using the fact that the length of each trajectory is bounded in high probability (similar to Zahavy et al.). The second approach can only guarantee that ‖gt − Eg̃t‖ ≤ H (Syed & Schapire, 2008). Therefore, using the robustness of MDA to adversarial noise (Zinkevich, 2003), we get that MDA converges with an additional error of H , i.e., Llin ( 1 T ∑T t=1 wt ) − Llin(w∗) ≤ O ( 1√ T ) + H . While this sampling mechanism comes with the cost of a controlled bias, usually it is more practical, in particular when the trajectories are given as a set demonstrations (offline data).

heading: 3.2 EVOLUTION STRATEGIES FOR COIRL
text: To minimize Eq. (1), we also design a derivative free algorithm (Algorithm 2) that is based on Evolution Strategies (Salimans et al., 2017, ES). For convex optimization problems, ES is a gradientfree descent method that is based on computing finite differences (Nesterov & Spokoiny, 2017), whose sample complexity is provided below in Theorem 2. The Theorem is given in terms of the Lipschitz constant, which is upper bounded by 2 √ dk
1−γ (Section 3.1). While this approach has looser upper-bound guarantees compared to MDA (Theorem 1), Nesterov & Spokoiny (2017) observed that in practice, it often outperforms subgradient based methods. Thus, we test this method empirically and compare it with the subgradient method (Section 3.1). ES is also known to perform well in practice, even with nonconvex objectives. Specifically, Salimans et al. (2017) has shown that ES can be used to optimize the parameters of a DNN to solve challenging high dimensional RL tasks like playing Atari.
Theorem 2 (ES Convergence Rate (Nesterov & Spokoiny, 2017)). Let Llin(W ) be a non-smooth convex function with Lipschitz constant L, such that ||x0 − x∗|| ≤ D, step size of αt =
D (dk+4) √ T+1L and ν ≤ 2L √ dk
then in T = 4(dk+4) 2D2L2
2 ES finds a solution which is bounded by EUT−1 [Llin(x̂T )]− Llin(x∗) ≤ , where UT = {u0, . . . , uT } denotes the random variables of the algorithm up to time T and x̂T = arg mint=1,...,T Llin(xt).

heading: 3.3 ELLIPSOID ALGORITHMS FOR COIRL
text: Algorithm 3 Ellipsoid algorithm for COIRL Initialize: Θ0 ← B∞(0, 1) ={x ∈ Rd·k : ||x||∞ ≤ 1} Θ1 ← MVEE(Θ0) for t = 1, 2, . . . do
Observe ct, let W t be the center of Θt Play episode using π̂t = arg maxπ V
π cTt Wt
if V ∗ cTt W ∗ − V π̂tcTt W∗ > then µ(π∗ct) is revealed Let at = ct ( µ(π∗ct)− µ(π̂t)
) Θt+1← MVEE ({ θ ∈ Θt : θTat ≥WTt at
}) else
Θt+1 ← Θt The final algorithm we consider is an ellipsoid method, introduced to the IRL setting by Amin et al. (2017). In this section we extend it to the contextual setting, specifically, we focus on finding a linear mapping W and further assume that W = {W : ||W ||∞ ≤ 1}, and that W ∗ ∈ W . The algorithm maintains an ellipsoidshaped feasibility set that contains W ∗. At any step, the current estimation Wt of W ∗ is defined as the center of the ellipsoid, and the agent acts optimally w.r.t. this estimation. If the agent performs suboptimally, the expert provides a demonstration in the form of the optimal feature expectations for ct, µ(π∗ct). The feature expectations are used to generate a linear constraint (hyperplane) on the ellipsoid that is crossing its center. Under this constraint, we construct a new feasibility set that is half of the
previous ellipsoid, and still contains W ∗. For the algorithm to proceed, we compute a new ellipsoid that is the minimum volume enclosing ellipsoid (MVEE) around this ”half-ellipsoid” 3. These updates are guaranteed to gradually reduce the volume of the ellipsoid (a well-known result (Boyd & Barratt, 1991)) until its center is a mapping which induces -optimal policies. Theorem 3 shows that this algorithm achieves a polynomial upper bound on the number of sub-optimal time-steps. Finally, note that in Algorithm 3 we use an underline notation to denote a ”flattening” operator for matrices, and to denote a composition of an outer product and the flattening operator. The proofs in this section are provided in the supplementary material, and are adapted from (Amin et al., 2017). Theorem 3. In the linear setting where R∗c(s) = cTW ∗φ(s), for an agent acting according to Algorithm 1, the number of rounds in which the agent is not -optimal is O(d2k2 log( dk(1−γ) )). Remark 4. Note that the ellipsoid method presents a new learning framework, where demonstrations are only provided when the agent performs sub-optimally. Thus, the theoretical results in this section cannot be directly compared with those of the descent methods. We further discuss this in the experiments and discussion sections. Remark 5. The ellipsoid method does not require a distribution over contexts - an adversary may choose them. MDA can also be easily extended to the adversarial setting via known regret bounds on online MDA (Hazan, 2016).
Practical ellipsoid algorithm: In many real-world scenarios, the expert cannot evaluate the value of the agent’s policy and cannot provide its policy or feature expectations. To address these issues, we follow Amin et al. (2017) and consider a relaxed approach, in which the expert evaluates each of the individual actions performed by the agent rather than its policy, and provides finite rollouts instead of a policy or feature expectations (see the supplementary material (Algorithm 4) for pseudo code). We define the expert criterion for providing a demonstration to be Q∗
cTt W ∗(s, a) + < V ∗ cTt W ∗(s) for each state-action pair (s, a) in the agent’s trajectory.
Near-optimal experts: In addition, we relax the optimality requirement of the expert and instead assume that, for each context ct, the expert acts optimally w.r.t. W ∗t which is close to W
∗; the expert also evaluates the agent w.r.t. this mapping. This allows the agent to learn from different experts, and from non-stationary experts whose judgment and performance slightly vary over time. If a sub-optimal action w.r.t. W ∗t is played at state s, the expert provides a roll-out of H steps from s to the agent. As this roll-out is a sample of the optimal policy w.r.t. W ∗t , we aggregate n examples to assure that with high probability, the linear constraint that we use in the ellipsoid algorithm does not exclude W ∗ from the feasibility set. Note that these batches may be constructed across different contexts, different experts, and different states from which the demonstrations start. Theorem 4 below upper bounds the number of sub-optimal actions that Algorithm 4 chooses.4
Theorem 4. For an agent acting according to Algorithm 4 , with probability of at least 1 − δ, for H = d 11−γ log( 8k (1−γ) )e and n = d 512k2 (1−γ)2 2 log(4dk(dk + 1) log( 16k √ dk
(1−γ) )/δ)e, if ∀t : W ∗t ∈ B∞(W
∗, (1−γ) 8k ) ∩ Θ0 the number of rounds in which a sub-optimal action is played is O ( d2k4 (1−γ)2 2 log ( dk (1−γ)δ log( dk (1−γ) ) )) .

heading: 4 EXPERIMENTS
text: The simulations in this section include two domains: (1) an autonomous driving simulation (Abbeel & Ng, 2004), that we adapted to the contextual setup and (2) a medical treatment regime, constructed from a data set of expert (clinician) trajectories for treating patients with sepsis5. In each of these domains we compare the algorithms in two setups: the ellipsoid learning framework and an offline framework. All the results are averaged across 10 random seeds in Section 4.1 and 5 seeds in Section 4.2 (we report the mean and the standard deviation). Due to space considerations we present the simulations in the ellipsoid framework only for the car domain, and the simulations in the offline framework only in the dynamic treatment regime. Complementary simulations can be found in the supplementary material.
3This procedure follows a sequence of linear algebra operations which we explain in the appendix. 4MDA also works with near optimal experts due to the robustness of MDA. The analysis of this case is
identical to the analysis of biased trajectories, as we discuss in the end of Section 3.1. 5The data, code and implementation of our algorithms can be found in github.com/CIRLMDP/CIRL.

heading: 4.1 DRIVING SIMULATION – THE ELLIPSOID FRAMEWORK
text: In the ellipsoid framework, an expert evaluates the agent policy. If the agent’s policy is sub-optimal, the expert provides the agent its feature expectations; otherwise, no demonstration is given. The algorithm performs learning in between demonstrations. This setup enables a proper comparison with the ellipsoid algorithm, which requires the additional expert supervision. We measure performance w.r.t. the following criteria: (1) # demonstrations – the amount of contexts on which each algorithm requested an expert demonstration (y-axis) as a function of time, i.e., the total number of contexts (x-axis). (2) Value – the difference in value, between the agent policy and the expert policy w.r.t. the true reward mapping, i.e., ∑ c∈Ctest fW∗(c) · ( µ(π̂c(W ))− µ(π∗c ) ) , where Ctest is a holdout (test) set of contexts. The x-axis measures the amount of demonstrations given.
Setup. This domain simulates a three-lane highway with two visible cars - cars A and B (illustration provided in the appendix). The agent, controlling car A, can drive both on the highway and off-road. Car B drives on a fixed lane, at a slower speed than car A. Upon leaving the frame, car B is replaced by a new car, appearing in a random lane at the top of the screen. The feature vector φ(s) is composed of 3 features: (1) a speed feature, (2) a collision feature, which is valued 0 in case of a collision and 0.5 otherwise, and (3) an off-road feature, which is 0.5 if the car is on the road and 0 otherwise.
In this task, the context vector implies different priorities for the agent; should it prefer speed or safety? Is going off-road to avoid collisions a valid option? For example, an ambulance will prioritize speed and may allow going off-road as long as it goes fast and avoids collisions, while a bus will prioritize avoiding both collisions and off-road driving as safety is its primary concern. To demonstrate the effectiveness of our solutions, the mapping f : C 7→ [−1, 1]k is constructed in a way that induces different behaviors for different contexts, making generalization a challenging task. We provide additional details on the domain as well as the hyper parameter selection in the appendix.
Linear: The optimal behavior is defined using a linear mapping W ∗. In this setting, all three approaches obtain competitive results, in terms of generalization, although the ES is capable of obtaining these results faster, as is seen through the regret and number of required demonstrations.
Nonlinear: For the nonlinear task, we consider two reward coefficient vectors r1 and r2, and define the mapping by f∗(c) = r1 if ||c||∞ ≥ 0.55, and r2 otherwise - an illustration is provided in the appendix. In order to learn the nonlinear mapping, we represent fW (c) using a DNN, a multi-layered perceptron, which maps from context to reward vector. DNNs have proven to be capable of extracting meaningful features from complex high-dimensional data, e.g., images - in these scenarios, the linear assumption no longer holds, yet DNNs often overcome such issues. In this setting, the superiority of the descent methods rises; as the linear assumption in the ellipsoid algorithm is not met, it fails to generalize and keeps requiring new demonstrations. We believe these results to be crucial when considering real-life applications, in which the problem is not necessarily linear. Such cases highlight the strength of the descent methods, which, as Fig. 2 shows, are capable of scaling to nonlinear high dimensional mappings.

heading: 4.2 DYNAMIC TREATMENT REGIME – THE OFFLINE FRAMEWORK
text: In the offline framework, we focus on the ability to learn from previously collected data. A data set of previously collected trajectories is given, such that a single trajectory of finite length is observed for each context and no context is observed more than once. We measure performance w.r.t. the following
criteria: (1) Value – as in the ellipsoid framework above, but here the x-axis corresponds to the amount of iterations. Each iteration corresponds to a single subgradient step, where the subgradient is computed from a mini batch of 10 contexts. (2) Loss – as in Eq. (1). (3) Accuracy % – the percent of actions on which the expert and the agent agree on. All these criteria are evaluated on a holdout set.
Setup. In the dynamic treatment regime, there is a clinician which acts to improve a sick patient’s medical condition. The context (static information) represents patient features, which do not change during treatment, such as age and gender. The state summarizes the dynamic measurements of the patient, e.g., blood pressure and EEG readouts. The actions are the forms of intervention a clinician may take, including combinations of various treatments provided in parallel. Dynamic treatment regimes are particularly useful for managing chronic disorders and fit well into the broader paradigm of personalized medicine (Komorowski et al., 2018; Prasad et al., 2017).
The agent needs to choose the right treatment for a patient that is diagnosed with sepsis. We use the MIMIC-III data set (Johnson et al., 2016) and follow the data processing steps that were taken in Jeter et al. (2019). As performing off-policy evaluation is not possible using this data set, due to it not satisfying basic requirements (Gottesman et al., 2018; 2019), we designed a simulator of a CMDP. The simulator is based on this data set; a complete overview and explanation on how it was created is provided in the appendix. The mapping W ∗ is linear, W ∗ ∈ R8×42, which we constructed from the data. In the simulator, the expert acts optimally w.r.t. this W ∗.
Specifically, when treating a sepsis patient, the clinician has several decisions to make, such as whether or not to provide a patient with vasopressors, drugs which are commonly provided to restore and maintain blood pressure in patients with sepsis. However, what is regarded as healthy blood pressure differs based on the age and weight of the patient (Wesselink et al., 2018). In our setting, W captures this information - as it maps from contextual (e.g., age) and dynamic information (e.g., blood pressure) to reward.
Results. Fig. 3 presents the ability of the descent methods to generalize to unseen contexts by learning from offline data (without supervision). The data is composed of a set of trajectories, i.e., offline data, that were collected from experts (clinicians treating patients). In each iteration, we sample a minibatch of 10 contexts, i.i.d, from the context distribution. For each context, there is a corresponding expert trajectory of length H = 40. Performance is measured on a holdout set of 300 contexts (that are sampled from the same context distribution) according to Theorem 1. We can see that both ES and PSGD attain near-optimal performance using only previously collected expert trajectories.
Looking at Fig. 3a, we can see that all the algorithms manage to minimize the loss to roughly the same error. The small bias is explained by the fact that we use truncated trajectories (as we discussed in the practical MDA paragraph) where in the ellipsoid framework experiments we used feature expectations. We can also see that minimizing the loss leads to policies that attain −optimal value w.r.t. the true reward Fig. 3b. Finally, in Fig. 3 we can see that all the algorithms reach around 70% accuracy with the expert policy. We emphasize here that 100% accuracy should not be expected for two reasons: (i) different policies may have the same feature expectations (hence the same value) but make different decisions (ii) there exists reward for which there is more than one optimal policy. Nevertheless, Fig. 3 suggests that accuracy is correlated with minimizing the COIRL loss (Eq. (1)).
Finally, we present results in the nonlinear setting. Here, there is a non-linear function of the context that determines which one of the two reward coefficient vectors is used, i.e., f∗(c) = r1 if age > 0.1, and r2 otherwise. Where age refers to the normalized age of the patient, which is an element
of the context vector. We use a DNN to learn the mapping and follow Section 3.1 (PSGD). As seen in Fig. 4, the PSGD algorithm minimizes the loss and achieves a value that is close to that of the expert. In addition, similarly to Fig. 7, accuracy and performance do not necessarily correlate one to another.

heading: 5 RELATED WORK
text: We begin with a short discussion on contextual policies, i.e., a policy that is a function of both the state and the context. While there is empirical evidence that learning such a policy may perform well in practice (e.g., Xu et al. (2018); Fu et al. (2019)), from a theoretical point of view, there exist hardness results in this setting. Specifically, given an MDP with k + 1 states, there is a reduction from the training problem of the union of k hyperplanes to the policy (see appendix E for proof).
Alternatively, one may consider applying an AL algorithm on a single, large MDP that includes all the states and the contexts. For a concrete example, consider a reduction from the CMDP model to a large MDP where each state is expanded by concatenating the context to it. The new states are s′ = (s, c), and the new features are φ(s′) = c φ(s). Generally speaking, applying an AL algorithm to this large MDP will give the same scalability and sample complexity as COIRL. However, as the large MDP has |S ′| = |C||S| states, computing the optimal policy in each iteration of the algorithm will require at least |C| times more time. To illustrate this problem we conducted a simple grid world experiment on a 3× 4 grid world MDP (with 12 states) and one-hot features (φ(si) = ei ∈ R12). The dynamics are deterministic, and the actions correspond to going up, down, left and right (with cyclic transitions on the borders). The contexts correspond to ”preferences” on the grid; mathematically, each context is sampled from a uniform distribution over the simplex. The mapping W ∗ is set to be I12·12, and for AL, we let w∗ to be a flattened version of W ∗.
We compare the performance of PSGD with the projection algorithm of (Abbeel & Ng, 2004) in Fig. 5. We measure performance by three metrics: run-time, value, and accuracy. Inspecting the results, we can see that AL in the large MDP requires significantly more time to run as the number of contexts grows, while the run time of PSGD (COIRL) is not affected by the number of contexts. We can also see that both methods achieve roughly the same performance: COIRL performs slightly better in terms of value while AL performs slightly better in terms of accuracy. To conclude, AL on a large MDP does not scale to problems with large context spaces. In addition, this construction is
only possible when there is a finite number of context and does not provide generalization results. We avoided all of these issues in the COIRL framework.

heading: 6 SUMMARY AND DISCUSSION
text: In this work, we formulated and studied the COIRL problem. We presented two types of algorithms to solve it: (1) cutting plane methods (ellipsoid) and (2) iterative descent approaches (MDA and ES). We summarize the theoretical guarantees of the different algorithms in Table 1.
We can see that the iterative descent approaches have better dependence in dk than the ellipsoid method, i.e., they scale better with the dimensions of the problem. In particular, the EW algorithm has a logarithmic dependence in dk, which makes it computationally comparable to standard IRL/AL (on a single, noncontextual MDP). In addition, the iterative methods extend naturally to the more general scenario where the mapping from contexts to rewards is not linear, and fW is modeled as a DNN. As Sutton (2019) puts it: ”The biggest lesson that can be read from 70 years of AI research is that general methods that leverage computation are ultimately the most effective, and by a large margin”.
The ellipsoid method has better sample complexity (as a function of ) than the descent methods in the deterministic setting. However, both methods attain the same complexity in the more realistic, stochastic setting. Our empirical findings suggest that the iterative methods always outperform the ellipsoid algorithm. Among these methods, we found the ES method to perform better than the MDA method. Similar findings were reported in (Nesterov & Spokoiny, 2017) for other convex problems.
The iterative methods have another advantage over the ellipsoid method – they can learn from previously collected demonstrations (i.e., offline learning). The ellipsoid framework, on the other hand, requires expert supervision throughout the entire learning process.
Finally, an attractive property of the ellipsoid learning framework is its safety, i.e., an IRL algorithm that is being supervised by an expert will never perform sub-optimally. In each step, either that the agent performs -optimally or that the expert acts on its behalf (provides a demonstration). This property is appealing in mission-critical domains where errors have a high cost; for instance, in health-care, a failure may result in a loss of lives. In the experimental section, we have seen that we can use this learning framework for the iterative methods as well while enjoying improved efficiency.

heading: Appendices

heading: CONTENTS

heading: A Complementary simulations 14
text: A.1 Autonomous driving simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
A.2 Dynamic treatment regime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

heading: B Experimental details 15
text: B.1 Autonomous driving simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
B.2 Dynamic treatment regime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

heading: C Proofs for Section 3 18

heading: D Proofs & pseudo code for Section 3.3 20
text: D.1 Ellipsoid Algorithm for trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . 20
D.2 MVEE computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
D.3 Proof of Theorem 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
D.4 Proof of Theorem 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
E Hardness of learning contextual policies 23

heading: A COMPLEMENTARY SIMULATIONS

heading: A.1 AUTONOMOUS DRIVING SIMULATION
text: Similar to Section 4.2, we compare the various methods in the offline framework. We can see that all the algorithms manage to minimize the loss, and achieve near optimal value. We can also see that they achieve high accuracy with the expert policy but not 100%.

heading: A.2 DYNAMIC TREATMENT REGIME
text: Similar to Section 4.1, we compare the various methods in the ellipsoid framework . We observe that ES outperforms the ellipsoid method. Additionally, we compare the accuracy, i.e., how often does the policy derived from W match the expert’s policy, which is derived from W ∗. As IRL is not a supervised learning problem, we observe that while there is a correlation between the success in the task and the ability to act similarly to the experts policy - this correlation is not strict, in the sense that the agent is capable of finding near-optimal policies with a relatively high miss rate (accuracy of approximately 80%). For more intuition, see the proof of Lemma 1.

heading: B EXPERIMENTAL DETAILS
text: In this section, we describe the technical details of our experiments, including the hyper-parameters used. To solve MDPs, we use value iteration. Our implementation is based on a stopping condition with a tolerance threshold, τ , such that the algorithm stops if |Vt−Vt−1| < τ. In the driving simulation we used τ = 10−4 and in the sepsis treatment we use τ = 10−3.

heading: B.1 AUTONOMOUS DRIVING SIMULATION
text: The environment is modeled as a tabular MDP that consists of 1531 states. The speed is selected once, at the initial state, and is kept constant afterward. The other 1530 states are generated by 17 X-axis positions for the agent’s car, 3 available speed values, 3 lanes and 10 Y-axis positions in which car B may reside. During the simulation, the agent controls the steering direction of the car, moving left or right, i.e., two actions.
In these experiments, we define our mappings in a way that induces different behaviours for different contexts, making generalization a more challenging task. Specifically, for the linear setting we use W ∗ = ( −1 0.75 0.75
0.5 −1 1 0.75 1 −0.75
) , before normalization. For our nonlinear mapping, contexts with
||c||∞ > 0.55 are mapped to reward coefficients vector (1,−1,−0.05), otherwise they are mapped to (−0.01, 1,−1), which induce the feature expectations (9.75, 3.655, 5), (5.25, 5, 2.343) respectively. The decision regions for the nonlinear mapping are visualized in Appendix B.1. The contexts are sampled uniformly in the 2-dimensional simplex. We evaluate all algorithms on the same sequences of contexts, and average the results over 20 such sequences. The algorithms were modified to fit the Ellipsoid framework; instead of iterating over the whole data, the algorithms iterate over the given expert feature expectations, one at a time, until convergence, i.e. every timestep a new demonstration is presented.
Hyper-parameter selection and adjustments:
Ellipsoid Framework: For the linear model the algorithms maintained a 3× 3 matrix to estimate W ∗. Ellipsoid: By definition, the ellipsoid algorithm is hyper-parameter free and does not require tuning.
PSGD: The algorithm was executed with with the parameters: α0 = 0.3, αt = 0.9tαt−1, and iterated for 40 epochs. An outer decay on the step size was added for faster convergence, the initial α0
becomes 0.94 · α0 every time a demonstration is presented. The gradient, gt is normalized to be gt = gt
gt ||gt||∞ and the calculated step is taken if: cWt ( µ(π̂tc)−µ(π∗c ) ) > cWt+1 ( µ(π̂t+1c )−µ(π∗c ) ) ,
where π̂tc denotes the optimal policy for a context c according to Wt.
ES: The algorithm was executed with the parameters: σ = 10−3,m = 250, α = 0.1 with decay rate of 0.95, for 50 iterations which didn’t iterate randomly over one the contexts, but rather used the entire training set (all of the observed contexts and expert demonstrations up to the current time-step) for each step. The matrix was normalized according to || · ||2, and so was the step calculated by the ES algorithm, before it was multiplied by α and applied.
For the nonlinear setting, the model used for the ES method was a fully connected DNN, with layers of sizes 15, 10, 5, 3. The activation function used was the leaky ReLU function, with a parameter of α = 0.1. Note that we can’t normalize the parameters here as in the linear case; therefore an L2-normalization layer is added to the output. The same parameters were used as in the linear case, except with 120 iterations over the entire training set. They were originally optimized for this model and setting and worked as-is for the linear environment. As we aim to estimate the gradient, a small σ = 10−3 was used and performed best. The number of points,m = 250, was selected as fewer points produced noisy results. The step size, decay rate and the number of iterations were 0.1, 0.96, 120 respectively, and were selected in a way that produced fast yet accurate convergence of the loss. Note that here the steps were also normalized before application, and the normalization was applied per layer. For PSGD, a similar network was used. Specifically, it had layer sizes 14, 10, 6, 3, and the same leaky ReLU activation function was used in this network. In parallel to the normalization used for the ES model, here we used batch normalization and gradient clipping. The learning rate was set to 0.1 · 0.98t and 120 iterations were performed. For this result, as with the ES method, the batch consisted of all available training data.
Offline Framework: In the offline framework we compute the subgradients using expert trajectories of length 40, instead of the feature expectations. In this framework at every iteration we sample a mini batch of 10 contexts (from a finite set) and their corresponding trajectories (sampled from the expert policy and dynamics) then taking one descent step according to them. The generalization is measured over 80 holdout contexts that are referred as the test set, where the W that is used to calculate the feature expectations of the agent is fitted to the EW algorithm requirement to be in the dk− 1 simplex. The PSGD and the EW algorithms are configured as the theory specifies, where each descent step is calculated from the whole batch. The ES algorithm is applied with the parameters σ = 10−3,m = 500, α = 0.1 with decay rate 0.95, for every iteration.

heading: B.2 DYNAMIC TREATMENT REGIME
text: The environment we describe in 4.2 simulates a decision-making process for treating sepsis. Sepsis is a life-threatening severe infection, where the treatment applied to a sepsis patient is crucial for saving his life. To create a sepsis treating simulator, we leverage the MIMIC-III data set (Johnson et al., 2016). This data set includes data from hospital electronic databases, social security, and archives from critical care information systems, that had been acquired during routine hospital care. We follow the data processing steps that were taken in Jeter et al. (2019) to extract the relevant data in a form of normalized measurements of sepsis patients during their hospital admission and the treatments that were given to each patient. The measurements include dynamic measures, e.g., heart rate, blood pressure, weight, body temperature, blood analysis standard measures (glucose, albumin, platelets count, minerals, etc.), as well as static measures such as age, gender, re-admission (of the patient), and more.
The processed data from Jeter et al. (2019) consists of 5366 trajectories, each representing the sequential treatment provided by a clinician to a patient. At each time-step, the available information for each patient consists of 8 static measurements and 41 dynamic measurements. In addition, each trajectory contains the reported actions performed by the clinician (the number of fluids and vasopressors given to a patient at each time-step and binned to 25 different values), and there is a mortality signal which indicates whether the patient was alive 90 days after his hospital admission.
In order to create a tabular CMDP from the processed data, we separate the static measurements of each patient and keep them as the context. We cluster the dynamic measurements using K-means (MacQueen et al., 1967). Each cluster is considered a state and the coordinates of the cluster centroids are taken as its features φ(s). We construct the transition kernel between the clusters using the
empirical transitions in the data given the state and the performed actions. Two states are added to the MDP and the feature vector is extended by 1 element, corresponding to whether or not the patient died within the 90 days following hospital release. This added feature receives a value of 0 on all non-terminal states, a value of −0.5 for the state representing the patient’s death and 0.5 for the one representing survival. In addition, as the data is limited, not all state-action pairs are available. In order to ensure the agent does not attempt to perform such an action for which the outcome is unknown, we add an additional terminal state. At this state, all features are set to−1 to make it clearly distinguishable from all other states in the CMDP.
In our simulator, we used the same structure as the raw data, i.e., we used the same contexts prevalent in the data and the same initial state distribution. Each context is projected onto the simplex and the expert’s feature expectations for each context are attained by solving the CMDP. While we focus on a simulator, as it allows us to analyze the performance of the algorithms, our goal is to have a reward structure which is influenced by the data. Hence, we produce W ∗ by running the ellipsoid algorithm on trajectories obtained from the data. As done in the Autonomous Driving Simulation, the algorithms were modified to fit the Ellipsoid framework.
Hyper-parameter selection and adjustments:
Ellipsoid Framework: ES algorithm: the same method as in the autonomous driving is applied, with the parameters σ = 10−4,m = 1000, α = 0.25 with decay rate 0.95, for 80 iterations over the entire training set.
Offline Framework: In the offline framework for the linear setting, we compute the subgradients using expert trajectories of length 40, instead of the exact feature expectations. At every iteration, we sample a mini batch of 10 contexts (from a finite set) and their corresponding trajectories (sampled from the expert policy and dynamics) and perform a single descent step. Generalization is measured over a set of 300 holdout contexts, referred to as the test set, where W is fitted to be in the dk − 1 simplex. The PSGD and the EW algorithms are configured as specified by the theory, where each descent step is calculated from the entire batch. The ES algorithm is applied with the parameters σ = 10−4,m = 1000, α = 0.3 with decay rate 0.95, for every iteration.
For the nonlinear setting, the model used for the ES method was a fully connected DNN, with layers of sizes 24, 42, 42 which include a bias term. The activation function used was the leaky ReLU function, with a parameter of α = 0.1. In this setting we use trajectories of length 80, mini batch of 32 contexts. The PSGD algorithm uses the gradient computed by back-propagation of the loss function value on the DNN, where each descent step is calculated from the entire batch. The algorithm is applied with modified parameters; α = 0.1 with decay rate 0.98, for every iteration.

heading: C PROOFS FOR SECTION 3
text: Definition 1 (Bregman distance). Let ψ :W → R be strongly convex and continuously differential in the interior ofW . The Bregman distance is defined byDψ(x, y) = ψ(x)−ψ(y)−(x−y) ·∇ψ(y), where ψ is strongly convex with parameter σ. Definition 2 (Conjugate function). The conjugate of a function ψ(y), denoted by ψ∗(y) is maxx∈W {x · y − ψ(x)}.
Example: let ‖ · ‖ be a norm on Rn. The associated dual norm, denoted ‖ · ‖∗, is defined as ‖z‖∗ = sup{zᵀx | ‖x‖ ≤ 1}. The dual norm of ‖ · ‖2 is ‖ · ‖2, and the dual norm of ‖ · ‖1 is ‖ · ‖∞.
Before we begin with the proof of Lemma 2, we make the following observation. By definition, π̂c(W ) is the optimal policy w.r.t. cTW ; i.e., for any policy π we have that
cTW · µ(π̂c(W )) ≥ cTW · µ(π). (2)

heading: Proof of Lemma 1.
text: Fix W . For any context c, we have that µ(π̂c) is the optimal policy w.r.t. reward fW (c), thus, fW (c) · ( µ(π̂c(W )) − µ(π∗c ) ) ≥ 0. Therefore we get that Loss(W ) ≥ 0. For W ∗, we have that µ(π̂c(W )) = µ(π ∗ c ), thus Loss(W ∗) = 0.
For the second statement, note that Loss(W ) = 0 implies that ∀c, fW (c)· ( µ(π̂c(W ))−µ(π∗c ) ) = 0. This can happen in one of two cases. (1) µ(π̂c(W )) = µ(π∗c ), in this case π ∗ c , π̂c(W ) have the same feature expectations. Therefore, they are equivalent in terms of value. (2) µ(π̂c(W )) 6= µ(π∗c ), but fW (c) · ( µ(π̂c(W ))− µ(π∗c ) ) = 0. In this case, π∗c , π̂c(W ) have different feature expectations, but still achieve the same value w.r.t. reward fW (c). Since π̂c(W ) is an optimal policy w.r.t. this reward, so is π∗c .

heading: Proof of Lemma 2.
text: 1. We need to show that ∀W1,W2 ∈ W,∀λ ∈ [0, 1], we have that
Llin(λW1 + (1− λ)W2) ≤ λLlin(W1) + (1− λ)Llin(W2)
Llin(λW1 + (1− λ)W2) = Ec [ cT (λW1 + (1− λ)W2) · ( µ ( π̂c ( λW1 + (1− λ)W2 )) − µ(π∗c ) )] = λEc [ cTW1 · ( µ ( π̂c ( λW1 + (1− λ)W2 )) − µ(π∗c )
)] + (1− λ)Ec [ cTW2 · ( µ ( π̂c ( λW1 + (1− λ)W2 )) − µ(π∗c )
)] ≤ λEc [ cTW1 · ( µ ( π̂c ( W1 )) − µ(π∗c ) )] + (1− λ)Ec [ cTW2 · ( µ ( π̂c ( W2 )) − µ(π∗c )
)] = λLlin(W1) + (1− λ)Llin(W2),
where the inequality follows from Eq. (2).
2. Fix z ∈ W. We have that Llin(z) = Ec [ cT z · ( µ(π̂c(z))− µ(π∗c ) )] ≥ Ec [ cT z · ( µ(π̂c(W ))− µ(π∗c )
)] = Llin(W ) + (z −W ) · Ec [ c ( µ(π̂c(W ))− µ(π∗c ) )] ,
where the inequality follows from Eq. (2).
3. Recall that a bound on the dual norm of the subgradients implies Lipschitz continuity for convex functions. Thus it is enough to show that ∀W ∈ W, ‖g(W )‖p =
‖Ec [ c ( µ(π̂c(W ))− µ(π∗c ) )] ‖p ≤ L. Let p =∞, we have that
‖g(W )‖∞ = ∥∥∥Ecc (µ(π̂c(W ))− µ(π∗c ))∥∥∥∞
≤ Ec‖c ( µ(π̂c(W ))− µ(π∗c ) ) ‖∞ (Jensen inequality)
≤ Ec‖c‖∞ ‖µ(πi)− µ(πj)‖∞ ≤ 2
1− γ . (3)
where in Eq. (3) we used the fact that ∀π we have that ‖µ(π)‖∞ ≤ 1 1−γ , thus, for any πi, πj , ‖µ(πi)− µ(πj)‖∞ ≤ 2 1−γ Therefore, L = 21−γ w.r.t. ‖·‖∞. Since ‖·‖2 ≤ √ dk ‖·‖∞ we get that L = 2 √ dk 1−γ w.r.t. ‖·‖2 .

heading: D PROOFS & PSEUDO CODE FOR SECTION 3.3

heading: D.1 ELLIPSOID ALGORITHM FOR TRAJECTORIES
text: Algorithm 4 Batch ellipsoid algorithm for COIRL
Initialize: Θ0 ← B∞(0, 1) = {x ∈ Rd·k : ||x||∞ ≤ 1} Θ1 ← MVEE(Θ0) i← 0, Z̄ ← 0, Z̄∗ ← 0 for t = 1, 2, 3, ... do
ct is revealed, Let W t be the center of Θt Play episode using π̂t = arg maxπ V
π cTt Wt
Θt+1 ← Θt if a sub-optimal action a is played at state s then
Expert provides H-step trajectory (sE0 = s, s E 1 , ..., s E H). Let x̂ ∗,H i be the H-step sample of
the expert’s feature expectations for ξ′i = 1s: x̂ ∗,H i = ∑H h=0 γ
hφ(sEh ) Let xi be the agent’s feature expectations for ξ′i : Eξ′i,P,πt [ ∑∞ h=0 γ hφ(sh)] Denote zi = ct xi, ẑ∗,Hi = ct x̂ ∗,H i i← i+ 1, Z̄ ← Z̄ + zi, Z̄∗ ← Z̄∗ + ẑ∗,Hi if i = n then
Θt+1← MVEE ({ θ ∈ Θt : ( θ −W t )T · ( Z̄ ∗ n − Z̄ n ) ≥ 0 }) i← 0, Z̄ ← 0, Z̄∗ ← 0

heading: D.2 MVEE COMPUTATION
text: This computation is commonly found in optimization lecture notes and textbooks. First, we define an ellipsoid by {x : (x − c)Q−1(x − c) ≤ 1} for a vector c, the center of the ellipsoid, and an invertible matrix Q. Our first task is computing Θ1- the MVEE for the initial feasibility set Θ0 = B∞(0, 1) = {x ∈ Rd·k : ||x||∞ ≤ 1}. The result is of course a sphere around 0: c1 = 0, Q1 = dkI .
For the update Θt+1 ← MVEE ({ θ ∈ Θt : ( θ −W t )T · at ≥ 0 }) , we define ãt = −1√
aTt Qtat at
and calculate the new ellipsoid by ct+1 = ct − 1dk+1Qãt , Qt+1 = d2k2 d2k2−1 (Qt − 2 dk+1Qtãtã T t Qt).

heading: D.3 PROOF OF THEOREM 3
text: For simpler analysis, we define a ”flattening” operator, converting a matrix to a vector: Rd×k → Rd·k by W = [w1,1, . . . , w1,k, . . . , wd,1, . . . , wd,k]. We also define the operator to be the composition of the flattening operator and the outer product: u v = [u1v1, . . . , u1vk, . . . , udv1, . . . , udvk]. Therefore, the value of policy π for context c is given by V πcTW∗ = c
TW ∗µ(π) = W ∗T (c µ(π)), where ||W ∗||∞ ≤ 1, ||c µ(π)||1 ≤ k1−γ .
Lemma 3 (Boyd & Barratt (1991)). If B ⊆ RD is an ellipsoid with center w, and x ∈ RD\{0}, we define B+ = MVEE({θ ∈ B : (θ − w)Tx ≥ 0}), then: V ol(B
+) V ol(B) ≤ e − 1 2(D+1) .
Proof of Theorem 3. We prove the theorem by showing that the volume of the ellipsoids Θt for t = 1, 2, ... is bounded from below. In conjunction with Lemma 3, which claims there is a minimal rate of decay in the ellipsoid volume, this shows that the number of times the ellipsoid is updated is polynomially bounded. We begin by showing that W ∗ always remains in the ellipsoid. We note that in rounds where V π ∗
cTt W ∗ − V π̂tcTt W∗ > , we have W
∗T ( ct ( µ(π∗ct) − µ(π̂t) )) > . In addition, as the agent acts
optimally w.r.t. the reward rt = cTt Wt, we have that W T t ( ct ( µ(π∗ct)−µ(π̂t) )) ≤ 0 . Combining
these observations yield:
(W ∗ −W t) T · ( ct ( µ(π∗ct)− µ(π̂t) )) > > 0 . (4)
This shows that W ∗ is never disqualified when updating Θt . Since W ∗ ∈ Θ0 this implies that ∀t : W ∗ ∈ Θt. Now we show that not only W ∗ remains in the ellipsoid, but also a small ball surrounding it. If θ is disqualified by the algorithm: (θ −W t)T · ( ct ( µ(π∗ct) − µ(π̂t) )) < 0 .
Multiplying this inequality by -1 and adding it to (4) yields: (W ∗−θ)T · ( ct ( µ(π∗ct)−µ(π̂t) )) > .
We apply Hölder inequality to LHS: < LHS ≤ ||W ∗ − θ||∞ · || ( ct ( µ(π∗ct) − µ(π̂t) )) ||1 ≤
2k 1−γ ||W ∗ − θ||∞. Therefore for any disqualified θ: ||W ∗ − θ||∞ > (1−γ) 2k , thus B∞(W ∗, (1−γ) 2k ) is never disqualified. This implies that ∀t : vol(Θt) ≥ vol(Θ0 ∩ B∞(W ∗, (1−γ) 2k )) ≥ vol(B∞(W ∗, (1−γ) 4k )). Finally, let MT be the number of rounds by T in which V π∗ cTt W ∗ − V π̂t cTt W ∗ > . Using Lemma 3 we get that: MT 2(dk+1) ≤ log ( vol(Θ1) ) − log ( vol(ΘT+1) ) ≤
log ( vol ( MVEE(B∞(0, 1)) )) − log ( vol(B∞(0, (1−γ) 4k )) ) ≤ log ( vol ( MVEE(B2(0, √ dk)) )) −
log ( vol(B2(0, (1−γ) 4k )) ) ≤ log ( ( 4k √ dk (1−γ) ) dk ) ≤ dk log 4k √ dk
(1−γ) . Therefore MT ≤ 2dk(dk + 1) log 4k √ dk
(1−γ) = O(d 2k2 log( dk(1−γ) )) .

heading: D.4 PROOF OF THEOREM 4
text: Lemma 4 (Azuma’s inequality). For a martingale {Si}ni=0, if |Si − Si−1| ≤ b a.s. for i = 1, ..., n:
P ( |Sn − S0| > b √ 2n log( 2δ ) ) < δ
Proof of Theorem 4. We first note that we may assume that for any t: ||W ∗−Wt||∞ ≤ 2. If Wt 6∈ Θ0, we update the ellipsoid by Θt ← MVEE ({ θ ∈ Θt : ( θ −W t )T · ej ≶ 0 }) where ej is the
indicator vector of coordinate j in which Wt exceeds 1, and the inequality direction depends on the sign of (Wt)j . If Wt 6∈ Θ0 still, this process can be repeated for a finite number of steps until Wt ∈ Θ0, as the volume of the ellipsoid is bounded from below and each update reduces the volume (Lemma 3). Now we have Wt ∈ Θ0, implying ||W ∗ −Wt||∞ ≤ 2. As no points of Θ0 are removed this way, this does not affect the correctness of the proof. Similarly, we may assume ||W ∗t −Wt||∞ ≤ 2 as W ∗t ∈ Θ0. We denote Wt which remains constant for each update in the batch by W . We define t(i) the timesteps corresponding to the demonstrations in the batch for i = 1, ..., n. We define z∗,Hi to be the expected value of ẑ∗,Hi , and z ∗ i to be the outer product of ct(i) and the feature expectations of the expert policy for W ∗t(i), ct(i), ξ ′ t(i) . We also denote W ∗ t(i) by W ∗ i . We bound the following term from below, as in Theorem 3:
(W ∗ −W )T · ( Z̄ ∗
n − Z̄ n ) = 1 n n∑ i=1 (W ∗ −W )T · (ẑ∗,Hi − zi) =
1
n n∑ i=1 (W ∗ −W )T · (z∗i − zi) + 1 n n∑ i=1 (W ∗ −W )T · (z∗,Hi − z ∗ i )+
1
n n∑ i=1 (W ∗ −W )T · (ẑ∗,Hi − z ∗,H i ) =
1
n n∑ i=1
(W ∗i −W )T · (z∗i − zi)︸ ︷︷ ︸ (1)
+ 1
n n∑ i=1
(W ∗ −W ∗i )T · (z∗i − zi)︸ ︷︷ ︸ (2) +
1
n n∑ i=1 (W ∗ −W )T · (z∗,Hi − z ∗ i )︸ ︷︷ ︸
(3)
+ 1
n n∑ i=1 (W ∗ −W )T · (ẑ∗,Hi − z ∗,H i )︸ ︷︷ ︸
(4)
(1): Since the sub-optimality criterion implies a difference in value of at least for the initial distribution which assigns 1 to the state where the agent errs, we may use identical arguments to the previous proof. Therefore, the term is bounded from below by .
(2): By assumption ||W ∗ −W ∗i ||∞ ≤ (1−γ) 8k thus since ||(z ∗ i − zi)||1 ≤ 2k1−γ by Hölder’s inequality the term is bounded by 4 .
(3): We have ||x∗,Hi − x∗i ||1 ≤ kγH 1−γ from definitions, thus ||z ∗,H i − z∗i ||1 ≤ kγH
1−γ since c ∈ ∆d−1. As mentioned previously we may assume ||W ∗ −Wt||∞ ≤ 2, therefore by Hölder’s inequality the term is bounded by 4 due to our choice of H: γ H = (1 − (1 − γ))H ≤ ( (1 − (1 − γ)) 1 1−γ )(1−γ)H
=( (1− (1− γ)) 1 1−γ )log( 8k (1−γ) ) ≤ e− log( 8k
(1−γ) ) = (1−γ) 8k .
(4): The partial sums ∑N i=1(W ∗−W )T · (z∗,Hi − ẑ ∗,H i ) for N = 0, ..., n form a martingale sequence.
Note that ||z∗,Hi ||1 ≤ k1−γ , and ||ẑ ∗,H i ||1 ≤ k1−γ . Also, we have that ||W ∗ − Wt||∞ ≤ 2, thus, we can apply Azuma’s inequality (Lemma 4) with b = 4k(1−γ) and with our chosen n this yields:∑n
i=1(W ∗ −W )T · (z∗,Hi − ẑ ∗,H i ) ≤ n 4 with probability of at least 1−
δ
2dk(dk+1) log( 16k √ dk (1−γ) ) .
Thus (W ∗ − W )T · ( Z̄ ∗
n − Z̄ n ) > 4 and as in Theorem 3 this shows B∞(W ∗, (1−γ) 8k ) is never disqualified, and the number of updates is bounded by 2dk(dk + 1) log(16k √ dk
(1−γ) ), and multiplied by n this yields the upper bound on the number of rounds in which a sub-optimal action is chosen. By union-bound, the required bound for term (4) holds in all updates with probability of at least 1− δ.

heading: E CONTEXTUAL POLICIES
text: Consider the following problem. Given a complete specification of an MDP, and a hypothesis class H , for each state s assign a hypothesis ha : C → A such that the return is maximized. The following theorem shows that it is NP-complete to find such a policy. We will use the class of linear separators. We define the following contextual MDP training problem. We are given an MDP with only transitions and a sample of m contexts and their reward (namely for each context ci we specify the rewards for each state and action). Theorem 5. There is a reduction from training problem of the union of k hyperplanes to the contextual MDP training problem of k + 1 states.
Proof. Consider the following MDP, which has two parameters r0 and r1 which will define the rewards for each context. The MDP composed from a line of k internal states i, and one sink state. Each internal state i ∈ [1, k − 1] has two actions: action 1 leads to the sink state with reward r1 and the action 0 lead to the next internal state i+ 1 with reward 0. In internal state k action 1 leads to the sink state with reward r1 and action 0 leads to the sink state with reward r0. In the sink state there is a singe action that stays in the sink state and has reward 0.
The context isC = Rd. The hypotheses class includes all hyperplanes, each hyperplane is characterize by a weight vector w ∈ Rd, and if w>c ≥ 0 then the action is 1 and otherwise it is 0. Given a sample of size m to the training problem (ci, yi), we generate rewards for the contextual MDP training problem by having the context ci and specifying the parameters r0 and r1. Specifically, we will set r0 = 1− yi and r1 = yi. Given k hyperplanes w1, . . . , wk, and m example, assume that they make e errors. By using those k hyperplanes in the k internal nodes for each context ci we have a reward 1 iff the union of the k hyperplanes classify it correctly. This implies that we have reward m− e. This is since each example that is labeled correctly will get a reward of 1 and each incorrect example will get a reward of 0.
Given an assignment of k hyperplanes to the internal nodes, which gets a return of m− e, we output as a hypothesis the union of the k hyperplanes. Again, the number of errors on the sample is exactly e errors.

Review summaries:
This paper introduces a formulation for the contextual inverse reinforcement learning (COIRL) problem and proposed three algorithms for solving the proposed problem. Theoretical analysis of scalability and sample complexity are conducted for cases where both the feature function and the context-to-reward mapping function are linear. Experiments were conducted in both a simulated driving domain and a medical treatment domain to compare the three proposed algorithms empirically. Empirical results for using a deep network as the contextual mapping function is also provided.

This work focuses on the problem of 'contextual' inverse reinforcement learning, where the reward is a function of the current state of the MDP, and a set of context features, which remain constant within each episode.  The primary contribution of this work is the formulation of inverse reinforcement learning (for restricted spaces of context-dependent reward functions) as a convex optimization problem.  Based on this formulation, the paper describes several IRL algorithms based on approaches to solving convex and non-convex optimization problems, including variations of mirror descent, and evolution strategies (in principle allowing for the optimization of reward functions with arbitrary parametric representations).  The algorithms presented in this work all assume that computing an optimal policy for a specific reward function is a relatively inexpensive subroutine, which limits their applicability to domains where such planning is straightforward.  Experimental results are presented for a simple highway driving domain, as well as a simulated patient treatment domain constructed from real-world clinical data.