Title: BLACK-BOX OFF-POLICY ESTIMATION FOR INFINITE-HORIZON REINFORCEMENT LEARNING

Abstract: Off-policy estimation for long-horizon problems is important in many real-life applications such as healthcare and robotics, where high-fidelity simulators may not be available and on-policy evaluation is expensive or impossible. Recently, Liu et al. (2018) proposed an approach that avoids the curse of horizon suffered by typical importance-sampling-based methods. While showing promising results, this approach is limited in practice as it requires data be drawn from the stationary distribution of a known behavior policy. In this work, we propose a novel approach that eliminates such limitations. In particular, we formulate the problem as solving for the fixed point of a certain operator. Using tools from Reproducing Kernel Hilbert Spaces (RKHSs), we develop a new estimator that computes importance ratios of stationary distributions, without knowledge of how the off-policy data are collected. We analyze its asymptotic consistency and finite-sample generalization. Experiments on benchmarks verify the effectiveness of our approach.
Consider a Markov decision process (MDP) (Puterman, 1994) M = S, A, P, R, p 0 , γ , where S and A are the state and action spaces, P is the transition probability function, R is the reward function,

1. INTRODUCTION:

As reinforcement learning (RL) is increasingly applied to crucial real-life problems like robotics, recommendation and conversation systems, off-policy estimation becomes even more critical. The task here is to estimate the average long-term reward of a target policy, given historical data collected by (possibly unknown) behavior policies. Since the reward and next state depend on what action the policy chooses, simply averaging rewards in off-policy data does not estimate the target policy's long-term reward. Instead, proper correction must be made to remove the bias in data distribution.

One approach is to build a simulator that mimics the reward and next-state transitions in the real world, and then evaluate the target policy against the simulator (e.g., Fonteneau et al., 2013; Ie et al., 2019) . While the idea is natural, building a high-fidelity simulator could be extensively challenging in numerous domains, such as those that involve human interactions. Another approach is to use propensity scores as importance weights, so that we could use the weighted average of rewards in off-policy data as a suitable estimate of the average reward of the target policy. The latter approach is more robust, as it does not require modeling assumptions about the real world's dynamics. It often finds more success in short-horizon problems like contextual bandits, but its variance often grows exponentially in the horizon, a phenomenon known as "the curse of horizon" (Liu et al., 2018) .

To address this challenge, Liu et al. (2018) proposed to solve an optimization problem of a minimax nature, whose solution directly estimates the desired propensity score of states under the stationary distribution, avoiding an explicit dependence on horizon. While their method is shown to give more accurate predictions than previous algorithms, it is limited in several important ways:

• The method requires that data be collected by a known behavior policy. In practice, however, such data are often collected over an extended period of time by multiple, unknown behavior policies. For example, observational healthcare data typically contain patient records, whose treatments were provided by different doctors in multiple hospitals, each following potentially different procedures that are not always possible to specify explicitly. • The method requires that the off-policy data reach the stationary distribution of the behavior policy. In reality, it may take a very long time for a trajectory to reach the stationary distribution, which may be impractical due to various reasons like costs and missing data.

In this paper, we introduce a novel approach for the off-policy estimation problem that overcome these drawbacks. The main contributions of our work are three-fold:

• We formulate the off-policy estimation problem into one of solving for the fixed point of an operator. Different from the related, and similar, Bellman operator that goes forward in time, this operator is backward in time.

• We develop a new algorithm, which does not have the aforementioned limitations of Liu et al. (2018) , and analyze its generalization bounds. Specifically, the algorithm does not require that the off-policy data come from the stationary distribution, or that the behavior policy be known. • We empirically demonstrate the effectiveness of our method on several classic control benchmarks. In particular, we show that, unlike Liu et al. (2018) , our method is effective even if the off-policy data has not reached the stationary distribution.

In the next section, we give a brief overview of recent and related works. We then move to describing the problem setting that we have used in the course of the paper and our off-policy estimation approach. Finally, we present several experimental results to show the effectiveness of our method.

Notation. In the following, we use ∆(X ) to denote the set of distributions over a set X . The 2 norm of vector x is x . Given a real-valued function f defined on some set X , let f 2 := X f (x) 2 dx. Finally, we denote by [n] the set {1, 2, . . . , n}, and 1{A} the indicator function.

2. RELATED WORKS:

Our work focuses on estimating a scalar (average long-term reward) that summarizes the quality of a policy and has extensive applications in practice. This is different from value function or policy learning from off-policy data (e.g., Precup et al., 2001; Maei et al., 2010; Sutton et al., 2016; Munos et al., 2016; Metelli et al., 2018) , where the major goal is to ensure stability and convergence. Yet, these two problems share numerous core techniques, such as importance reweighting and doubly robustness. Off-policy estimation and evaluation can also be used as a component for policy optimization (e.g., Jiang & Li, 2016; Gelada & Bellemare, 2019; Liu et al., 2019; Zhang et al., 2019) .

