Bayesian Nonparametric Feature Construction for Inverse Reinforcement Learning

1 Introduction

Inverse reinforcement learning (IRL) aims to recover the expert's underlying reward function from her demonstrations and the environment model [Russell, 1998]. IRL is applied to increasingly various research areas, such as robotics [Argall et al., 2009], computer animation [Lee and Popovi, 2010], preference learning [Erkin et al., 2010; Ziebart et al., 2008b] and cognitive science [Baker et al., 2009].

In the last decade, a number of studies on IRL algorithms have appeared in the literature [Ng and Russell, 2000; Abbeel and Ng, 2004; Ratliff et al., 2006; Syed et al., 2008; Ziebart et al., 2008a]. Most of them, especially when dealing with large state spaces, assume that the reward function is a linear function of some pre-defined features. The problem then reduces to learning the unknown weights of the linear function, but the result is highly dependent on the selection of features. Hence, it is desirable to automatically construct the features that compactly describe the structure of the reward function [Abbeel and Ng, 2004].

FIRL [Levine et al., 2010] is one of the few IRL algorithms that construct the features for the reward function. It is assumed that the domain expert supplies the set of all potentially relevant features as atomic features, and the algorithm constructs the logical conjunctions of these features for the reward function. GPIRL [Levine et al., 2011] is another IRL algorithm that represents the reward function as a nonlinear function of the atomic features using the Gaussian process (GP). The reward features are implicitly learned as the hyperparameters of the GP kernel.

In this paper, we propose a Bayesian nonparametric approach to constructing the reward function features in IRL. We extend Bayesian IRL [Ramachandran and Amir, 2007] by defining a prior on the composite features, which are defined to be the logical conjunctions of the atomic features. Since the number of composite features is not known a priori, we define the prior using the Indian buffet process (IBP) to infer the features and the number of features. Our approach has a number of advantages: First, it learns an explicit representation of the reward function features in the form of logical formulas, which are readily interpretable. Second, it can robustly learn from a noisy behaviour data since it uses a probabilistic model of the data. Third, it can incorporate the domain knowledge on the reward features and their weights into the prior distribution since it is a Bayesian framework.

2 Preliminaries

We use the following notations throughout the paper: x is a column vector with elements xi and x-i = x\xi. X is a matrix with elements Xi,j, X:,j is the j-th column of X, X-(i,j) =X\Xi,j, andX-i,j =X:,j\Xi,j.

We denote a discrete-state Markov decision process (MDP) [Puterman, 1994] as a tuple M = 〈S,A, {T a }a∈A, r, γ〉 where S is the set of states, A is the set of actions, T a is the state transition probability such that T a s,s′ = P(s′ |s, a), r is the (state-dependent) reward function, and γ ∈ [0, 1) is the discount factor.

A policy is defined as a mapping π : S → A. The value of the policy π is defined as the expected discounted cumulative rewards of executing π, which satisfies the equation vπ = r + γT π vπ where T π s,s′ = P(s′ |s, π(s)). Similarly, the Q-function is defined by the equation Qπ:,a = r+γT avπ . Given an MDPM, an optimal policy π∗ maximizes the value function for all the states. The value of an optimal policy should satisfy the Bellman optimality equation: v∗ = maxa∈A(r + γT av∗ ).

The IRL refers to the problem of inferring the reward function of an expert given an MDP\r and the expert's behaviour data D. We assume that D = {τ1, . . . , τN} is generated by executing an optimal policy with respect to the unknown reward vector r, where the n-th trajectory τn is an H-step sequence of state-action pairs, i.e., τn = {(sn,1, an,1), . . . , (sn,H, an,H)}.

2.1 BIRL: Bayesian Framework for IRL

Ramachandran and Amir [2007] proposed a Bayesian framework for IRL by modeling the compatibility of the reward with the behaviour data as the likelihood and the preference on the rewards as the prior. The likelihood is defined as an independent exponential, or softmax, distribution:

where η is the parameter representing the confidence of actions being optimal, and Q∗ (r) is the optimal Q-function computed using r. Assuming that the reward entries are independently distributed, the prior on the reward function is defined as P(r) = ∏s∈S P(rs). Various distributions such as the uniform, normal, or Laplace distribution can be used as the prior depending on the domain knowledge on the reward function. The reward function is then inferred from the posterior distribution, using Markov Chain Monte Carlo (MCMC) for estimating posterior mean [Ramachandran and Amir, 2007] or gradient ascent for estimating the maximuma-posteriori (MAP) [Choi and Kim, 2011].

