Computer Science Mathematics
Computer Science-2021.gz:143146
Computer Science
8990c62f00fc31c3dabdba049eb0a1d888c3ee40
Non-asymptotic Confidence Intervals of Off-policy Evaluation: Primal and Dual Bounds

Off-policy evaluation (OPE) is the task of estimating the expected reward of a given policy based on offline data previously collected under different policies. Therefore, OPE is a key step in applying reinforcement learning to real-world domains such as medical treatment, where interactive data collection is expensive or even unsafe. As the observed data tends to be noisy and limited, it is essential to provide rigorous uncertainty quantification, not just a point estimation, when applying OPE to make high stakes decisions. This work considers the problem of constructing non-asymptotic confidence intervals in infinite-horizon off-policy evaluation, which remains a challenging open question. We develop a practical algorithm through a primal-dual optimization-based approach, which leverages the kernel Bellman loss (KBL) of Feng et al.(2019) and a new martingale concentration inequality of KBL applicable to time-dependent data with unknown mixing conditions. Our algorithm makes minimum assumptions on the data and the function class of the Q-function, and works for the behavior-agnostic settings where the data is collected under a mix of arbitrary unknown behavior policies. We present empirical results that clearly demonstrate the advantages of our approach over existing methods.


Introduction
Off-policy evaluation (OPE) seeks to estimate the expected reward of a target policy in reinforcement learning (RL) from observational data collected under different behavior policies (e.g., Murphy et al., 2001;Fonteneau et al., 2013;Jiang and Li, 2016;Liu et al., 2018a). OPE plays a central role in applying reinforcement learning (RL) with only observational data and has found important applications in areas such as medical treatments, autonomous driving, where interactive "on-policy" data is expensive or even infeasible to collect.
As the observed data is often noisy and limited, applying OPE to assist high-stakes decision-making entails a rigorous, non-asymptotic approach to quantify the uncertainty in the estimation. In this way, the decision makers are informed of the degree of the risk and can avoid the dangerous case of being overconfident in making costly and/or irreversible decisions. This work addresses the uncertainty quantification in OPE by constructing non-asymptotically correct confidence intervals that contains the true expected reward of the target policy with a probability no smaller than a user-specified confidence level.
To date, Off-policy evaluation per se has remained a key technical challenge in the literature (e.g., Precup, 2000;Thomas and Brunskill, 2016;Jiang and Li, 2016;Liu et al., 2018a), let alone gaining rigorous confidence estimation of it. This is especially true when 1) the underlying RL problem is long or infinite horizon (a.k.a. the curse of horizon (Liu et al., 2018a)); and 2) the data is collected under arbitrary and unknown algorithms (a.k.a. behavior-agnostic). In such cases, the collected data can exhibit complex and unknown time-dependency structure, which makes constructing rigorous non-asymptotic confidence bounds particularly challenging. Traditionally, the main approach of deriving non-asymptotic confidence bounds in OPE is to combine importance sampling (IS) with certain concentration inequality (e.g., Thomas et al., 2015a). However, the classical IS-based methods tend to degenerate for long or infinite horizon problems (Liu et al., 2018a. Moreover, they cannot be applied to the behavior-agnostic settings and also fail to handle the complicated time-dependency structure inside individual trajectories. In fact, the IS-based methods only work for short-horizon problems with a large number of independently collected trajectories drawn under known policies. In this work, we provide a practical algorithm with theoretical guarantee for Behavior-agnostic, Off-policy, Infinite-horizon, Non-asymptotic, Confidence intervals based on arbitrarily Dependent data (BONDIC). This method is motivated by a recently proposed optimization-based approach to estimating OPE confidence bounds (Feng et al., 2020), which leverages a tail bound of kernel Bellman loss (Feng et al., 2019). Our approach achieves a new dual bound that is both an order-of-magnitude tighter and more computationally efficient than that of Feng et al. (2020). These improvements are based on two pillars: 1) a new primal-dual perspective on the non-asymptotic confidence bounds of infinite-horizon OPE; and 2) a new martingale concentration inequality on the kernel Bellman loss that applies to behavior-agnostic off-policy data with arbitrary time-dependency between transition pairs. Our simple and practical method can provide reliable and informative bounds on a variety of benchmarks.
Outline The rest of the paper is organized as follows. We introduce the problem setting in Section 2 and discuss related works in Section 3. In Section 4, we give an overview on two dual approaches to infinite-horizon OPE that are closely related to our method (but do not consider non-asymptotic error bounds). We then present a martingale concentration bound for kernel Bellman loss in Section 5, which is used in our main approach described in Section 6. Finally, we perform empirical studies in Section 7. The proofs and additional discussions can be found in Appendix.