Importance reweighting, or inverse propensity scoring, has been used for off-policy RL (e.g., Precup et al., 2001; Murphy et al., 2001; Li et al., 2015; Munos et al., 2016; Hanna et al., 2017; Xie et al., 2019) . Its accuracy can be improved by various techniques (Jiang & Li, 2016; Thomas & Brunskill, 2016; Guo et al., 2017; Farajtabar et al., 2018) . However, these methods typically have a variance that grows exponentially with the horizon, limiting their application to mostly short-horizon problems like contextual bandits (Dudík et al., 2011; Bottou et al., 2013) .

There have been recent efforts to avoid the exponential blow-up of variance in basic inverse propensity scoring. A few authors explored the alternative to estimate the propensity score of a state's stationary distribution (Liu et al., 2018; Gelada & Bellemare, 2019) , when behavior policies are known. Later, Nachum et al. (2019) extended this idea to situations with unknown behavior policies. However, their approach only works for the discounted reward criterion. In contrast, our work considers the more general and challenging undiscounted criterion. In the next section, we briefly mention the setting under which we study this problem and then present our black-box off-policy estimator.

Our black-box estimator is inspired by previous work for black-box importance sampling (Liu & Lee, 2017) . Interestingly, the authors show that it is beneficial to estimate propensity scores from data without using knowledge of the behavior distribution (called proposal distribution in that paper), even if it is available; see also Henmi et al. (2007) for related arguments. Similar benefits may exist for our black-box off-policy estimator developed here, although a systematic study is outside the scope of this paper. p 0 ∈ ∆(S) is the initial state distribution, and γ ∈ [0, 1] is the discount factor. A policy π maps states to a distribution over actions: π : S → ∆(A), and π(a|s) is the probability of choosing action a in state s by policy π. With a fixed π, a trajectory τ = (s 0 , a 0 , r 0 , s 1 , a 1 , r 1 , . . .) is generated as follows:1 

s 0 ∼ p 0 (•), a t ∼ π(•|s t ), r t = R(s t , a t ), s t+1 ∼ P (•|s t , a t ), ∀t ≥ 0 .

Given a target policy π, we consider two reward criteria, undiscounted (γ = 1) and discounted (γ < 1), where E π [•] indicates the trajectory τ is controlled by policy π:

EQUATION

EQUATION

In the above, d π is the stationary distribution over S × A, which exists and is unique under certain assumptions (Levin & Peres, 2017) .

The γ < 1 case can be reduced to the undiscounted case of γ = 1, but not vice versa. Indeed, one can show that the discounted reward in equation 2 can be interpreted as the stationary distribution of an induced Markov process, whose transition function is a mixture of P and the initial-state distribution p 0 . We refer interested readers to Appendix A for more details. Accordingly, in the following and without the loss of generality, we will merely focus on the more general undiscounted criterion in equation 1, and suppress the unnecessary dependency on p 0 and γ.

In the off-policy estimation problem, we are interested in estimating ρ π for a given target policy π. However, instead of having access to on-policy trajectories generated by π, we have a set of n transitions collected by some unknown (i.e., "black-box" or behavior-agnostic (Nachum et al., 2019) ) behavior mechanism π BEH :

D := {(s i , a i , r i , s i )} 1≤i≤n .

Therefore, the goal of off-policy estimation is to estimate ρ π based on D, for a given target policy π.

The setting we described above is quite general, covering a number of situations. For example, π BEH might be a single policy and D might consist of one or multiple trajectories collected by π BEH . In this special case, s i = s i+1 for 1 < i < n; this is the off-policy RL scenario widely studied (e.g., Precup et al., 2001; Sutton et al., 2016; Munos et al., 2016; Liu et al., 2018; Gelada & Bellemare, 2019) . Furthermore, if π BEH = π, we recover the on-policy setting. On the other hand, π BEH and D can consist of multiple policies and their corresponding trajectories. In this situation, unlike the single policy case s i and s i+1 might originate from two distinct policies. In general, one can consider π BEH as a distribution over S × A where (s i , a i ) in D are sampled from. Having introduced the general setting of the problem, we will describe our estimation approach in the next section.

4. BLACK-BOX ESTIMATION:

Our estimator is based on the following operator defined on functions over S × A. For discrete state-action spaces, given any

EQUATION

While we will develop the rest of the paper using the discrete version above for simplicity, the continuous version can be similarly obtained without affecting the estimator and results:

EQUATION

where P is now interpreted as the transition kernel.

We should note that B π is indeed different from the Bellman operator (Puterman, 1994) ; although they have some similarities. In particular, given some state-action pair (s, a), the Bellman operator is defined using next state s of (s, a), while B π is defined using previous state-actions (ξ, α) that transition to s. It is in this sense that B π is backward (in time). Furthermore, as we will show later, d has the interpretation of a distribution over S × A. Therefore, B π describes how visitation flows from (ξ, α) to (s, a) and hence, we call it the backward flow operator. Note that similar forms of B π have appeared in the literature, usually used to encode constraints in a dual linear program for an MDP (e.g., Wang et al., 2007; Wang, 2017; Dai et al., 2018) . However, the application of B π for the off-policy estimation problem as considered here appears new to the best of our knowledge.