2.2 Indian Buffet Process

The Indian buffet process (IBP) [Ghahramani et al., 2007] defines a distribution over binary matrices with infinitely many columns. This process can be explained by a culinary metaphor: The customers enter an Indian buffet one after another and choose dishes from an infinite number of dishes. The first customer chooses dishes whose number is distributed by Poisson(α). Thereafter, the i-th customer chooses the dish j with a probability of ζj/i where ζj is the number of customers who have previously chosen the dish j. In addition, the i-th customer chooses new dishes whose number is distributed by Poisson(α/i). Let Z be an M ×K binary matrix where M and K is the total number of customers and dishes, whose elements Zij = 1 indicates that the i-th customer chooses the dish j. By taking the limit K → ∞, we obtain the probability distribution induced by the IBP:

where ζk = ∑Mm=1 Zmk. The parameter α controls the sparsity of the matrix.

3 BNP-FIRL: Bayesian Nonparametric Feature Construction for IRL

We assume that the reward vector is linearly parameterized so that r = Φw, where Φ = [φ1, . . . , φK] is a reward feature matrix whose k-th column φk is an |S|-dimensional binary vector representing the k-th reward feature, and w is an Kdimensional weight vector.

As we have mentioned in the introduction, most of the IRL algorithms assume that the reward features are already given, and just compute the weights [Ng and Russell, 2000; Abbeel and Ng, 2004; Ratliff et al., 2006; Syed et al., 2008; Ziebart et al., 2008a]. In contrast, the goal of our work is to learn the reward features as well as to compute the weights. A small assumption made here is that a domain expert provides a set of (binary) atomic features, which includes all the relevant indicators that may possibly influence determining the rewards in some unknown way. Since the true reward function is hardly a linear function of the atomic features, we should learn the reward features as (non-linear but hopefully not too complex) functions of atomic features, so that we can still represent the reward function as a linear function of reward features. We refer to the reward features as composite features in order to distinguish them from atomic features.

In the following sections, we describe our Bayesian nonparametric approach to the feature construction in IRL, where the composite features are constructed by taking the logical conjunctions of the atomic features. We first define a nonparametric prior over the logical conjunctions of the atomic features using the IBP. We then derive the posterior over the composite features and the reward weights by extending BIRL. Finally, we present an MCMC algorithm for the posterior inference.

3.1 Prior for Feature Construction

Given a matrix of atomic features, defined as Ψ = [ψ1, . . . ,ψM] where the m-th column ψm is an |S|- dimensional binary vector representing the m-th atomic feature, we define the constructed composite feature matrix as Φ = f (x, Z, U ; Ψ) with:

where x is an M-dimensional binary vector indicating the participation of the m-th atomic feature in any of the composite features, Z is an M × K binary matrix indicating the participation of the m-th atomic feature in the k-th composite feature, and U is also an M × K binary matrix indicating whether the m-th atomic feature is negated in the formula for the k-th composite feature so that ˜ ψ(k) m = ψm if Um,k = 1 and ˜ ψ(k) m = ¬ψm otherwise.

We define the prior on Z, U, and x as follows: first, we start with Z that indicates which atomic features are used for each composite feature. Each composite feature generally consists of more than one exchangeable atomic features. In addition, we do not know a priori the number K of columns in Z, which is the number of the composite features to be constructed. Thus, analogous to the latent feature modeling, we use the IBP as the prior on Z:

Note that, by using the IBP, we inherently assume that K is unbounded even though it is at most 3M . The rationale behind our assumption is that Z can have duplicated columns as well as infinitely many zero columns. Second, we use the Bernoulli distribution for every entry of the binary matrix U, so that the atomic features are negated with p = 0.5:

Finally, we define the prior on the binary vector x. Since xm = 0 implies that the atomic feature ψm is not used at all for any composite feature, we use the prior that favors x being sparse. Specifically, we use the Bernoulli distribution with p = κ where κ is Beta distributed with β = [β1, β2]:

This is analogous to controlling the row-wise sparsity in the IBP proposed by Rai and Daum´e III [2008].

By combining Eqns (3), (4), and (5), the prior is defined as

3.2 Posterior Inference

BNP-FIRL extends BIRL by using the prior defined in the above and treating the behaviour D = {τ1, . . . , τN} as being drawn from the generative process as follows:

Fig 1 shows the graphical model of BNP-FIRL. Note that, as in BIRL, the reward weights wk's are assumed to be independently distributed so that P (w) = ∏Kk=1 P (wk).

The posterior over the composite features for the reward function and the associated weights is then formulated as

where Θ = {η, α, β,Ψ}, P(D|Φ,w, η) is the BIRL likelihood defined in Eqn (1) with r = Φw, P(w) is the prior on the weights, and P(Φ|α, β,Ψ) is the prior on the composite features defined in Eqn (6). Note that α controls the number of composite features to be constructed and β controls the total number of the unique atomic features used in

Algorithm 1 MCMC algorithm for BNP-FIRL

Initialize x, Z, U , w. for t = 1 to T do for m = 1 to M do Sample xm according to Eqn (8). end for for m = 1 to M do for k = 1 to K do Sample Um,k according to Eqn (9). end for end for for m = 1 to M do for k = 1 to K do if ∑i̸=m Zi,k > 0 then Sample Zm,k according to Eqn (10). end if end for Propose ξ = 〈K+, Z+, U +, w+〉: K+ ∼ Poisson(α/M). Accept ξ with probability min{1, ρZ}. end for for k = 1 to K do Propose w′ ∼ N(wk, λ). Accept w′ with probability min{1, ρw}. end for end for

the construction. In other words, α and β are the parameters that control the column-wise and the row-wise sparsity of the composite feature matrix Φ, respectively.

We infer the composite features and the reward weights from the posterior using an MCMC algorithm described as follows (Algorithm 1): we first update xm by sampling from the probability distribution conditioned on the rest of the random variables,

Note in the above that, instead of drawing the Bernoulli distribution parameter κ from Beta(β = [β1, β2]) and then drawing xm, we have collapsed κ for an efficient inference.

Next, in order to sample Um,k which indicates whether we should negate the atomic feature ψm in the composite feature φk, we sample it from the likelihood

since the prior on Um,k is uniformly random in {0, 1}.

We then sample Zm,k, which indicates whether the atomic feature ψm appears in the composite feature φk. For the sake of exposition, we rephrase the IBP culinary metaphor into our context, the atomic features being the customers and the composite features being the dishes. The first atomic feature chooses Poisson(α) composite features. The m-th atomic feature chooses the k-th composite feature (already chosen by preceding atomic features) with probability ∑m-1 i=1 Zi,k/m and additionally chooses Poisson(α)/m new composite features, where Zi,k = 1 indicates that the i-th atomic feature ψi have chosen the k-th composite feature φk.

This metaphor and the exchangeability of the atomic features leads to the update method for Zm,: [Rai and Daum´e III, 2008] as follows: first, for the composite features that are already chosen (i.e., ∑i̸=m Zi,k > 0), we update Zm,k according to the conditional distribution

Second, for choosing the new composite features, we sample Zm,k using the Metropolis-Hastings (MH) update: we first sample ξ = 〈K+,Z+,U +,w+〉 where K+ ∼ Poisson(α/m), Z+ is the m×K+ binary matrix whose entries in the m-th row are set to 1 and others are set to 0, and the m × K+ binary matrix U + and K+-dimensional vector w+ are drawn from the corresponding priors. We then accept ξ with probability min{1, ρZ} where

In the above equation, [X,Y ] denotes horizontal concatenation, and [x; y] denotes vertical concatenation.

Finally, we sample the weights wk, again using the MH update. We first sample w′ ∼ N(wk, λ) and then accept it with probability min{1, ρw} where

with wnew formed by taking wnew k = w′ and wnew -k =w-k. The posterior mean is commonly used for inferring the reward function since it is known to minimize the square error LSE(r, ˆ r) = ||r - ˆ r||2 [Ramachandran and Amir, 2007]. We thus estimate the posterior mean of the reward function. On the other hand, when we show the learned composite features, we choose the sample with the maximum posterior, i.e., 〈 ˆΦMAP, ˆwMAP〉 = argmax〈Φ(t) ,w(t) 〉 P(Φ(t) ,w(t) |D,Θ) and ˆ rMAP = ˆΦMAP ˆwMAP. This is because the sample mean is ill-defined for ˆΦ(t) 's with different dimensions.

4 Experimental Results