Background, Data Assumption, Problem Setting
Consider an agent acting in an unknown environment. At each time step t, the agent observes the current state s t in a state space S, takes an action a t ∼ π(· | s t ) in an action space A according to a given policy π; then, the agent receives a reward r t and the state transits to s t = s t+1 , following an unknown transition/reward distribution (r t , s t+1 ) ∼ P(· | s t , a t ). Assume the initial state s 0 is drawn from an known initial distribution D 0 . Let γ ∈ (0, 1) be a discount factor. In this setting, the expected reward of π is defined as J * := J π,P := lim which is the expected total discounted rewards when we execute π starting from D 0 for H steps. The focus of this work is the infinite-horizon case with H → +∞, although our method can also be applied to a finite-horizon, time-inhomogeneous problem by converting it into an infinite-horizon and time-homogeneous problem that incorporates time t as a special state variable (Appendix I). The state-action space S × A can be any domain on which we can define a reproducing kernel Hilbert space (RKHS), either discrete or continuous.
Assume the model P is unknown, but we can observe a set of transition pairsD n = (s i , a i , r i , s i ) n i=1 that obeys P in that (r i , s i ) ∼ P(· | s i , a i ). The goal is to construct an interval [Ĵ − (D n ),Ĵ + (D n )] based on the observed data such that the probability that it contains the true J * is no smaller than a user-specified confidence level. We consider the challenging case where the data is off-policy, behavior-agnostic, and time-dependent, as formally specified in the following.
Assumption 2.1 (Data Assumption). Assume the dataD n = (s i , a i , r i , s i ) n i=1 is drawn from an unknown joint distribution D 1:n on (S × A × R × S) n that satisfies D 1:n (r i , s i | s i , a i ;D <i ) = P(r i , s i | s i , a i ), ∀i = 1, . . . , n, whereD <i := (s j , a j , r j , s j ) j<i .
Goal Given a datasetD n , a target policy π, and a confidence level δ ∈ (0, 1), we want to construct an interval [Ĵ − (D n ),Ĵ + (D n )] ⊂ R, such that where Pr(·) is w.r.t. the randomness of the dataD n drawn from the distribution D 1:n .
The data assumption provides a partial specification of the joint distribution D 1:n , which only requires that (r i , s i ) should be generated from P(· | s i , a i ) given (s i , a i ) ∪D <i at each step, while imposing no constraints on how (s i , a i ) is generated based onD <i . This provides great flexibility in terms of the data collection procedure. For example,D n can be either collected from a single MDP governed by an unknown, time-varying, non-Markovian behavior policy (in which case (s i , a i ) = (s i+1 , a i+1 )), or be a combination of many short MDP segments (in which case (s i , a i ) may not equal (s i+1 , a i+1 )). Note that this significantly relaxes the data assumptions in recent works (Liu et al., 2018a;Mousavi et al., 2020;, which require (s i , a i ) n i=1 to be independent or i.i.d. Eq (3) is a correctness requirement of the confidence interval. In addition, we also want to make the interval estimation as tight as possible. Specifically, it is desirable that the length (Ĵ + (D n ) −Ĵ − (D n )) of the interval vanishes to zero as n → ∞, ideally with a fast converge rate (e.g., O(n −1/2 )). While assumption 2.1 allows us to construct bounds that are correct and non-vacuous, it is too mild to give meaningful a priori bound on the tightness of the interval estimation. For example, under Assumption 2.1, it is possible that the whole dataset is the replication of a single transition pair, i.e., (s i , a i , r i , s i ) = (s 1 , a 1 , r 1 , s 1 ), in which case we cannot have a 1 − δ confidence interval whose length vanishes to zero as n → ∞.
A primary advantage of having an interval estimation is that we know a posteriori length of the confidence interval once we construct it, which provides an instance-based uncertainty estimation without making any additional assumption beyond Assumption 2.1. If a stronger assumption is imposed (e.g., (s i , a i , r i , s i ) n i=1 are i.i.d. or weakly dependent in a proper sense for different i), we can derive a priori estimation of the state action reward state-action next state-action next state-action + reward Table 1: Summary of the notation. We assume we observe pairs of (s i , a i , r i , s i ), and for each i, we augment the data with the next action a i as part of the algorithm by drawing from π(· | s i ), where π is the known target policy. The data assumption in Eq (2) is equivalent to saying that y i ∼ P π (· | x i ); see Eq (4).
convergence rate of the tightness (Ĵ + (D n ) −Ĵ − (D n )) for the confidence interval as shown in Section 6.3 and Appendix D.
A crucial fact is that the data assumption above implies a martingale structure on the empirical Bellman residual operator of the Q-function. As we will show in Section 5, this enables us to derive a key concentration inequality underpinning our non-asymptotic confidence bound.

Notation
We define a few notations that will simplify the presentation in the rest of work. First of all, to facilitate the definition of the empirical Bellman operator in the sequel (see Section 4.1 and Remark 4.1), we append each (s i , a i , r i , s i ) with a random action a i drawn from the target policy π(· | s i ) following s i . This can be done for free as long as the target policy π is given. Also, we write Correspondingly, define X = S × A to be the state-action space. Then the observed data can be written as pairs of {x i , y i } n i=1 . We equalize the dataD n with its empirical measurê Then Assumption 2.1 is equivalent to saying that y i ∼ P π (· | x i ). In this way, our setting can be viewed as a supervised learning problem on {x i , y i } with an unknown P π . The difference is that we are interested in estimating J π,P , which is a nonlinear functional of P π , rather than the whole model P π (which would yield a model-based method; see discussions in Section 3).
We use the star notation to label quantities that depend on the unknown parameters (such as J * ), and the hat notation to denote quantities that depend on observed data (such asD n andĴ ± (D n )). The prime notation always denotes the next state (such as s i ).

Related Work
We give an overview of different approaches for uncertainty estimation in OPE.
Finite-Horizon Importance Sampling (IS) Assume the data is collected by rolling out a known behavior policy π 0 up to a trajectory length T , then we can estimate the finite horizon reward by changing E π,P [·] to E π 0 ,P [·] with importance sampling(e.g., Precup, 2001;Thomas et al., 2015a,b). Taking the trajectory-wise importance sampling as an example, assume we collect a set of independent trajectories τ i := {s i t , a i t , r i t } T −1 t=0 , i = 1, . . . , m up to a trajectory length T by unrolling a known behavior policy π 0 . When T is large, we can estimate J * by a weighted averaging: One can construct non-asymptotic confidence bounds based onĴ IS using variants of concentration inequalities (Thomas, 2015;Thomas et al., 2015b). Unfortunately, a key problem with this IS estimator is that the importance weight ω(τ i ) is a product of the density ratios over time, and hence tends to cause an explosion in variance when the trajectory length T is large. Although improvement can be made by using per-step and self-normalized weights (Precup, 2001), or control variates (Jiang and Li, 2016;Thomas and Brunskill, 2016), the curse of horizon remains to be a key issue to the classical IS-based estimators (Liu et al., 2018a).
Moreover, due to the time dependency between the transition pairs inside each trajectory, the non-asymptotic concentration bounds can only be applied on the trajectory level and hence decay with the number m of independent trajectories in an O(1/ √ m) rate, though m can be small in practice. We could in principle apply the concentration inequalities of Markov chains (e.g., Paulin, 2015) to the time-dependent transition pairs, but such inequalities require to have an upper bound of certain mixing coefficient of the Markov chain, which is unknown and hard to construct empirically. Our work addresses these limitations by constructing a non-asymptotic bound that decay with the number n = mT of transitions pairs, while without requiring known behavior policies and independent trajectories.  (2020), as well as the DICE-family (e.g., Nachum et al., 2019a,b;Zhang et al., 2020a;Wen et al., 2020;. These methods are based on either estimating the value function, or the stationary visitation distribution, which is shown to form a primal-dual relation (Tang et al., 2020a;Uehara et al., 2020;) that we elaborate in depth in Section 4.
Besides Feng et al. (2020) which directly motivated this work, there has been a recent surge of interest in interval estimation under infinite-horizon OPE (e.g., Liu et al., 2018b;Duan et al., 2020;Feng et al., 2020;Tang et al., 2020b;Lazic et al., 2020). For example,  develop an asymptotic confidence bound (CoinDice) for DICE estimators with an i.i.d assumption on the off-policy data; Duan et al. (2020) provide a data dependent confidence bounds based on Fitted Q iteration (FQI) using linear function approximation when the off-policy data consists of a set of independent trajectories; Jiang and Huang (2020) provide a minimax method closely related to our method but do not provide analysis for data error; Tang et al. (2020b) propose a fixed point algorithm for constructing deterministic intervals of the true value function when the reward and transition models are deterministic and the true value function has a bounded Lipschitz norm.
Model-Based Methods Since the model P is the only unknown variable, we can construct an estimatorP of P using maximum likelihood estimation or other methods, and plug it into (1) to obtain a plug-in estimator J = J π,P . This yields the model-based approach to OPE (e.g., Jiang and Li, 2016;Liu et al., 2018b). One can also estimate the uncertainty in J π,P by propagating the uncertatinty inP (e.g., Asadi et al., 2018;Duan et al., 2020), but it is hard to obtain non-asymptotic and computationally efficient bounds unlessP is assumed to be simple linear models. In general, estimating the whole model P can be an unnecessarily complicated problem as an intermediate step of the possibly simpler problem of estimating J π,P .
Bootstrapping, Bayes, Distributional RL As a general approach of uncertainty estimation, bootstrapping has been used in interval estimation in RL in various ways (e.g., White and White, 2010;Hanna et al., 2017;Kostrikov and Nachum, 2020;Hao et al., 2021). Bootstrapping is simple and highly flexible, and can be applied to time-dependent data (as appeared in RL) using variants of block bootstrapping methods (e.g., Lahiri, 2013; White and White, 2010). However, bootstrapping typically only provides asymptotic guarantees; although non-asymptotic bounds of bootstrap exist (e.g., Arlot et al., 2010), they are sophistic and difficult to use in practice and would require to know the mixing condition for the dependent data. Moreover, bootstrapping is time consuming since it requires to repeat the whole off-policy evaluation pipeline on a large number of resampled data.
Bayesian methods (e.g., Engel et al., 2005;Ghavamzadeh et al., 2016;Yang et al., 2020a) offer another general approach to uncertainty estimation in RL, but require to use approximate inference algorithms and do not come with non-asymptotic frequentist guarantees. In addition, distributional RL (e.g., Bellemare et al., 2017) seeks to quantify the intrinsic uncertainties inside the Markov decision process, which is orthogonal to the epistemic uncertainty that we consider in off-policy evaluation.

Two Dual Approaches to Infinite-Horizon Off-Policy Estimation
The main idea of infinite-horizon OPE is to transform the estimation of the expected reward into estimating either the Q-function or the visitation distribution (or its related density ratio) by exploiting the stationary property of the MDP under π. This section discusses these two tightly connected methods, which form a primal-dual relation and together lay out a foundation for our main confidence bounds.
The Q-function associated with policy π and model P is defined as where the expectation is taken when we execute π under model P initialized from a fixed state-action pair x 0 = (s 0 , a 0 ). Let D π,t be the distribution of (x t , y t ) = (s t , a t , s t , a t , r t ) when executing policy π starting from s 0 ∼ D 0 for t steps. The visitation distribution of π is defined as Note that D * integrates to 1/(1 − γ), although we still treat it as a probability measure in the notation.
The expected reward J * can be expressed using either q * or D * as follows: where r ∼ D * (resp. x ∼ D π,0 ) denotes sampling from the r-(resp. x-) marginal distribution of D * (resp. D π,0 ). Eq. (6) transforms the estimation of J * into estimating either q * or D * and plays a key role in infinite-horizon OPE.

Value Estimation via Q Function
Because D π,0 (x) = D 0 (s)π(a|s) is known, we can estimate J * by E x∼D π,0 [q(x)] with any estimationq of the true Q-function q * ; the expectation under x ∼ D π,0 can be estimated to any accuracy with Monte Carlo.
To estimate q * , we consider the empirical and expected Bellman residual operator: It is well-known that q * is the unique solution of the Bellman equation Rq = 0. Since y i ∼ P π (·|x i ) for each data point inD n , if q = q * , thenRq(x i , y i ) is a zero-mean random variable conditional on x i . Let ω be any function from X to R, then iR q(x i , y i )ω(x i ) also has zero mean. This motivates the following functional Bellman loss (Feng et al., 2019(Feng et al., , 2020Xie and Jiang, 2020), where W is a set of functions ω : X → R. To ensure that the sup is positive and finite, W is typically set to be a unit ball of some normed function space W, that is, Feng et al. (2019) consider the simple case when W is the unit ball K of the reproducing kernel Hilbert space (RKHS) with a positive definite kernel k : X × X → R, for which the loss has a simple closed form: Note that the RHS of Eq. (9) is the square root of the kernel Bellman V-statistics in Feng et al. (2019). Feng et al. (2019) showed that, when the support of the data distributionD n covers the state-action space X (which requires an infinite data size when the domain size is infinite) and k is an integrally strictly positive definite kernel, we have L K (q;D n ) = 0 iff q = q * . Therefore, one can estimate q * by minimizing L K (q,D n ).
Remark 4.1. The empirical Bellman residual operatorR can be extended tô where {a } m i=1 are i.i.d. drawn from π(·|s ). As m increases, this gives an unbiased estimator of Rq with lower variance. If m = +∞, we havê which coincides with the operator used in the expected SARSA (Sutton and Barto, 1998). Therefore, the choice of m provides a trade-off between accuracy and computational cost. Without any modification, all results in this work can be applied to theRq in (10)-(11) for any positive integer m. We use m = 5 in most of our experiments.

Value Estimation via Visitation Distribution
Another way to estimate J * in Eq. (6) is to approximate D * with a weighted empirical measure of the data (Liu et al., 2018a;Nachum et al., 2019a;Mousavi et al., 2020;Zhang et al., 2020a). The key idea is to assign an importance weight ω(x i ) to each data point x i inD n . We can choose the function ω : X → R properly such that D * and hence J * can be approximated by the ω-weighted empirical measure ofD n (and reward) as follows: Intuitively, ω can be viewed as the density ratio between D * andD n , although the empirical measureD n may not have well-defined density. Liu et al. (2018a); Mousavi et al. (2020) proposed to estimate ω by minimizing a discrepancy measure betweenD ω n and D * . To see this, note thatD ω n = D * if ∆(D ω n , q) = 0 for any function q, with ∆ defined as where we use the fact that Liu et al., 2018a). Also note that the RHS of Eq. (13) can be practically calculated given any ω and q without knowing D * . Let Q be a set of functions q : X → R. One can define the following loss for ω: Similar to L K (q;D n ), when Q equals the unit ballK of the RKHS of a positive definite kernelk(x,x), Eq. (14) can be expressed into a quadratic closed form shown in Mousavi et al. (2020): , is obtained by applyingTx andT x sequentially by treating k as a function ofx and then of x.

Primal and Dual Deterministic Bounds (Infinite Data Case)
The two types of estimators above are connected in a primal-dual fashion, which can be used to build deterministic upper and lower bounds of J * that can be viewed as the limit of our empirical bounds with an infinite data size n (Tang et al., 2020a;. We provide an overview of this framework and derive a number of bounds whose empirical version will be discussed in depth in Section 6.
Let Q be a function set that includes the true Q-function q * . It is known that the J * can be represented using the following optimization problem: which holds because J * = E D π,0 [q * ] and q = q * is the unique solution of Rq(x) = 0, ∀x ∈ X .
Using Lagrange duality, Eq (16) is equivalent to where ω(x)D ∞ (x) serves as the Lagrange multiplier: D ∞ is a fixed distribution whose support is X and ω : X → R is optimized in the set of all functions such that the objective is defined.
Although D ∞ was technically introduced as a part of the Lagrange multiplier, it should be intuitively viewed as a population data distribution, the limit of the empirical dataD n as n → ∞; in practice, we replace D ∞ withD n when using the bounds in this section. However, note that Assumption 2.1 does not prescribe the existence of such D ∞ from the dataD n (since it does not consider the limit when n → ∞). We need additional assumptions, such as when ( or follows an ergodic Markov chain to relateD n with a well defined limit D n . Such additional assumption is needed when ensuring a-priori bounds on the length of the confidence interval that we construct; see Theorem D.1. We can obtain an upper bound of J * from the minimax representation in Eq. (18) by constraining the optimization domain of ω to a normed function space W whose unite ball is W, where the first inequality is due to the constraint on W; the second inequality is due to exchanging the order of sup and inf, which turns to equality if strong duality holds.
Importantly, the Lagrange function M (q, ω; D ∞ ) is connected to both L W (q, D ∞ ) and I Q (ω; D ∞ ), which allows us to derive a pair of primal and dual bounds whose non-asymptotic version will be derived in this work.
By recognizing λ as a scalar Lagrange multiplier, we have which can be seen as a relaxation of Eq. (16) since L W (q, D ∞ ) ≤ 0 (equivalent to L W (q, D ∞ ) = 0) can be a weaker constraint than {Rq(x) = 0, ∀x ∈ X } if W is a small function set.
On the other hand, plugging Eq. (21) into Eq. (19), we have Here Eq. (22) and Eq. (23) are dual to each other and provide bounds of J * in terms of q and ω, respectively. In Section 6, we provide empirical variants of Eq. (22) and Eq. (23) which replace D ∞ withD n while providing non-asymptotic bounds.
Doubly Robust Estimation and the Lagrangian Another related key feature of the Lagrange function M (q, ω; D ∞ ) above is the following "double robustness" property: where q * is the true Q-function, and ω * is the density ratio between D * and D ∞ such that D ω ∞ = D * . The double robustness in Eq. (24) says that M (q, ω; D ∞ ) equals J * if either q = q * or ω = ω * holds. Therefore, letq andω be estimations of q * and ω * respectively, thenĴ dr := M (q,ω;D n ) yields doubly robust estimation of J * in the sense thatĴ dr forms an accurate estimation of J * if eitherq orω is accurate; see more discussion in Tang et al. (2020a).
All the bounds in this Section depend on D ∞ which need to be replaced by the empirical dataD n in practice. One difficulty, however, is that Assumption 2.1 does not directly imply thatD n converges to a fixed limit distribution D n as n → ∞. Therefore, special care is taken in Section 5 and Section 6 to sidestep the introduction of D ∞ as we construct non-asymptotic confidence bounds dependent onD n .

Concentration Inequality of Kernel Bellman Loss
To quantify the data error, we establish in this section a concentration inequality to bound the deviation of the kernel Bellman loss (KBL) L K (q * ;D n ) away from zero, which serves as a key building block of the non-asymptotic bounds in Section 6. This inequality holds under the mild data assumption 2.1 and does not require additional i.i.d. or independence assumption thanks to a martingale structure implied Assumption 2.1.
We first introduce the following semi-expected kernel Bellman loss (KBL) which L K (q;D n ) concentrates around for any q: where we replace the empirical Bellman residual operatorRq in Eq. (9) with its expected counterpart Rq, but still keep the empirical average over {x i } n i=1 inD n . For a more general function set W, we can similarly define L * W (q;D n ) by replacingRq with Rq in Eq. (8). Note that we have L * W (q;D n ) = 0 when q = q * for anyD n and any W.
Theorem 5.1 below shows that L K (q;D n ) concentrates around L * K (q;D n ) with an O(n −1/2 ) error under Assumption 2.1. At a first glance, it may seem surprising that the concentration bound can hold even without any independence assumption between {x i }. An easy way to make sense of this is by recognizing that the randomness in y i conditional on x i is aggregated through averaging, even if {x i } are deterministic.
Theorem 5.1. Assume K is the unit ball of RKHS with a positive definite kernel k(·, ·). Let c q,k := sup x,y (Rq(x, y) − Rq(x)) 2 k(x, x) < ∞. Under Assumption 2.1, for any δ ∈ (0, 1), with at least probability 1 − δ, we have In particular, when q = q * , we have c q * ,k = sup x,y (Rq * (x, y)) 2 k(x, x), and An upper bound of the coefficient c q,k can be calculated easily in practice; see Feng et al. (2020) and Appendix C.1.
Intuitively, to see why we can expect an O(n −1/2 ) bound, note that L K (q,D n ) consists of the square root of the product of twoRq terms, each of which contributes an O(n −1/2 ) error w.r.t. Rq. Technically, the proof is made possible by observing that Assumption 2.1 ensures that The proof also leverages a special property of RKHS and applies a Hoeffding-like inequality by Pinelis (1992) on Hilbert spaces. See Appendix C for details. For other more general function sets W, we establish in Appendix F a similar bound by using Rademacher complexity, although it requires to know an upper bound of the Rademacher complexity of a function set associated with W and yields a less tight bound than Eq. (26) if W = K.
IfD n weakly converges to a limit D ∞ as n → ∞, we can expect that L K (q,D n ) converges to L K (q, D ∞ ). However, D ∞ is not implied from Assumption 2.1. Theorem 5.1 sidesteps the introduction of D ∞ thanks to the semi-expected KBL L * K (q,D n ). Also, as Assumption 2.1 does not impose any (weak) independence between {x i }, without introducing further assumptions, we cannot establish that L K (q;D n ) concentrates around the full expectation ED n [L K (q;D n ) 2 ] 1/2 whereD n is averaged w.r.t. the underlying data generation distribution D 1:n .

Primal-Dual Non-Asymptotic Confidence Bounds
We are ready to extend the bounds in Eq. (22)-Eq. (23) to finite data case. To avoid introducing D ∞ , we start with building an empirical counterpart of bound (16) in Section 6.1, and then proceed to derive the non-asymptotic counterparts of bound (22) in Section 6.2 and bound (23) in Section 6.3.

A Data-Dependent Oracle Bound
Let Q be a function set that contain the true Q-function q * , that is, q * ∈ Q.
Given a datasetD n , we have the following upper bound of J * that generalizes (16): By definition, this is the tightest upper bound given only the information inD n through the empirical Bellman operator, as in this case q and q * would look indistinguishable ifRq(x i , y i ) =Rq * (x i , y i ) for all i. The lower boundĴ − Q, * can be defined analogously by replacing sup q∈Q with inf q∈Q . Note thatĴ + Q, * andĴ − Q, * are still not practically computable because they depends on the unknown q * .
Note that this result holds trivially because q * is included in the optimization domain in Eq. (28). It does not require any assumption on the dataD n = {x i , y i }. That is, J * ∈ [Ĵ − Q, * ,Ĵ + Q, * ] remains to be true (but may not be useful) even ifD n is filled with random numbers irrelevant to P and π; the bound becomes useful if we have sufficiently number of data points satisfying y i ∼ P π (· | x i ).
Introducing Q is Necessary It is necessary to introduce the function Q to ensure a finite bound unless S is a small discrete set. Removing the constraint of q ∈ Q in Eq. (28) would lead to an infinite upper/lower bound, unless the {s i , s i } n i=1 pairs from the dataD n almost surely covers the whole state space S.
Proposition 6.2. Define Q π,∞ = {q :Rq(x i , y i ) =Rq * (x i , y i ), ∀i = 1, . . . , n}. Then, unless Pr s∼D π,0 (s Note that Pr s∼D π,0 (s / ∈ {s i , s i } n i=1 ) = 0 can hold only when the data size n is no smaller than the cardinality of the state space S, which is infinite when S is a continuous domain. Therefore, it is necessary to introduce a Q that is smaller than Q π,∞ unless S is a discrete set with a small number of elements.
Note that every element of the Q π,∞ defined in Proposition 6.2 is indistinguishable with q * from the information accessible through the empirical Bellman operatorR. Thus, unless Q = Q π,∞ , which would make it too large to be useful, we can not provably guarantee that q * ∈ Q, because every element in Q π,∞ \ Q can have a chance to be the true Q-function q * , where \ denotes the set minus operator. Therefore, the correctness of our bound will ultimately relies on an un-checkable model assumption, which is unavoidable in many statistical estimation problems in general.
On the other hand, it is possible to empirically reject a poorly chosen Q by hypothesis testing when Q ∩ Q π,∞ = ∅. Consider the following test: Null: q * ∈ Q vs. Alternative : q * ∈ Q.
We can reject q * ∈ Q with a false positive error δ if inf q∈Q L W (q;D n ) ≥ ε n , where ε n satisfies Pr(L W (q * ;D n ) ≤ ε n ) ≥ 1 − δ. This is because When Q is a finite ball in RKHS, we can practically solve min q∈Q L W (q;D n ) by using the finite representer theorem of RKHS (Scholkopf and Smola, 2018); see related discussion in Section 6.2 and Appendix E.

Empirical Primal Bound via Functional Bellman Loss
To make use of the concentration inequality in Section 5, we relax Eq. (28) to the following bound that depends on a function set W: andĴ − Q,W, * is defined analogously. Obviously,Ĵ + Q,W, * is an upper bound ofĴ + Q, * because the optimization domain in Eq. (29) is larger than that of Eq. (28).    Table 2.Ĵ + Q,W andĴ ++ Q,W can be calculated from empirical data, whileĴ + Q, * andĴ + Q,W, * cannot since they depend on the unknown q * .
Note thatĴ + Q,W, * is still not practically computable because it depends on the unknown L W (q * ;D n ). However, because L W (q * ;D n ) concentrates around zero whenD n satisfies Assumption 2.1, we can construct a fully empirical bound as follows: where ε n is a positive number properly chosen such that We set ε n = 2c q * ,k log(2/δ)/n when W = K following Theorem 5.1. The lower boundĴ − Q,W can be defined analogously. Obviously, Eq. (31) is the empirical counterpart of Eq. (22). Proposition 6.4. Assume Eq. (32) holds. Then for any set Q, we have Further, if q * ∈ Q, we have Computation ofĴ + Q,W in Eq. (31) If Q is taken to be an RKHS ball, the optimization in Eq. (31) can be shown to reduce to a finite dimensional convex optimization, and hence can be solved in practice. Precisely, if Q is a finite ball of the RKHS associated with a positive definite kernelk(·, ·) (which should be distinguished with kernel k(·, ·) of K), then by the finite representer theorem of RKHS (Scholkopf and Smola, 2018), the global optimum of (31) can be achieved by a function of form
Plugging this into Eq. (31), the optimization can be shown to reduce to a convex optimization on with a linear objective and quadratic inequality constraint.
Unfortunately, when the data size n is large, solving Eq. (31) still leads to a high computational cost. Importantly, because the guarantee in Eq. (33) only holds when the maximization is solved to global optimality, fast approximation methods, such as random feature approximation (Rahimi and Recht, 2007), should not be used in principle. We address this problem by considering the dual form of Eq. (31), which avoids to solve the challenging global optimization in Eq. (31). Moreover, the dual form enables us to better understand the tightness of the confidence interval and issues regarding the choices of Q and W.

The Dual Bound
To derive the dual bound, let us plug the definition of L W (q;D n ) into Eq. (8) and introduce a Lagrange multiplier λ:Ĵ where we assume that W is the unit ball of the normed space W and hence can write any ω in W into ω(x) = λh(x) such that h ∈ W and λ = ||w|| W . Exchanging the order of min/max and some further derivation yields the following main result.
Theorem 6.5. Let W be the unit ball of a normed function space W. We have Here −Q = {−q : q ∈ Q} and hence I −Q (ω;D n ) = I Q (ω;D n ) if Q = −Q.
Further, the bound is tight, that is,Ĵ + Q,W = inf ω∈ WF + Q (ω) andĴ − Q,W = sup ω∈ WF − Q (ω), if Q is convex and there exists a function q ∈ Q that satisfies the strict feasibility condition that L W (q;D n ) < ε n . Therefore, when Eq. (32) holds, for any function set Q, and any function ω + , ω − ∈ W (the choice of Q, ω + , ω − can depend onD n arbitrarily), we have Theorem 6.5 transforms the original bound in Eq. (31), framed in terms of q and L W (q;D n ), into a form that involves the density-ratio ω and the related loss I Q (ω;D n ). The bounds in Eq. (36) can be interpreted as assigning an error bar around the ω-based estimatorĴ ω = ED ω n [r] in Eq. (12), with the error bar of I ±Q (ω;D n ) + ε n ω W . Specifically, the first term I ±Q (ω;D n ) measures the discrepancy betweenD ω n and D * as discussed in Eq. (14), whereas the second term captures the randomness in the empirical Bellman residual operatorRq * .
Compared with Eq. (31), the global maximization on q ∈ Q is now transformed inside the I Q (ω;D n ) term, which yields a simple closed form solution when Q is a finite ball in RKHS. We can optimize ω + and ω − by minimizing/maximizingF + Q (ω) andF − Q (ω) to obtain the tightest possible bound (and hence recover the primal bound). However, it is not necessary to find the exact globally optimal solutions for practical purpose. When W is an RKHS, by the standard finite representer theorem (Scholkopf and Smola, 2018), the optimization on ω reduces to a finite dimensional optimization, which can be approximately solved with any practical technique without sacrificing the correctness of the bound (although the optimization quality of ω impacts the tightness of the bound). We elaborate on this in Appendix E.
Length of the Confidence Interval The form in Eq. (36) also makes it much easier to analyze the tightness of the confidence interval. Suppose ω = ω + = ω − and Q = −Q, the length of the optimal confidence interval isĴ Given that ε n is O(n −1/2 ), we can make the overall length of the optimal confidence interval also O(n −1/2 ) if W is rich enough to include a good density ratio estimator ω * that satisfies I Q (ω * ;D n ) = O(n −1/2 ) and has a bounded norm ω * W . Assumption 2.1 does not ensure the existence of such ω * . However, we can expect to have such a ω * if 1) Q has an O(n −1/2 ) sequential Rademacher complexity (Rakhlin et al., 2015) (which holds if Q is a finite ball in RKHS); and 2)D n is collected following a Markov chain with a strong mixing condition and weakly converges to some limit distribution D ∞ whose support is X ; in this case, we can define ω * as the density ratio between D * and D ∞ . See Appendix D for more discussions. Indeed, our experiments show that the lengths of practically constructed confidence intervals do tend to decay with an O(n −1/2 ) rate approximately.
Using Data-Dependent W and Q To ensure Eq. (33) holds, the choice of W cannot depend on the datâ D n , since it may introduce additional dependency and hence make the concentration inequality in Theorem 5.1 invalid. Therefore, if we want to use a data-dependent W, we need to either base the construction of W on separate holdout data, or introduce a generalization bound to account the dependence of W on the dataD n .
Interestingly, Eq. (33) holds even if we take Q = Q(D n ) to be an arbitrary function of the dataD n . This is because Q is irrelevant to the concentration inequality (Eq. (32)). This justifies that one can construct Q adaptively based on the data to get tighter confidence interval. For example, we can make Q an RKHS ball centering around an estimatorq ≈ q * given by a state-of-the-art method (e.g., fitted iteration or modelbased methods). A caveat is that we will need to ensure that Pr(q * ∈ Q(D n )) = 1 in order to ensure Pr J * ∈ Ĵ − Q(Dn),W ,Ĵ + Q(Dn),W ≥ 1 − δ (see Proposition 6.4), which, as we discussed in Section 6.1, can not be provably guaranteed based on only empirical observation.

Experiments
We compare our method with a variety of existing algorithms for obtaining asymptotic and non-asymptotic bounds on a number of benchmarks. We find our method can provide confidence interval that correctly covers the true expected reward with probability larger than the specified success probability 1 − δ (and is hence safe) across the multiple examples we tested. In comparison, the non-asymptotic bounds based on IS provide much wider confidence intervals. On the other hand, the asymptotic methods, such as bootstrap, despite giving tighter intervals, often fail to capture the true values with the given probability in practice.

Environments and Dataset Construction
We test our method on three environments: Inverted-Pendulum and CartPole from OpenAI Gym (Brockman et al., 2016), and a Type-1 Diabetes medical treatment simulator. 1 We follow a similar procedure as Feng et al. (2020) to construct the behavior and target policies. more details on environments and data collection procedure are included in Appendix H.1.

Algorithm Settings
We test the dual bound described in Algorithm 1. Throughout the experiment, we always set W = K, the unit ball of the RKHS with positive definite kernel k, and set Q = r QK , the ball of radius r Q in the RKHS with another kernelk. We take both kernels to be Gaussian RBF kernel and choose r Q and the bandwidths of k andk using the procedure in Appendix H.2. We use a fast approximation method to optimize ω in F + Q (ω) and F − Q (ω) as shown in Appendix E. Once ω is found, we evaluate the bound in Eq. (36) exactly to ensure that the theoretical guarantee holds.

Baseline Algorithms
We compare our method with four existing baselines, including the IS-based nonasymptotic bound using empirical Bernstein inequality by Thomas et al. (2015b), the IS-based bootstrap bound of Thomas (2015), the bootstrap bound based on fitted Q evaluation (FQE) by Kostrikov and Nachum (2020), and the bound in Feng et al. (2020) which is equivalent to the primal bound in (31) but with looser concentration inequality (they use a ε n = O(n −1/4 ) threshold). Results Figure 1 shows our method obtains much tighter bounds than Feng et al. (2020), which is because we use a much tighter concentration inequality, even the dual bound that we use can be slightly looser than the primal bound used in Feng et al. (2020). Our method is also more computationally efficient than that of Feng et al. (2020) because the dual bound can be tightened approximately while the primal bound requires to solve a global optimization problem. Figure 1 (b) shows that we provide increasingly tight bounds as the data size n increases, and the length of the interval decays with an O(n −1/2 ) rate approximately. Figure 1 (c) shows that when we increase the significance level δ, our bounds become tighter while still capturing the ground truth. Figure 1 (d) shows the percentage of times that the interval fails to capture the true value in a total of 100 random trials (denoted asδ) as we vary δ. We can see thatδ remains close to zero even when δ is large, suggesting that our bound is very conservative. Part of the reason is that the bound is constructed by considering the worse case and we used a conservative choice of the radius r Q and coefficient c qπ,k in Eq. (27) (See Appendix H.2).
In Figure 2 we compare different algorithms on more examples with δ = 0.1. We can again see that our method provides tight and conservative interval that always captures the true value. Although FQE (Bootstrap) yields tighter intervals than our method, it fail to capture the ground truth much more often than the promised δ = 0.1 (e.g., it fails in all the random trials in Figure 2 (a)).
We conduct more ablation studies on different hyper-parameter and data collecting procedure. See Appendix H.2 and H.3 for more details.

Conclusion and Future Directions
We develop a practical algorithm for constructing non-asymptotic confidence intervals for infinite-horizon off-policy evaluation with a very mild data assumption that holds for behavior-agnostic and time-dependent data. Our work opens a number of future direction: how to apply our bounds to develop new methods for safe policy optimization and safe exploration? how to develop rigorous and practical approaches to select Q and W adaptively based on the observed data, including their kernel, bandwidth and the radius of ball? how to extend our method to obtain bounds for different target policies π (see e.g., ) and initial distributions D 0 simultaneously and computationally efficiently, without repeatedly solving the optimization in each case? The oracle bound in Eq. (28) provides a notion of oracle bound given the information drawn from the empirical Bellman operator. Is it possible to obtain even tighter bounds by exploiting other information from data, by e.g., model-based methods?
A Proof of Lemma 4.2 Proof. By assumption, we can write any ω 0 ∈ W into ω 0 (x) = λω(x) where λ ≥ 0 and ω ∈ W. Therefore, following Eq (18), we have where we used the definition that L W (q, For Eq (21), we note that where we use the definition of ∆(D ω ∞ , q) in Eq (13). Therefore, B Proof of the Dual Bound in Theorem 6.5 Proof. Note that we assumed that W is the unit ball of the normed space W. Therefore, we can write any ω in W into ω(x) = λh(x) with λ ∈ R, h ∈ W and ω W = λ.
Using Lagrange multiplier, the bound in Eq. (31) is equivalent tô Then we have The lower bound follows analogously. The strong duality holds when the Slater's condition (Nesterov, 2003) is satisfied, which amounts to saying that the primal problem in Eq. (31) is convex and strictly feasible; this requires that Q is convex and there exists at least one solution q ∈ Q that satisfies that constraint strictly, that is, L W (q;D n ) < ε n ; note that in our case the objective function Q is linear on q and the constraint function L W (q;D n ) is always convex on q (since it is the sup a set of linear functions on q following Eq. (8)).

C Proof of Concentration Bound in Theorem 5.1
Our proof requires the following Hoeffding inequality on Hilbert spaces by Pinelis (Theorem 3, 1992); see also Section 2.4 of Rosasco et al. (2010).
Lemma C.1. (Theorem 3, Pinelis, 1992) Let H be a Hilbert space and {f i } n i=1 is a martingale difference sequence in H that satisfies sup i f i H ≤ σ almost surely. We have for any > 0, Therefore, for δ ∈ (0, 1), with probability at least 1 − δ, we have 1 Lemma C.2. Let k(x, x ) be a positive definite kernel whose RKHS is H k . Define Assume Assumption 2.1 holds, then {f i } n i=1 is a martingale difference sequence in H k w.r.t. T <i := (x j , y j ) j<i ∪ (x i ). That is, E [f i+1 (·) | T <i ] = 0 for i = 1, . . . , n. In addition, and f i 2 H k ≤ c q,k for ∀i = 1, . . . , n.
Proof of Theorem 5.1. Following Lemma C.1 and Lemma C.2, since {f i } n i=1 is a martingale difference sequence in H k with f i 2 H k ≤ c q,k almost surely, we have with probability at least 1 − δ, Using Lemma C.3 below, we have This completes the proof.
Lemma C.3. Assume k(x, x ) is a positive definite kernel. We have Proof. Defineĝ Rq(x i )k(x i , ·).
Then we have The result then follows the triangle inequality ĝ H k − g H k ≤ ĝ − g H k .
C.1 Calculation of c q * ,k The practical calculation of the coefficient c q * ,k in the concentration inequality was discussed in Feng et al. (2020), which we include here for completeness.
Lemma C.4. (Feng et al. (2020) Lemma 3.1) Assume the reward function and kernel function is bounded with sup x |r(x)| ≤ r max and sup x,x |k(x, x )| ≤ K max , we have: In practice, we evaluate K max from the kernel function that we choose (e.g., K max = 1 for RBF kernels), and r max from the knowledge of the environment.

D More Discussion on the Tightness of the Confidence Interval
The benefit of having both upper and lower bounds is that we can empirically access the tightness of the bound by checking the length of the interval [F − Q (ω − ),F + Q (ω + )]. However, from the theoretical perspective, it is desirable to know a priori that the length of the interval will decrease with a fast rate as the data size n increases. We now show that this is the case if W is chosen to be sufficiently rich so that it includes a ω ∈ W such thatD ω n approximates D * closely in a proper sense.
Theorem D.1. Assume W is sufficiently rich to include a "good" ω * in W withD ω * n ≈ D * in the sense that where c and α are two positive coefficients. Then we have Assumption (39) holds ifD n is collected following a Markov chain that has certain strong mixing condition and weakly converges to a limit continuous discussion D ∞ whose support is X , and the density ratio between D * and D ∞ , denoted by ω * , is included in W. In this case, if Q is a finite ball in RKHS, then we can achieve Eq. (39) with α = 1/2, and yields the overall bound of rate O(n −1/2 ). For more general function classes, α depends on the martingale Rademacher complexity of the function setRQ = {Rq(x, y) : q ∈ Q} (Rakhlin et al., 2015). In our empirical reults, we observe that the lengths of the practically constructed confidence intervals do tend to follow the O(n −1/2 ) rate approximately; see Figure 1(b) Proof. Note that γq(x ) − q(x) = −Rq(x, y) − r and hene Because ω * ∈ W, we havê The case of lower bound follows similarly.

E Practical Optimization on W
Consider the optimization of ω in W, Assume W is the RKHS of kernel k(x,x). By the finite representer theorem of RKHS (Scholkopf and Smola, 2018), the optimization of ω ∈ W can be reduced to a finite dimensional optimization problem. Specifically, the optimal solution of (40) can be written into a form of ω( The optimization of ω reduces to solving the following optimization on α: When Q is a finite ball of an RKHS (that is different from W), we can calculate I Q (Kα;D n ) using Eq. (15). This computation can be still expensive when n is large. Fortunately, our confidence bound holds correctly for any ω; better ω only gives tighter bounds, but it is not necessary to find the exact global optimal ω. Therefore, one can use any approximation algorithm to find ω, which provides a trade-off of tightness and computational cost. We discuss two methods: 1) Approximating α The length of α can be too large when n is large. To address this, we assume α i = g(x i , θ), where g is any parametric function (such as a neural network) with a parameter θ; assume θ has a much lower dimension than α. We can then optimize θ with stochastic gradient descent, by approximating the empirical averaging 1 n n i=1 (·) in the objective with averages over small mini-batches; this would introduce biases in gradient estimation, but it is not an issue when the goal is only to get a reasonable approximation of the optimal ω.
2) Replacing kernel k Assume the kernel k(x,x) yields a random feature expansion of form k(x,x) = E β∼π [φ(x, β)φ(x, β)], where φ(x, β) is a feature map with parameter β and π is a distribution of β. We draw {β i } m i=1 i.i.d. from π, where m is taken to be much smaller than n. We replace k withk(x,x) = 1 m m i=1 φ(x, β i )φ(x, β i ) and let W to be the RKHS of kernelk. Then, we consider to solvê It is known that any function ω in W can be represented as ω( In this way, the problem reduces to optimizing an m-dimensional vector w, which can be (approximately) solved by standard convex optimization techniques.

F Concentration Inequality of General Functional Bellman Losses
Theorem 5.1 provides the concentration inequality of L W (q;D n ) when W is the unit ball of an RKHS. When W is a general function set, one can still obtain a general concentration bound using Radermacher complexity. DefineRq • W := {h(x, y) =Rq(x, y)ω(x) : ω ∈ W}. Using the standard derivation in Radermacher complexity theory in conjunction with Martingale theory (Rakhlin et al., 2015), we have where c q,W = sup ω∈W sup x,y (Rq(x, y) − Rq(x)) 2 ω(x) 2 and Rad(Rq • W) is the sequential Radermacher complexity (Rakhlin et al., 2015). A triangle inequality yields Therefore, When W equals the unit ball K of the RKHS of kernel k, we have c q,k = c q,W and hence this bound is strictly worse than the bound in Theorem 5.1.

G More on the Oracle Bound and its Dual Form
The oracle bound (28) provides another starting point for deriving optimization-based confidence bounds. We start with deriving the dual form of (28). Using Lagrangian multiplier, the optimization in Eq. (28) can be rewritten intoĴ where and ω is optimized in the set of all functions and serves as the Lagrangian multiplier here. Definê Then by the weak duality, we have J * Q,+ ≤F + Q, * (ω), ∀ω.
The derivation follows similarly for the lower bound. So for any ω ∈ W, we have Here the first two terms ofF + Q, * (ω) can be empirically estimated (it is the same as the first two terms of Eq. (36)), but the third term R(ω, q * ) depends on the unknown q * and hence need to be further upper bounded. Different upper bounds of R(ω, q * ) may yield different practical confidence intervals.
Our method can be viewed as constraining ω in the unit ball W of a normed function space W, and applying a worst case bound: for any ω ∈ W, we havê where we note that L W (q * ,D n ) = sup h∈W R(h, q * ) and the last step applies the high probability bound that With the same derivation on the lower bound counterpart, we have G.1 Proof of Proposition 6.2 Proof. Let Q null be the set of functions that are zero on {s i , s i } n i=1 , that is, Then we haveR (q * + g)(x i , y i ) =Rq * (x i , y i ), ∀i = 1, . . . , n. and Taking g(s, a) = zI(s / ∈ {s i , s i } n i=1 ), where z is any real number. Then we have Because Pr s∼D π,0 (s / ∈ {s i , s i } n i=1 ) = 0, we can take z to be arbitrary value to make E D π,0 [q * + g] to take arbitrary value.

H.1 Experimental Details
Environments and Dataset Construction We test our method on three environments: Inverted-Pendulum and CartPole from OpenAI Gym (Brockman et al., 2016), and a Type-1 Diabetes medical treatment simulator. For Inverted-Pendulum we discretize the action space to be {−1, −0.3, −0.2, 0, 0.2, 0.3, 1}. The action space of CartPole and the medical treatment simulator are both discrete.

Policy Construction
We follow a similar setup as Feng et al. (2020) to construct behavior and target policies. For all of the environments, we constraint our policy class to be a softmax policy and use PPO (Schulman et al., 2017) to train a good policy π, and we use different temperatures of the softmax policy to construct the target and behavior policies (we set the temperature τ = 0.1 for target policy and τ = 1 to get the behavior policy, and in this way the target policy is more deterministic than the behavior policy). We consider other choices of behavior policies in Section H.3.
For horizon lengths, We fix γ = 0.95 and set horizon length H = 50 for Inverted-Pendulum, H = 100 for CartPole, and H = 50 for Diabetes simulator.
Algorithm Settings We test the bound in Eq.(36)-(37). Throughout the experiment, we always set W = K, a unit ball of RKHS with kernel k(·, ·). We set Q = r QK , the zero-centered ball of radius r Q in an RKHS with kernelk(·, ·). We take both k andk to be Gaussian RBF kernel. The bandwidth of k andk are selected to make sure the function Bellman loss is not large on a validation set. The radius is selected to be sufficiently large to ensure that q * is included in Q. To ensure a sufficiently large radius, we use the data to approximate aq so that its functional Bellman loss is small than n . Then we set r Q = 10 * q K . We optimize ω using the random feature approximation method described in Appendix E. Once ω + and ω − are found, we evaluate the bound in Eq. (36) exactly, to ensure the theoretical guarantee holds.

H.2 Sensitivity to Hyper-Parameters
We investigate the sensitivity of our algorithm to the choice of hyper-parameters. The hyper-parameter mainly depends on how we choose our function class Q and W.
Radius of Q Recall that we choose Q to be a ball in RKHS with radius r Q , that is, whereK is the unit ball of the RKHS with kernelk. Ideally, we want to ensure that r Q ≥ q * K so that q * ∈ Q.
Since it is hard to analyze the behavior of the algorithm when q * is unknown, we consider a synthetic environment where the true q * is known. This is done by explicitly specifying a q * insideK and then infer the corresponding deterministic reward function r(x) by inverting the Bellman equation: r(x) := q * (x) − γE x ∼Pπ(·|x) [q * (x )].
Here r is a deterministic function, instead of a random variable, with an abuse of notation. In this way, we can get access to the true RKHS norm of q * : For simplicity, we set both the state space S and action space A to be R and set a Gaussian policy π(a|s) ∝ exp(f (s, a)/τ ), where τ is a positive temperature parameter. We set τ = 0.1 as target policy and τ = 1 as behavior policy. Figure 3 shows the results as we set r Q to be ρ * , 10ρ * and 100ρ * , respectively. We can see that the tightness of the bound is affected significantly by the radius when the number n of samples is very small. However, as the number n of samples grow (e.g., n ≥ 2 × 10 3 in our experiment), the length of the bounds become less sensitive to the changing of the predefined norm of Q.

Similarity Between Behavior Policy and Target Policy
We study the performance of changing temperature of the behavior policy. We test on Inverted-Pendulum environment as previous described. Not surprisingly, we can see that the closer the behavior policy to the target policy (with temperature τ = 0.1), the tighter our confidence interval will be, which is observed in Figure 4(a).
Bandwidth of RBF kernels We study the results as we change the bandwidth in kernel k andk for W and Q, respectively. Figure 4(b) shows the length of the confidence interval when we use different bandwidth pairs in the Inverted-Pendulum environment. We can see that we get relatively tight confidence bounds as long as we set the bandwidth in a reasonable region (e.g., we set the bandwidth of k in [0.1, 0.5], the bandwidth ofk in [0.5, 3]).

H.3 Sensitivity to the Data Collection Procedure
We investigate the sensitivity of our method as we use different behavior policies to collect the datasetD n .
Varying Behavior Policies We study the effect of using different behavior policies. We consider the following cases: 1. Data is collected from a single behavior policy of form π α = απ + (1 − α)π 0 , where π is the target policy and π 0 is another policy. We construct π and π 0 to be Gaussian policies of form π(a|s) ∝ exp(f (s, a)/τ ) with different temperature τ , where temperature for target policy is τ = 0.1 and temperature for π 0 is τ = 1.
We show in Figure 5(a) that the length of the confidence intervals by our method as we vary the number n of transition pairs and the mixture rate α. We can see that the length of the interval decays with the sample size n for all mixture rate α. Larger α yields better performance because the behavior policies are closer to the target policy.
Varying Trajectory Length T inD n As we collectD n , we can either have a small number of long trajectories, or a larger number of short trajectories. In Figure 5(b)-(c), we vary the length T of the trajectories as we collectD n , while fixing the total number n of transition pairs. In this way, the number of trajectories in eachD n would be m = n/T . We can see that the trajectory length does not impact the results significantly, especially when the length is reasonably large (e.g., T ≥ 20).

I Finite Horzion and Time-varying MDPs
In this paper we mainly focus on infinite-horizion and time-homogeneous MDPs. Here we consider how to extend our method to finite-horizon and time-inhomogeneous MDPs, where the transition probability and reward are time-dependent and we have a finite horizon length H < +∞. That is, we consider cases when the next state and reward follow an unknown, time-dependent transition distribution, (r t , s t+1 ) ∼ P(·|s t , a t ; t), and we want to estimate the following finite-horizon expected reward J π,P,H = E π,P H−1 t=0 γ t r t | s 0 ∼ D 0 , where the horizon length H is finite but can be large. We should distinguish H (which is a part of the problem definition) with the trajectory length T of the data inD n (which is a part of the data collection procedure).
We can approach this problem by transforming it into an infinite-horizon and time-homogeneous problem and then apply our method. This can be done by incoporating the time t as a part of the state vector. To be concrete, assume we collect a set of transition pairs (s i , a i , r i , s i , t i ) n i=1 , where in addition to (s i , a i , r i , s i ), we also record the time t i when the transition occurs. We define an augmented states i = [s i , t i ] and s i = [s i , t i + 1] which include the time t i as a part. Our method can be then applied without modification on the augmented dataset {s i , a i , r i } n i=1 , once we ensure that q([s, t], a) = 0 for all t ≥ H when defining the function set Q. To build an RKHS Q that satisfy q([s, t], a) = 0, ∀t ≥ H, q ∈ Q, we just need to use a kernelk ([s, a, t], [s , a , t ]) such thatk([s, a, t], [s , a , t ]) = 0 whenever t ≥ H or t ≥ H.