An important property of B π is that, under certain assumptions, the stationary distribution d π is the unique fixed point that lies in ∆(S × A) (Levin & Peres, 2017) :

EQUATION

)

This property is the key element we use to derive our estimator as we describe in the following.

4.1. BLACK-BOX ESTIMATOR:

In most cases, off-policy estimation involves a weighted average of observed rewards r i in D. We therefore aim to directly estimate these (non-negative) weights which we denote by w

= {w i } ∈ ∆([n]); that is, w i ≥ 0 for i ∈ [n] and n i=1 w i = 1.

Note that the normalization of w may be ensured by techniques such as self-normalized importance sampling (Liu, 2001) . Once such a w is obtained, the estimated reward is given by:

EQUATION

Effectively, any w ∈ ∆([n]) defines an empirical distribution which we denote by d w over S × A:

EQUATION

)

Equation 6 is equivalent to ρπ = E (s,a)∼dw [r] . Comparing it to equation 1, we naturally want to optimize w so that d w is close to d π . Therefore, inspired by the fixed-point property of d π in equation 5, the problem naturally becomes one of minimizing the discrepancy between d w and B π d w .

In practice, w is often represented in a parametric way:

EQUATION

where W (.) is a parametric model, such as neural networks, with parameters ω ∈ Ω. We have now reached the following optimization problem:

EQUATION

where D(• •) is some discrepancy function between distributions. In practice, B is unknown, and must be approximated by samples in the dataset D:

Bπ d w (s, a) := π(a|s) n i=1 w i 1{s i = s} .

Clearly, Bπ d w is a valid distribution over S × A induced by w and D, and the black-box estimator solves for w by minimizing D(d w Bπ d w ).

4.2. BLACK-BOX ESTIMATOR WITH MMD:

There are different choices for D(• •) in equation 9, and multiple approaches to solve it (e.g., Nguyen et al., 2010; Dai et al., 2017) . Here, we describe one such algorithm based on Maximum Mean Discrepancy (MMD) (Muandet et al., 2017a) . For simplicity, the discussion in this subsection assumes S × A is finite, but the extension to continuous S × A is immediate.

Let k(•, •) be a positive definite kernel function defined on (S × A) 2 . Given two real-valued functions, f and g, defined on S × A we define the following bilinear functional

EQUATION

Clearly, we have k [f ; f ] ≥ 0 for any f due to the positive definiteness of k. In addition, k is called

strictly integrally positive definite if k [f ; f ] = 0 implies f 2 = 0.

Let H be the reproducing kernel Hilbert space (RKHS) associated with the kernel function k. This is the unique Hilbert space that includes functions that can be expressed as a sum of countably many terms:

f (•) = i u i k((s i , a i ), •)

, where {u i } ⊂ R, and {(s i , a i )} ⊆ S × A. The space H is equipped with an inner product defined as follow: given f, g ∈ H such that

f = i u i k((s i , a i ), •) and g = i v i k((s i , a i ), •), the inner product is f, g H := i,j u i v j k(s i , a i ), (s j , a j )), which

induces the norm defined by f H := i,j u i u j k((s i , a i ), (s j , a j )).

Given H, the maximum mean discrepancy between two distributions, µ 1 and µ 2 , is defined by

D k (µ 1 µ 2 ) := sup f ∈H {E µ1 [f ] -E µ2 [f ], s.t. f H ≤ 1} .

Here, f may be considered as a discriminator, playing a similar role as the discriminator network in generative adversarial networks (Goodfellow et al., 2014) , to measure the difference between µ 1 and µ 2 . A useful property of MMD is that it admits a closed-form expression (Gretton et al., 2012) : w i w j k((s i , a i ), (s j , a j ))

D k (µ 1 µ 2 ) = k [µ 1 -µ 2 ; µ 1 -µ 2 ] = k [µ 1 ; µ 1 ] -2k [µ 1 ; µ 2 ] + k [µ 2 ; µ 2 ] , where k [•; •] is

K (0) i,j k d w ; Bπ d w = i,j w i w j a π(a |s j )k((s i , a i ), (s j , a )) K (1) i,j k Bπ d w ; Bπ d w = i,j w i w j a i ,a j π(a i |s i )π(a j |s j )k((s i , a i ), (s j , a j )) K (2) i,j

.

Defining K i,j := K (0) i,j -2K (1) i,j + K (2)

i,j , we can express the objective as a function of ω (since {w i } depends on ω; see equation 8):

(ω) = i,j

EQUATION

Remark 4.1. Mini-batch training is an effective approach to solve large-scale problems. However, the objective (ω) is not in a form that is ready for mini-batch training, as w i requires normalization (equation 8) that involves all data in D. Instead, we may equivalently minimize L(ω) := log (ω), which can be turned into a form that allow mini-batch training, using a trick that is also useful in other machine learning contexts (e.g., Jean et al., 2015) . See Appendix D for more details.

Algorithm 1 in Appendix E summarizes our estimator. We next present theoretical analysis of our approach. We show the consistency of our result and provide a sample complexity bound. 