In this section, we show the performance of BNP-FIRL in a number of problem domains and compare it to FIRL [Levine et al., 2010] and GPIRL [Levine et al., 2011]. 1 The performance was evaluated using the expected value difference (EVD) 1 |S| ||v∗ (r) - vπ(r′ ) (r)||1 where r is the expert's reward function (i.e., the ground truth), r′ is the learned reward function, and π(r′ ) is the optimal policy computed using r′ . The EVD thus can be seen as a measurement of the loss in the optimality incurred by using the policy from the learned reward function instead of the expert's reward function. In order to evaluate how well the composite features are learned, we computed the EVD on the same problem instance (i.e., the training problem instance) where the expert's behaviour data was generated, and on additional random problem instances (i.e., transfer problem instances) as well.

1 code available at http://graphics.stanford.edu/projects/gpirl

4.1 Objectworld Domain

The first set of experiments was performed on the objectworld domain [Levine et al., 2011], where the agent can move north, south, east, west, or stay in the current location in an N ×N grid. Each action has a failure probability of 0.3 which makes the agent move in a random direction. There are a number of colored objects randomly placed in the grid, each of them with one of the C ≥ 2 colors.

We prepared a total of CN atomic features of the form "d(s, ci) < j" indicating that the distance between the agent's location s and the nearest color i object is less than j, where i ∈ {1, . . . ,C} and j ∈ {1, . . . ,N}. The true reward function was set to

We generated 10 random training problem instances by sampling the locations of objects and their colors, and gathered trajectories of length 20. In order to measure how well the learned reward function generalizes to novel yet similar problem instances, we measured EVD on additional 10 transfer problem instances for each training problem instance, generated in the same way we prepared training problem instances.

The top row in Fig 2 shows the EVD performances when the trajectories are generated by the optimal policy. BNPFIRL outperformed FIRL in the experiments. This is mainly due to the fact that BNP-FIRL is allowed to have multiple composite features for each state, whereas FIRL can only have one because of the way the algorithm constructs the composite features. Tbl 1 shows the 4 reward features and the corresponding weights found by BNP-FIRL. In comparison, FIRL produced a total of 21 features. On the other hand, BNP-FIRL was on par in performance with GPIRL. Nonetheless, one of the advantages of using BNP-FIRL is that the learned features are explicitly represented and thus readily interpretable.

The bottom row in Fig 2 shows the EVD performances when the trajectories are generated from an ϵ-greedy policy. It is interesting to observe that BNP-FIRL outperforms GPIRL and FIRL. We conjecture that this is due to the likelihood in BNP-FIRL being particularly robust to noisy behaviour data.

4.2 Simulated-highway Domain

The second set of experiments was done on the simulatedhighway domain [Levine et al., 2011], where the agent drives a vehicle by moving one lane left or right at speeds {1, 2, 3, 4}

on a three-lane highway. The actions that change the lane or speed fail with probability 0.3. There are other vehicles driving at speed 1, each of them being either a civilian or a police vehicle and either a car or a motorcycle (a total combination of 4 categories). We prepared an expert that prefers to drive as fast as possible but avoids driving at speeds 3 or 4 within a distance of 2 from a police vehicle. We also prepared 3 types of atomic features: The first indicate the current speed (4 features) and the second indicate the current lane (3 features). The third indicate the distance to the nearest car from each category, of the form "d(s, c) ≤ j" where c ∈ {police, civilian}×{car, motorcycle} and j ∈ {0, . . . , 5} (24 features). We generated the trajectory data by executing the optimal policy for 200 time steps.

Fig 3 shows that BNP-FIRL again performs better than FIRL. On the other hand, BNP-FIRL performs slightly worse than GPIRL although, as shown in Fig 4, the learned reward functions from the two algorithms were very similar to the true one. On the other hand, FIRL was not able to capture good features for the reward function and as a result, the learned reward function is very different from the true one.

4.3 Taxi Driver Behaviour Prediction

The final set of experiments was done on learning the taxi drivers' preference using the GPS trace data collected in San Francisco [Piorkowski et al., 2009]. We modeled the road network as a graph, which consists of 20531 road segments obtained from the OPENSTREETMAP. 2 As in Ziebart et al. [2008b], we assumed that the taxi drivers try to reach a destination in a trip by taking the road segments according to their preferences. Note that each trip typically has a different destination. We used the goal-oriented MDP Mg to define the environment model, where the states are road segments, and the state transitions correspond to taking one of the road segments at intersections. The rewards were assumed to be negative everywhere except at the destination, which was represented as an absorbing goal state g with a zero reward.