D k (d || B π d) = D k(d || d π ) ,

where k(x, x) is a positive definite kernel, defined by

k(x, x) = E π [k(x, x) -k(x, x ) -k(x , x) + k(x , x ) | (x, x)],

where the expectation is with respect to 

x ∼ d π (• | x) and x ∼ d π (• | x),

D k (d ŵ || B π d ŵ) -D k (d w * || B π d w * ) ≤ 16r max R n (W) + 16r 2 max + r 2 max 8 log(1/δ) √ n ,

where R n (W) denotes the expected Rademacher complexity of W with data size n, and

r max := max( W ∞ , sup x k(x, x)), with W ∞ := sup{ W ∞ : W ∈ W}. This suggests a gener- alization error of O(1/ √ n) if R n (W) = O(1/ √ n)

, which is typical for parametric families of functions.

5. EXPERIMENTS:

In this section, we present experiments to compare the performance of our proposed method with other baselines on the off-policy evaluation problem. In general and for each experiment, we use a behavior policy π BEH to generate trajectories of length T BEH . We then use these generated samples from a behavior policy to estimate the expected reward of a given target policy π. To compare our approach with other baselines, we use the root mean squared error (RMSE) with respect to the average long-term reward of the target policy π. The latter is estimated using a trajectory of length T TAR 1. In particular, we compare our proposed black-box approach with the following baselines:

• naive averaging baseline in which we simply estimate the expected reward of a target policy by averaging the rewards over the trajectories generated by the behavior policy. • model-based baseline where we use the kernel regression technique with a standard Gaussian RBF kernel. We set the bandwidth of the kernel to the median (or 25 th or 75 th percentiles) of the pairwise euclidean distances between the observed data points. • inverse propensity score (IPS) baseline introduced by Liu et al. (2018) .

We will first use a simple MDP from Thomas & Brunskill (2016) to highlight the IPS drawback we previously mentioned in Section 1. We then move to classical control benchmarks. The RMSE of different methods as we change the length of horizon w.r.t the target policy reward. IPS depends on the horizon length while our method is independent of the horizon length.

5.1. TOY EXAMPLE:

The ModelWin domain first introduced in Thomas & Brunskill ( 2016) is a fully observable MDP with three states and two actions as denoted in Figure 1 (a). The agent always begins in s 1 and should choose between two actions a 1 and a 2 . If the agent chooses a 1 , then with probability of p and 1 -p it makes a transition to s 2 and s 3 and receives a reward of r = 1 and r = -1, respectively. On the other hand, if the agent chooses a 2 , then with probability of p and 1 -p it makes a transition to s 3 with the reward of r = -1 and s 2 with the reward of r = 1, respectively. Once the agent is in either s 2 or s 3 , it goes back to the s 1 in the next step without any reward. In our experiments, we set p = 0.4.

We define the behavior and target policies as the following. In the target policy, once the agent is in s 1 , it chooses a 1 and a 2 with the probability of 0.9 and 0.1, respectively. On the other hand and for the behavior policy, once the agent is in s 1 , it chooses a 1 and a 2 with the probability of 0.7 and 0.3, respectively. We calculate the average on-policy reward from samples based on running a trajectory of length T TAR = 50, 000 collected by the target policy. We estimate this on-policy reward using trajectories of length T BEH ∈ {4, 8, 16, 32, 64, 128} collected by the behavior policy. In each case, we set the number of trajectories such that the total number of transitions (i.e., T BEH times the number of trajectories) is kept constant. For example, for T BEH = 4 and T BEH = 8 we use 50,000 and 25,000 trajectories, respectively. Since the problem has finitely many state-actions, we use the tabular method and hence, equation 11 turns into a quadratic programming. We then report the result of each setting based on 10 Monte-Carlo samples.

As we can see in Figure 1 (b), the naive averaging method performs poorly consistently and independent of the length of trajectories collected by the behavior policies. On the other hand, IPS performs poorly when the collected trajectories have short-horizon and gets better as the horizon length of trajectories get larger. This is expected for IPS -as mentioned in Section 1, it requires data be drawn from the stationary distribution. In contrast, as shown in Figure 1 (b), our black-box approach performance is independent of the horizon length, and substantially better.

5.2. CLASSIC CONTROL:

We now focus on four classic control problems. We begin by briefly describing each problem and then compare the performance of our method with other approaches on these problems. Note that for these problems are episodic, we convert them into infinite-horizon problems by resetting the state to a random start state once the episode terminates.

Pendulum. In this environment, the goal is to control a pendulum in a vertical position. State variables are the pole angle θ and velocity θ. The action a is the torque in the set {-2, -1, 0, 1, 2} applied to the base. We set the reward function to -(θ 2 + 0.1 θ2 + 0.001a 2 ).

Mountain Car. For this problem, the goal is to drive up the car to top of a hill. Mountain Car has a state space of R 2 (the position and speed of the car) and three possible actions (negative, positive, or zero acceleration). We set the reward to +100 when the car reaches the goal and -1 otherwise. Cartpole. The goal here is to prevent an attached pole to a cart from falling by changing the cart's velocity. Cartpole has a state space of R 4 (cart position, velocity, pole angle and velocity) and two possible actions (moving left or right). Reward is -100 when the pole falls and +1 otherwise.

Acrobot. In this problem, our goal is to swing a 2-link pendulum above the base. Acrobot has a state space of R 6 (sin(.) and cos(.) of both angles and angular velocities) and three possible actions (applying +1, 0 or -1 torque on the joint). Reward is +100 for reaching the goal and -1 otherwise.

For each environment, we train a near-optimal policy π + using the Neural Fitted Q Iteration algorithm (Riedmiller, 2005) . We then set the behavior and target policies as π BEH = α 1 π + + (1 -α 1 )π -and π = α 2 π + + (1 -α 2 )π -, where π -denotes a random policy, and 0 ≤ α 1 , α 2 ≤ 1 are two constant values making the behavior policy distinct from the target policy. In our experiments, we set α 1 = 0.7 and α 2 = 0.9. In order to calculate the on-policy reward, we use a single trajectory collected by π with T TAR = 50, 000. For off-policy data, we use multiple trajectories collected by π BEH with T BEH = 200. In all the cases, we use a 3-layer (having 30, 20, and 10 hidden neurons) feed-forward neural network with the sigmoid activation function as our parametric model in equation 8. For each setting, we report results based on 20 Monte-Carlo samples.

Figure 2 shows the log of RMSE w.r.t. the target policy reward as we change the number of trajectories collected by the behavior policy. We should note that all methods except the naive averaging method have hyperparameters to be tuned. For each method, the optimal set of parameters might depend on the number of trajectories (i.e., size of the training data). However, in order to avoid excessive tuning and to show how much each method is robust to a change in the number of trajectories, we only tune different methods based on 50 trajectories and use the same set of parameters for other settings. As we can see, the naive averaging performance is almost independent of the number of trajectories. Our method outperforms other approaches on three environments and it is only the Acrobot in which IPS performs comparably to our black-box approach. In order to have a robust evaluation against outliers, we have plotted the median and error bars at 25 th and 75 th percentiles in Figure 3 . If we compare the Figures 2 and 3 , we notice that the trend of results is almost the same in both. In Appendix G, we have studied how changing α 1 of the behavior policy affects the final RMSE.

6. CONCLUSIONS:

In this paper, we presented a novel approach for solving the off-policy estimation problem in the long-horizon setting. Our method formulates the problem as solving for the fixed point of a "backward flow" operator. We showed that unlike previous works, our approach does not require the knowledge of the behavior policy or stationary off-policy data. We presented experimental results to show the effectiveness of our approach compared to previous baselines. In the future, we plan to use structural domain knowledge to improve the estimator and consider a random time horizon in episodic RL.

A REDUCTION FROM DISCOUNTED TO UNDISCOUNTED REWARD:

The same reduction is used in Liu et al. (2018) . For completeness, we give the derivation details here, for the case of finite state/actions. The derivation can be extended to general state-action spaces, with proper adjustments in notation.

Let τ = (s 0 , a 0 , r 0 , s 1 , a 1 , . . .) be a trajectory generated by π, and d t ∈ ∆(S × A) be the distribution of (s t , a t ). Clearly, d 0 (s, a) = p 0 (s)π(a|s)

d t+1 (s, a) = ξ,α

d t (ξ, α)P (s|ξ, α)π(a|s), ∀t > 0 .

Using a matrix form, the recursion above can be written equivalently as d t+1 = P T π d t , where P π is given by P π (s, a|ξ, α) = P (s|ξ, α)π(a|s) . The discounted reward of policy π is

ρ π,γ = (1 -γ)E π ∞ t=0 γ t r t = E (s,a)∼dπ,γ [R(s, a)] , with d π,γ := (1 -γ) d 0 + γd 1 + γ 2 d 2 + • • • .

Multiplying both sides of above by γP T π , we have

γP T π d π,γ = (1 -γ) γP T π d 0 + γ 2 P T π d 1 + γ 3 P T π d 2 + • • • = (1 -γ) γd 1 + γ 2 d 2 + γ 3 d 3 + • • • = d π,γ -(1 -γ)d 0 . Therefore, d π,γ = γP T π d π,γ + (1 -γ)d 0 = (γP π + (1 -γ)d 0 1 T ) T d π,γ .

Accordingly, d π,γ is the fixed point of an induced transition matrix given by P π,λ := γP π + (1γ)d 0 1 T . This completes the reduction from the discounted to the undiscounted case.

B PROOF OF THEOREM 4.1 Note that 

D k (d B π d) = k [d -B π d; d -B π d] .

D k (d B π d) = k [(d -B π d); (d -B π d)] = k [(d -B π d -(d π -B π d π )); (d -B π d -(d π -B π d π ))] = k [(δ w -B π δ w ); (δ w -B π δ w )] .

Recalling that B π d(x) = x0 P π (x|x 0 )d(x 0 ), we have