The composite features φ(s,a),k and the weights wk were assumed to be independent of the destination. In other words, given the destination g, the reward was calculated as rg = F(g) · w where wk < 0, F(g,a),k(g) = 0 and F(s,a),k(g) = φ(s,a),k for all s ∈ S\{g}. This change leads to a slight modification to the likelihood so that P(D|Φ,w, η) = ∏n P(τn|rgn , η;Mgn ) where gn is the destination in the trajectory τn. Note that, although we have different MDPs Mgn for each trajectory due to different destination, they all share the same reward features Φ and weights w.

We prepared state-dependent atomic features representing the properties of the road segments, such as the type, the speed limit, the number of lanes, and one-way. 3 We also prepared action-dependent atomic features representing the angle of the turn. 4

Among the traces of 500 taxis in the original data, we selected 500 trips from 10 taxis, visualized in Fig 5. The accumulated distance of the trips was 1137 miles. We then segmented the GPS traces and mapped them to the road segments using a preprocessing algorithm described in Lou et al. [2009].

2 http://www.openstreetmap.org

3 type ∈ {highway, primary street, secondary street, living street}, speed limit ∈ {below 20 mph, 20-30 mph, 30-40 mph, above 40 mph}, # lanes ∈ {1, 2, 3+}. These features were obtained using the OPENSTREETMAP

4 turn angle ∈ {hard left, soft left, straight, soft right, hard right, u-turn}.

Figure 5: GPS traces of 500 trips from 10 taxis collected in San Francisco.

We compared BNP-FIRL to two baseline methods and two IRL algorithms in the literature. As for the first baseline, we used the shortest path to the destination. As for the second baseline, we used the optimal policy from an MDP with heuristically set reward weight w′ . A simple heuristic would be counting the feature visitations in the trace data, and set the reward weights accordingly. Since we enforce the reward weights to be negative, we used the normalized counts of feature visitations offset by the maximum count value. As for the two IRL algorithms, we chose the maximum entropy IRL (MaxEntIRL) [Ziebart et al., 2008a] and the gradient ascent BIRL (MAP-BIRL) [Choi and Kim, 2011]. Since both algorithms depend on a pre-defined set of features, we supplied the set of all atomic features. We do not report the results from FIRL and GPIRL since it was not straightforward to modify them to appropriately handle multiple MDPs sharing the same reward features and weights.

Tbl 2 shows the average prediction accuracy and their standard errors. The prediction accuracy was measured in terms of the turn and the route predictions. For each trip, we computed the path from the origin to the destination using the algorithms. The correct actions taken at the intersections were counted in the turn prediction accuracy, and the ratio of the total distance taking the correct road segments was calculated for the route prediction accuracy. All the results were obtained by the 5-fold cross validation. Tbl 3 presents the weights and composite features learned by BNP-FIRL, which are fairly intuitive. For example, making hard left turns (φ1) is avoided with the highest penalty. Making turns on a highway is also highly undesirable except making a soft right turn to take an exit ramp (φ2). It was also found that taxi drivers prefer to take the primary street (φ9) or the road with 20-30 mph limit (φ8), rather than the highway (φ5) or the road with speed limit higher than 40 mph (φ3). This was because the trips were generally short, the average distance being 2.27 miles.

5 Conclusion

We presented BNP-FIRL, a Bayesian nonparametric approach to constructing reward features in IRL. We defined the reward features as the logical conjunctions of the atomic features provided by the domain expert, and formulated a nonparametric prior over the reward features using the IBP. We derived an MCMC algorithm to carry out the posterior inference on the reward features and the corresponding reward weights.

We showed that BNP-FIRL outperforms or performs on par with prior feature-learning IRL algorithms through experiments on synthetic domains. In comparison to FIRL, BNPFIRL produces more succinct sets of features with richer representations and learns better reward functions. In comparison to GPIRL, BNP-FIRL produces the set of features that are explicitly represented and readily interpretable. We also presented reward feature learning results on a real GPS trace data collected from taxi drivers, predicting their behaviour with higher accuracy than prior IRL algorithms.

Acknowledgments

This work was supported by National Research Foundation of Korea (Grant# 2012-007881), the Defense Acquisition Program Administration and Agency for Defense Development of Korea (Contract# UD080042AD), and the IT R&D program of MKE/KEIT (Contract# 10041678).