D k (d B π d) = k [(δ w -B π δ w ); (δ w -B π δ w )] = x,x k(x, x)(δ w (x) -B π δ w (x))(δ w (x) -B π δ w (x)) = x,x k(x, x) δ w (x) - x0 P π (x|x 0 )δ w (x 0 ) δ w (x) - x0 P π (x|x 0 )δ w (x 0 ) .

Define the adjoint operator of B π ,

P π f (x) := x P π (x |x)f (x ).

Denote by P x π the operator applied on k(x, x) in terms of variable x, that is, P x π k(x, x) := x P π (x |x)k(x , x). This gives

D k (d B π d) = x,x k(x, x ) δ w (x) - x0 P π (x|x 0 )δ w (x 0 ) δ w (x) - x0 P π (x|x 0 )δ w (x 0 ) = x,x δ w (x) k(x, x) -P x π k(x, x) -P x π k(x, x) + P x π P x π k(x, x) δ(x) = x,x δ w (x) kπ (x, x)δ w (x) ,

completing the proof.

C PROOF OF THEOREM 4.2

First, note that the error can be decomposed as follows: 

D k (d ŵ || B π d ŵ) ≤ D k (d ŵ || Bπ d ŵ) + D k ( Bπ d ŵ || B π d ŵ) ≤ D k (d w * || Bπ d w * ) + D k ( Bπ d ŵ || B π d ŵ) ≤ D k (d w * || B π d w * ) + D k ( Bπ d w * || B π d w * ) + D k ( Bπ d ŵ || B π d ŵ) ≤ D k (d w * || B π d w * ) + 2 sup w∈W D k ( Bπ d w || B π d w ) = D k (d w * || B π d w * ) + 2Z,

It remains to calculate:

Z := sup w∈W D k ( Bπ d w B π d w )). Recall from the definition of D k that D k ( Bπ d w B π d w )) = sup f ∈B k E x∼ Bπdw [f (x)] -E x∼Bπdw [f (x)] . Recall that {(x i , x i )} n i=1 is a set of transitions with x i ∼ d π (• | x i )

, where d π denotes the transition probability under policy π. Following the definitions of Bπ d w and B π d w , we have

Bπ d w (x) = n i=1 w(x i )1{x = x i } , B π d w (x) = n i=1 w(x i )d π (x | x i ) . Therefore, Z = sup w∈W D k ( Bπ d w B π d w )) = sup w∈W sup f ∈B k E x∼ Bπdw [f (x)] -E x∼Bπdw [f (x)] = sup w∈W,f ∈B k 1 n n i=1 w(x i ) f (x i ) -E x i ∼dπ [f (x i )|x i ] ,

where x i ∼ d π (• | x i ) is an independent copy of x i that follows the same distribution. Note that x i is introduced only for the sample complexity analysis. Note that Z is a random variable, and E[Z] denotes its expectation w.r.t. random data {x i , x i } n i=1 . We assume different (x i , x i ) are independent with each other. First, by McDiarmid inequality, we have

P (Z ≥ E[Z] + ) ≤ exp - n 2 2 W 2 ∞ B k 2 ∞

. This is because when changing each data point (x i , x i ), the maximum change on Z is at most

2 W ∞ B k ∞ /n. Therefore, we have Z ≤ E[Z] + 2 log(1/δ) W 2 ∞ B k 2 ∞ /n with proba- bility at least 1 -δ.

Accordingly, we now just need to bound E[Z]:

E[Z] = E X sup w∈W,f ∈B k 1 n n i=1 w(x i ) (f (x i ) -E X [f (x i )|x i ]) ≤ E X, X sup w∈W,f ∈B k 1 n n i=1 w(x i )(f (x i ) -f (x i ) = E X, X,σ sup w∈W,f ∈B k 1 n n i=1 σ i w(x i )(f (x i ) -f (x i ) (because {σ i } are i.i.d. Rademacher random variables) ≤ 2E sup w∈W,f ∈B k 1 n n i=1 σ i w(x i )f (x i ) = 2R n (W ⊗ B k ), where W ⊗ B k = {f (x)g(x ) : f ∈ W, g ∈ B k }.

By Lemma C.1 below, we have

E[Z] ≤ 2R n (W ⊗ B k ) ≤ 4 ( W ∞ + B k ∞ ) (R n (W) + R n (B k )) .

Combining the bounds, we have with probability 1 -δ,

2Z ≤ 4R n (W ⊗ B k ) + 8 log(1/δ) W 2 ∞ B k 2 ∞ n ≤ 8 ( W ∞ + B k ∞ ) (R n (W) + R n (B k )) + 8 log(1/δ) W 2 ∞ B k 2 ∞ n .

Plugging Lemma C.2, we have

2Z ≤ 8 ( W ∞ + r k ) R n (W) + 8r k W ∞ + r k + W ∞ log(1/δ)/8 √ n .

Assume r max = max( W ∞ , r k ). We have

2Z ≤ 16 r max R n (W) + 16r 2 max + r 2 max 8 log(1/δ) √ n . Lemma C.1. Denote by W ∞ = sup{ f ∞ : f ∈ W} the super norm of a function set W. Then, R n (W ⊗ B k ) ≤ 2 ( W ∞ + B k ∞ ) (R n (W) + R n (B k )) .

Proof. Note that f (x)g(x ) = 1 4 (f (x) + g(x )) 2 -1 4 (f (x) -g(x )) 2 .

Note that x 2 is 2(

W ∞ + B k ∞ )-Lipschitz on interval [-W ∞ -B k ∞ , W ∞ + B k ∞ ].

Applying Lemma C.1 of Liu & Wang (2018) , we have

R n (W ⊗ B k ) ≤ 2( W ∞ + B k ∞ )(R n (W ⊕ B k ),

where W ⊕ B k = {f (x) + g(x ) : f ∈ W, g ∈ B k }, and

R n (W ⊕ B k ) = E vz [ sup f ∈W,g∈B k i z i (f (x i ) + g(x i ))] ≤ E vz [ sup f ∈W i z i f (x i )] + E vz [ sup g∈B k i z i g(x i )] = R n (W) + R n (B k ) .

Remark The same result holds when w is defined as a function of the whole transition pair (x, x ), that is, w i = w(x i , x i ).

Lemma C.2. Let H be the RKHS with a positive definite kernel k(x, x ) on domain X . Let B k = {f ∈ H k : f H k ≤ 1} be the unit ball of H k , and r k = sup x∈X k(x, x ). Then,

B k ∞ ≤ r k , R n (B k ) ≤ r k √ n .

Proof. These are standard results in RKHS theory, and we give a proof for completeness. For B k ∞ , we just note that for any f ∈ B k and x ∈ X ,

f (x) = f, k(x, •) H k ≤ f H k k(x, •) H k ≤ k(x, •) H k = k(x, x) ≤ r k .

The Rademacher complexity can be bounded as follows: 

R n (B k ) = E X,σ sup f ∈B k 1 n i σ i f (x i ) ≤ E X,σ   sup f ∈B k f, 1 n n i=1 σ i k(x i , •) H k   = E X,σ   1 n n i=1 σ i k(x i , •) H k   ≤ E X,σ   1 n n i=1 σ i k(x i , •) 2 H k   1/2 = E X,σ   1 n 2 n i,j=1 σ i σ j k(x i , x j )   1/2 = E X 1 n 2 n i=1 k(x i , x i ) 1/2 ≤ r k √ n .

G SENSITIVITY:

Finally, in Figure 5 we measure how robust our approach is to changing the behavior policy compared to other methods. In particular, we vary α 1 that corresponds to the behavior policy to measure how the RMSE is affected. While α 2 is fixed to 0.9, in each experiment we choose α 1 from {0.7, 0.5, 0.3, 0.1}.

For each experiment, we use data from 50 trajectories (with T BEH = 200) collected by the behavior policy and report results based on 20 Monte-Carlo samples. According to Figure 5 , as α 1 diverges more from α 2 , the performance of all the methods degrade while our method is the least affected. It is worth mentioning that for the Mountain Car problem and α 1 = 0.1, the behavior policy is close to a random policy and hence the car has not been able to drive up to top of the hill. This means that all the methods have constantly received a reward of -1 during all the steps and hence the estimated on-policy reward has been -1 for all the methods as well. Therefore, the RMSE of all four methods are equal in this case.

:

For simplicity in exposition, we assume rewards are deterministic. However, everything in this work generalizes directly to the case of stochastic rewards.

Figures and Tables:

defined in equation 10, and we used the bilinear property k [•; •]. Interested readers are referred to surveys (e.g. Berlinet & Thomas-Agnan, 2011; Muandet et al., 2017b) for more background on RKHS and MMD. Applying MMD to our objective, we have D k (d w Bπ d w ) = k [d w ; d w ] -2k d w ; Bπ d w + k Bπ d w ; Bπ d w . In the above, both d w and Bπ d w are simply probability mass functions on a finite subset of S × A, consisting of state-actions encountered in D. It follows immediately from equation 10 that k [d w ; d w ] = i,j

THEORETICAL ANALYSISConsistency. The following theorem shows that the exact minimizer of equation 9 coincides with the fixed point of B π , and the objective function measures the norm of the estimation error in an induced RKHS. To simplify exposition, we assume x = (s, a) and x = (s , a ) a successive actionstate pair following x: x ∼ d π (• | x), where d π (x | x) is the transition probability from x to x , that is, d π (x | x) = P (s | s, a)π(a |s ). Similarly, we denote by (x, x ) = ((s, ā), (s , ā )) an independent copy of (x, x ). Theorem 4.1. Suppose k is strictly integrally positive definite, and d π is the unique fixed point of B π in equation 5. Then, for any d ∈ ∆(S × A), D k (d || B π d) = 0 ⇐⇒ d = d π . Furthermore, D k (d || B π d) equals an MMD between d and d π , with a transformed kernel:

with x and x drawn independently. Generalization. We next give a sample complexity analysis. In practice, the estimated weight ŵ is based on minimizing the empirical loss D k (d w || Bπ d w ), where B π is replaced by the empirical approximation Bπ . The following theorem compares the empirical weights ŵ with the oracle weight w * obtained by minimizing the expected loss D k (d w || B π d w ), with the exact transition operator B π . Theorem 4.2. Assume the weight function is decided by w i = W (s i , a i ; ω)/n. Denote by W = { W (•; ω) : ω ∈ Ω} the model class of W (•; ω). Assume ŵ is the minimizer of the empirical loss D k (d w || Bπ d w ) and w * the minimizer of expected loss D k (d w || B π d w ). Assume {x i } n i=1 are i.i.d. samples. Then, with probability 1 -δ we have

Figure 1: (a) ModelWin MDP from Thomas & Brunskill (2016). (b) The RMSE of different methods as we change the length of horizon w.r.t the target policy reward. IPS depends on the horizon length while our method is independent of the horizon length.

Figure 2: The RMSE of different methods w.r.t the target policy reward as we change the number of trajectories. Our black-box approach outperforms other methods on three problems.

Figure 3: The median and error bars at 25 th and 75 th percentiles of different methods w.r.t the target policy reward as we change the number of trajectories. The trend of results is similar to Figure 2.

Following the definition of the strictly integrally positive definite kernels, we have that D k (d B π d) = 0 implies d -B π d = 0, which in turn implies d = d π by the uniqueness assumption of d π . We have thus proved the first claim. For the second claim, define δ w = d -d π . Since d π -B π d π = 0, we have

where we define Z := D k ( Bπ d w || B π d w ). Therefore, we just need to bound Z. Denote by B k := {f : f ∈ H, f H ≤ 1} the unit ball of RKHS. Define B k := sup f ∈B k and R n (B k ) the expected Rademacher complexity of B k of data size n. From classical RKHS theory (see Lemma C.2 below), we know that B k ∞ ≤ r k and R n (B k ) ≤ r k / √ n, where r k := sup x∈X k(x, x ).

Figure 4: Bias and standard deviation of different methods for the ModelWin MDP (Fig. Our method has a smaller bias but larger variance compared to other algorithms.

Figure 5: The RMSE of different methods w.r.t the target policy reward as we change the behavior policy. Our method outperform other approaches on different behavior policies.

Peer review (iclr2020):
This paper proposes a new algorithm to the off-policy evaluation problem in reinforcement learning, based on stationary state visitation. The proposed algorithm does not need to know the behavior policy probabilities, which provide a broader application scenario compared with previous work. It is based on directly solving the fixed point problem of the resulting target stationary distribution, instead of solving a fixed point problem about the stationary ratios. In general, this paper is looking at an interesting problem with a lot of increasing focus these days. I think this paper contributes great idea and good theoretical analysis. I like the idea of matching the resulting target distribution instead of minimizing a loss over the ratio. However, several unclear places in the current version hurt the presentation of results. I would like to see them get improved and will increase my score accordingly if so. Detailed comments: 1) The algorithm part could be presented more clearly. So far I did not see where the empirical operator hat{B} is formally defined. The word *empirical* is also confusing to me in *B is approximated by empirical data* because B is not an expectation, but an *integral* which has no *empirical* opposite of it. For the equations on top of page 5, shouldn*t they be k[, ] about empirical operator instead of the expected operator since the RHS is already in sample-based form? 2) Related with the last one, B has an integral. To approximate the integral, we only have one sample from the transition probability actually, and the sample state is not uniformly distributed. It needs some explanation of why that would not cause a problem to approximate the integral. 3) The current loss function is invariant to the scale of w at all. Since the w is normalized, this is not a problem for the resulting estimator, ideally. However, that can be a numerical issue for float numbers. It*s possible that the output from function approximator w goes to 0 or infty. Both cases can lead to NaN of the output/function approximator update eventually. I*d like to hear if the author has met this problem in the experiment or not, and how can that be fixed. 4) I have to point out, as just a slight con of this paper, the technique used in this paper is not that much different from Liu et al 2018. Since it minimizes a loss function which is a supremum over an RKHS, and the resulting empirical loss also has a similar form. It*s nice to see the author provide some details of making it work with mini-batch. These details are important for function approximator as NN. Minor point: - On page 12 the equations after *We have by the definition of D_k*, I did not follow the second step of the equations. Suggestion: This paper study the similar settings (behavior-agnostic OPE), using similar method (on the stationary distribution) came out several months ago: https://arxiv.org/pdf/1906.04733.pdf. I knew it*s unfair to ask the author to compare it with a very recent prior/parallel work. However since they are in such a similar case, and they have code available, is it possible to directly compare with the result from their code? https://github.com/google-research/google-research/tree/master/dual_dice ======= After rebuttal ======= The author*s feedback clarified some of my concerns in the initial review. After reading the author*s feedback and other reviews, I think this paper has enough show contribution to the related work. Some of my previous concerns (point 2 and 3) seems true for many related works in this area in general. I partly agree that it is not very fair to ask this paper to fix them. The updated version also presents the theory section in a more clear way. So I*d like to raise my score. I*ve no problem with acceptance, but I won*t fully heart argue for acceptance.