Year: 2021

Abstract:
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 nonasymptotic 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.

heading: 1 INTRODUCTION
text: Off-policy evaluation (OPE) seeks to estimate the expected reward of a target policy in reinforcement learnings (RL) from observational data collected under different policies (e.g., Murphy et al., 2001; Fonteneau et al., 2013; Jiang & 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 medicine, self-driving, where interactive “on-policy” data is expensive or even infeasible to collect. A critical challenge in OPE is the uncertainty estimation, as having reliable confidence bounds is essential for making high-stakes decisions. In this work, we aim to tackle this problem by providing non-asymptotic confidence intervals of the expected value of the target policy. Our method allows us to rigorously quantify the uncertainty of the prediction and hence avoid the dangerous case of being overconfident in making costly and/or irreversible decisions.
However, off-policy evaluation per se has remained a key technical challenge in the literature (e.g., Precup, 2000; Thomas & Brunskill, 2016; Jiang & 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, and 2) the data is collected under arbitrary and unknown algorithms (a.k.a. behavior-agnostic). As a consequence, the collected data can exhibit arbitrary dependency structure, which makes constructing rigorous non-asymptotic confidence bounds particularly challenging. Traditionally, the only approach to provide non-asymptotic confidence bounds in OPE is to combine importance sampling (IS) with concentration inequalities (e.g., Thomas et al., 2015a;b), which, however, tends to degenerate for long/infinite horizon problems (Liu et al., 2018a). Furthermore,
∗Equal contribution.
neither can this approach be applied to the behavior-agnostic settings, nor can it effectively handle the complicated time dependency structure inside individual trajectories. Instead, it requires to use a large number of independently collected trajectories drawn under known policies.
In this work, we provide a practical approach for Behavior-agnostic, Off-policy, Infinite-horizon, Non-asymptotic, Confidence intervals based on arbitrarily Dependent data (BONDIC). Our method is motivated by a recently proposed optimization-based (or variational) approach to estimating OPE confidence bounds (Feng et al., 2020), which leverages a tail bound of kernel Bellman statistics (Feng et al., 2019). Our approach achieves a new bound that is both an order-of-magnitude tighter and computationally efficient than that of Feng et al. (2020). Our improvements are based on two pillars 1) developing a new primal-dual perspective on the non-asymptotic OPE confidence bounds, which is connected to a body of recent works on infinite-horizon value estimation (Liu et al., 2018a; Nachum et al., 2019a; Tang et al., 2020a; Mousavi et al., 2020); and 2) offering a new tight concentration inequality on the kernel Bellman statistics that applies to behavior-agnostic off-policy data with arbitrary dependency between transition pairs. Empirically, we demonstrate that our method can provide reliable and tight bounds on a variety of well-established benchmarks.
Related Work Besides the aforementioned approach based on the combination of IS and concentration inequalities (e.g., Thomas et al., 2015a), bootstrapping methods have also been widely used in off-policy estimation (e.g., White & White, 2010; Hanna et al., 2017; Kostrikov & Nachum, 2020). But the latter is limited to asymptotic bounds. Alternatively, Bayesian methods (e.g. Engel et al., 2005; Ghavamzadeh et al., 2016a) offers a different way to estimate the uncertainty in RL, but fails to guarantee frequentist coverage. In addition, Distributed RL (Bellemare et al., 2017) seeks to quantify the intrinsic uncertainties inside the Markov decision process, which is orthogonal to the estimation of uncertainty that we consider.
Our work is built upon the recent advances in behavior-agnostic infinite-horizon OPE, including Liu et al. (2018a); Feng et al. (2019); Tang et al. (2020a); Mousavi et al. (2020), as well as the DICE-family (e.g., Nachum et al., 2019a; Zhang et al., 2020a; Yang et al., 2020b). In particular, our method can be viewed as extending the minimax framework of the infinite-horizon OPE in the infinite data region by Tang et al. (2020a); Uehara et al. (2020); Jiang & Huang (2020) to the non-asymptotic finite sample region.
Outline For the rest of the paper, we start with the problem statement in Section 2 , and an overview on the two dual approaches to infinite-horizon OPE that are tightly connected to our method in Section 3. We then present our main approach in Section 4 and perform empirical studies in Section 5. The proof and an abundance of additional discussions can be found in Appendix.

heading: 2 BACKGROUND, DATA ASSUMPTION, PROBLEM SETTING
text: Consider an agent acting in an unknown environment. At each time step t, the agent observes the current state st in a state space S, takes an action at ∼ π(· | st) in an action space A according to a given policy π; then, the agent receives a reward rt and the state transits to s′t = st+1, following an unknown transition/reward distribution (rt, st+1) ∼ P(· | st, at). Assume the initial state s0 is drawn from an known initial distribution D0. Let γ ∈ (0, 1) be a discount factor. In this setting, the expected reward of π is defined as Jπ := Eπ [∑T t=0 γ trt | s0 ∼ D0 ] , which is the expected total discounted rewards when we execute π starting from D0 for T steps. In this work, we consider the infinite-horizon case with T → +∞. Our goal is to provide an interval estimation of Jπ for a general and challenging setting with significantly released constraints on the data. In particular, we assume the data is behavior-agnostic and off-policy, which means that the data can be collected from multiple experiments, each of which can execute a mix of arbitrary, unknown policies, or even follow a non-fixed policy. More concretely, suppose that the model P is unknown, and we have a set of transition pairs D̂n = (si, ai, ri, s′i) n i=1 collected from previous experiments in a sequential order, such that for each data point i, the (ri, s′i) is drawn from the model P(· | si, ai), while (si, ai) is generated with an arbitrary black box given the previous data points. We formalize both the data assumption and goal as below.
Assumption 2.1 (Data Assumption). Assume the data D̂n = (si, ai, ri, s′i)ni=1 is drawn from an arbitrary joint distribution, such that for each i = 1, . . . , n, conditional on D̂<i := (sj , aj , rj , s′j)j<i ∪ (si, ai), the subsequent local reward and next state (ri, s′i) are drawn from P(· | si, ai).
Goal Given a confidence level δ ∈ (0, 1), we want to construct an interval [Ĵ−, Ĵ+] ⊂ R based on the data D̂n, such that Pr(Jπ ∈ [Ĵ−, Ĵ+]) ≥ 1− δ, where Pr(·) is w.r.t. the randomness of the data. The partial ordering on the data points is introduced to accommodate the case that si+1 equals s′j for some j ≤ i. The data assumption only requires that (ri, s′i) is generated from P(· | si, ai), and imposes no constraints on how (si, ai) is generated. This provides great flexibility in terms of the data collection process. In particular, we do not require (si, ai)ni=1 to be independent as always assumed in recent works (Liu et al., 2018a; Mousavi et al., 2020).
A crucial fact is that our data assumption actually implies a martingale structure on the empirical Bellman residual operator of the Q-function, As we will show in Section 4.1, this enables us to derive a key concentration inequality underpinning our non-asymptotic confidence bounds.
Here, we summarize a few notations that will simplify the presentation in the rest of work. First of all, we append each (si, ai, ri, s′i) with an action a ′ i ∼ π(· | s′i) following s′i. This can be done for free as long as π is given (See the Remark in Section 3). Also, we write xi = (si, ai), x′i = (s ′ i, a ′ i), and yi = (x′i, ri) = (s ′ i, a ′ i, ri). Correspondingly, define X = S ×A to be the state-action space and Y = X ×R. Denote Pπ(y | x) = P(s′, r | x)π(a′ | s′). In this way, the observed data can be written as pairs of {xi, yi}ni=1, and Assumption 2.1 is equivalent to saying that yi ∼ Pπ(· | xi) given D̂<i, which is similar to a supervised learning setting. We equalize the data D̂n with its empirical measure D̂n = ∑n i=1 δxi,yi/n, where δ is the Delta measure.

heading: 3 TWO DUAL APPROACHES TO INFINITE-HORIZON OFF-POLICY ESTIMATION
text: The deficiency of the traditional IS methods on long-/infinite-horizon RL problems (a.k.a. the curse of horizon (Liu et al., 2018a)) has motivated a line of work on developing efficient infinite-horizon value estimation (e.g., Liu et al., 2018a; Feng et al., 2019; Nachum et al., 2019a; Zhang et al., 2020a; Mousavi et al., 2020; Tang et al., 2020a). The main idea is to transform the value estimation problem into estimating either the Q-function or the visitation distribution (or its related density ratio) of the policy π. This section introduces and reinterprets these two tightly connected methods, which serves to lay out a foundation for our main confidence bounds from a primal and dual perspective.
Given a policy π, its Q-function is defined as qπ(x) = Eπ [ ∑∞ t=0 γ
trt | x0 = x], where the expectation is taken when we execute π initialized from a fixed state-action pair (s0, a0) = x0 = x. Let Dπ,t be the distribution of (xt, yt) = (st, at, s′t, a ′ t, rt) when executing policy π starting from s0 ∼ D0 for
t steps. The visitation distribution of π is defined as Dπ = ∑∞ t=0 γ
tDπ,t. Note that Dπ integrates to 1/(1− γ), while we treat it as a probability measure in the notation. The expected reward Jπcan be expressed using either qπ or Dπ as follows:
Jπ := Eπ [ ∞∑ t=0 γtrt ] = Er∼Dπ [r] = Ex∼Dπ,0 [qπ(x)], (1)
where r ∼ Dπ (resp. x ∼ Dπ,0) denotes sampling from the r-(resp. x-) marginal distribution of Dπ (resp. Dπ,0). Eq. (1) plays a key role in the infinite-horizon value estimation by transforming the estimation of Jπ into estimating either qπ or Dπ .
Value Estimation via Q Function Because Dπ,0(x) = D0(s)π(a|s) is known, we can estimate Jπ by Ex∼Dπ,0 [q̂(x)] with any estimation q̂ 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:
R̂q(x, y) = q(x)− γq(x′)− r, Rπq(x) = Ey∼Pπ(·|x) [ R̂q(x, y) ] . (2)
It is well-known that qπ is the unique solution of the Bellman equation Rπq = 0. Since yi ∼ Pπ(·|xi) for each data point in D̂n, if q = qπ , then R̂q(xi, yi), i = 1, . . . , n are all zero-mean random variables.
Let ω be any function from X to R, then∑i R̂q(xi, yi)ω(xi) also has zero mean. This motivates the following functional Bellman loss (Feng et al., 2019; 2020; Xie & Jiang, 2020),
LW(q; D̂n) := sup ω∈W
{ 1
n n∑ i=1 R̂q(xi, yi)ω(xi)
} , (3)
whereW is a set of functions ω : X → R. To ensure that the sup is finite,W is typically set to be an unit ball of some normed function spaceWo, such thatW = {ω ∈ Wo : ‖ω‖Wo ≤ 1}. Feng et al. (2019) considers the simple case whenW is taken to be the unit ball K of the reproducing kernel Hilbert space (RKHS) with a positive definite kernel k : X × X → R, in which case the loss has a simple closed form solution:
LK(q; D̂n) = √√√√ 1 n2 n∑ ij=1 R̂q(xi, yi)k(xi, xj)R̂q(xj , yj). (4)
Note that the RHS of Eq. (4) 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 data distribution D̂n covers the whole space (which may require an infinite data size) and k is an integrally strictly positive definite kernel, LK(q; D̂n) = 0 iff q = qπ . Therefore, one can estimate qπ by minimizing LK(q, D̂n).
Remark The empirical Bellman residual operator R̂ can be extended to R̂q(x, y) = q(x)− r − γ 1m ∑m `=1 q(s
′, a′`), where {a′`}mi=1 are i.i.d. drawn from π(·|s′). As m increases, this gives a lower variance estimation of Rπq. If m = +∞, we have R̂q(x, y) = q(x)− r − γEa′∼π(· | s′)[q(s′, a′)], which coincides with the operator used in the expected SARSA (Sutton & Barto, 1998). In fact, without any modification, all results in this work can be applied to R̂q for any m.
Value Estimation via Visitation Distribution Another way to estimate Jπ in Eq. (1) 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 ω(xi) to each data point xi in D̂n. We can choose the function ω : X → R properly such that Dπ and hence Jπ can be approximated by the ω-weighted empirical measure of D̂n as follows:
Jπ ≈ Ĵω := ED̂ωn [r] = 1
n n∑ i=1 ω(xi)ri, Dπ ≈ D̂ωn := 1 n n∑ i=1 ω(xi)δxi,yi . (5)
Intuitively, ω can be viewed as the density ratio between Dπ and D̂n, although the empirical measure D̂n may not have well-defined density. Liu et al. (2018a); Mousavi et al. (2020) proposed to estimate ω by minimizing a discrepancy measure between D̂ωn and Dπ. To see this, note that D = Dπ if and only if ∆(D, q) = 0 for any function q, where
∆(D, q) = ED[γq(x′)− q(x)]− EDπ [γq(x′)− q(x)] = ED[γq(x′)− q(x)] + EDπ,0 [q(x)], (6)
using the fact that EDπ [γq(x′) − q(x)] = −EDπ,0 [q(x)] (Theorem 1, Liu et al., 2018a). Also note that the RHS of Eq. (6) can be practically calculated given any D and q without knowing Dπ . Let Q be a set of functions q : X → R. One can define the following loss for ω:
IQ(ω; D̂n) = sup q∈Q
{ ∆(D̂ωn , q) } . (7)
Similar to LW(q; D̂n), when Q is a ball in RKHS, IQ(ω; D̂n) also has a bilinear closed form analogous to Eq. (4); see Mousavi et al. (2020) and Appendix F. As we show in Section 4, IQ(ω; D̂n) and LW(q; D̂n) are connected to the primal and dual views of our confidence bounds, respectively.

heading: 4 MAIN APPROACH
text: Let Q be a large enough function set including the true Q-function qπ, that is, qπ ∈ Q. Following Feng et al. (2020), a confidence interval [Ĵ−Q,W , Ĵ + Q,W ] of Jπ can be constructed as follows:
Ĵ+Q,W = sup q∈Q
{ EDπ,0 [q] s.t. LW(q; D̂n) ≤ εn } , (8)
and Ĵ−Q,W is defined in a similar way by replacing sup on q ∈ Q with inf . The idea here is to seek the extreme q function with the largest (resp. smallest) expected values in set F := Q ∩ {q : LK(q; D̂n) ≤ εn}. Therefore, Eq. (8) would be a 1− δ confidence interval if qπ is included in F with at least probability 1− δ, which is ensured when qπ ∈ Q and
Pr(LW(qπ; D̂n) ≤ εn) ≥ 1− δ . (9) Feng et al. (2020) showed that in the RKHS case whenW = K, Eq. (9) can be achieved with
εn = √√√√2cqπ,k ( n− 1 n √ log(1/δ) n + 1 n ) , (10)
when n is an even number, where cqπ,k = supx,y R̂qπ(x, y) 2k(x, x). This was proved using Hoeffding’s inequality for U-statistics (Hoeffding, 1963). To solve Eq. (8) efficiently, Feng et al. (2020) took Q to be a ball in RKHS with random feature approximation. Unfortunately, this method as described by Eq. (8)-(10) has two major disadvantages:
1) Bound Needs to Be Tightened (Section 4.1) The bound of εn = O(n−1/4) in Eq. (10) is sub-optimal in rate. In Section 4.1, we improve it by an εn = O(n−1/2) bound under the mild Assumption 2.1, which gets rid of the independence requirement between the transition pairs. Our tightened bound is achieved by firstly noting a Martingale structure on the empirical Bellman operator under Assumption 2.1, and then applying an inequality in Pinelis (1992).
2) Dependence on Global Optimization (Section 4.2) The bound in Eq. (8) is guaranteed to be a 1− δ confidence bound only when the maximization in Eq. (8) is solved to global optimality. With a large n, this leads to a high computational cost, even when choosing Q as the RKHS. Feng et al. (2020) solved Eq. (8) approximately using a random feature technique, but this method suffers from a gap between the theory and practice. In Section 4.2, we address this problem by presenting a dual form of Eq. (8), which sidesteps solving the challenging global optimization in Eq. (8). Moreover, the dual form enables us to better analyze the tightness of the confidence interval and issues regarding the choices of Q andW .

heading: 4.1 A TIGHTER CONCENTRATION INEQUALITY
text: In this section, we explain our method to improve the bound in Eq. (10) by giving a tighter concentration inequality for the kernel Bellman loss in Eq. (4). We introduce the following semi-expected kernel Bellman loss:
L∗K(q; D̂n) = √√√√ 1 n2 n∑ ij=1 Rπq(xi)k(xi, xj)Rπq(xj) , (11)
in which we replace the empirical Bellman residual operator R̂q in Eq. (3) with its expected counterpart Rπq, but still take the empirical average over {xi}ni=1 in D̂n. For a more general function set W , we can similarly define L∗W(q; D̂n) by replacing R̂q with Rπq in Eq. (3). Obviously, we have L∗W(q; D̂n) = 0 when q = qπ .
Theorem 4.1 below shows that LK(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 is able to hold even without any independence assumption between {xi}. An easy way to make sense of this is by recognizing that the randomness in yi conditional on xi is aggregated through averaging, even when {xi} are deterministic. As Assumption 2.1 does not impose any (weak) independence between {xi}, we cannot establish that LK(q; D̂n) concentrates around its mean ED̂n [LK(q; D̂n)] (which is a full expected kernel bellman loss), without introducing further assumptions.
Theorem 4.1. Assume K is the unit ball of RKHS with a positive definite kernel k(·, ·). Let cq,k := supx∈X ,y∈Y(R̂q(x, y)−Rπq(x))2k(x, x) <∞. Under Assumption 2.1, for any δ ∈ (0, 1), with at
least probability 1− δ, we have∣∣∣LK(q; D̂n)− L∗K(q; D̂n)∣∣∣ ≤ √ 2cq,k log(2/δ)
n . (12)
In particular, when q = qπ , we have cqπ,k = supx,y(R̂qπ(x, y)) 2k(x, x), and LK(qπ; D̂n) ≤ √ 2cqπ,k log(2/δ)
n . (13)
Intuitively, to see why we can expect an O(n−1/2) bound, note that LK(q, D̂n) consists of the square root of the product of two R̂q terms, each of which contributes an O(n−1/2) error w.r.t. Rπq.
Technically, the proof is based on a key observation: Assumption 2.1 ensures that Zi := R̂q(xi, yi)− Rπq(xi), i = 1, . . . , n forms a martingale difference sequence w.r.t. {D̂<i : ∀i = 1, . . . , n}, in the sense that E[Zi | D̂<i] = 0, ∀i. See Appendix B for details. The proof also leverages a special property of RKHS and applies a Hoeffding-like inequality on the Hilbert spaces as in Pinelis (1992) (see Appendix B). For other more general function setsW , we establish in Appendix E a similar bound by using Rademacher complexity, although it yields a less tight bound than Eq. (12) when W = K.

heading: 4.2 DUAL CONFIDENCE BOUNDS
text: We derive a dual form of Eq. (8) that sidesteps the need for solving the challenging global optimization in Eq. (8). To do so, let us plug the definition of LW(q; D̂n) into Eq. (3) and introduce a Lagrange multiplier:
Ĵ+Q,W = sup q∈Q inf h∈W inf λ≥0
EDπ,0 [q]− λ ( 1
n n∑ i=1
h(xi)R̂q(xi, yi)− εn )
(14)
= sup q∈Q inf ω∈Wo
{ EDπ,0 [q]− 1
n n∑ i=1 ω(xi)R̂q(xi) + εn ‖ω‖Wo
} , (15)
where we use ω(x) = λh(x). Exchanging the order of min/max and some further derivation yields the following main result.
Theorem 4.2. I) LetW be the unit ball of a normed function spaceWo. We have Ĵ+Q,W ≤ F̂+Q (ω) := ED̂ωn [r] + IQ(ω; D̂n) + εn ‖ω‖Wo , ∀ω ∈ Wo , Ĵ−Q,W ≥ F̂−Q (ω) := ED̂ωn [r]− I−Q(ω; D̂n)− εn ‖ω‖Wo , ∀ω ∈ Wo ,
(16)
where −Q = {−q : q ∈ Q} and hence I−Q(ω; D̂n) = IQ(ω; D̂n) if Q = −Q. Further, we have Ĵ+Q,W = infω∈Wo F̂ + Q (ω) and Ĵ − Q,W = supω∈Wo F̂ − Q (ω) if Q is convex and there exists a q ∈ Q that satisfies the strict feasibility condition that LW(q; D̂n) < εn.
II) For D̂n and δ ∈ (0, 1), assume Wo and εn ∈ R satisfy Eq. (9) (e.g., via Theorem 4.1). Then for any function set Q with qπ ∈ Q, and any function ω+, ω− ∈ Wo (the choice of Q, ω+, ω− can depend on D̂n arbitrarily), we have
Pr ( Jπ ∈ [ F̂−Q (ω−), F̂ + Q (ω+) ]) ≥ 1− δ . (17)
Theorem 4.2 transforms the original bound in Eq. (8), framed in terms of q and LW(q; D̂n), into a form that involves the density-ratio ω and the related loss IQ(ω; D̂n). The bounds in Eq. (16) can be interpreted as assigning an error bar around the ω-based estimator Ĵω = ED̂ωn [r] in Eq. (5), with the error bar of I±Q(ω; D̂n) + εn ‖ω‖Wo . Specifically, the first term I±Q(ω; D̂n) measures the discrepancy between D̂ωn and Dπ as discussed in Eq. (7), whereas the second term captures the randomness in the empirical Bellman residual operator R̂qπ .
Compared with Eq. (8), the global maximization on q ∈ Q is now transformed inside the IQ(ω; D̂n) term, which yields a simple closed form solution in the RKHS case (see Appendix F). In practice, we can optimize ω+ and ω− to obtain the tightest possible bound (and hence recover the primal bound) by minimizing/maximizing F̂+Q (ω) and F̂ − Q (ω), but it is not necessary to solve the optimization to global optimality. WhenWo is an RKHS, by the standard finite representer theorem (Scholkopf & Smola, 2018), the optimization on ω reduces to a finite dimensional optimization, which can be sped up with any favourable approximation techniques. We elaborate on this in Appendix D.
Length of the Confidence Interval The form in Eq. (16) 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
length([Ĵ−Q,W , Ĵ + Q,W ]) = inf
ω∈Wo
{ 2IQ(ω; D̂n) + 2εn ‖ω‖Wo } .
Given εn is O(n−1/2), we can make the overall length of the optimal confidence interval also O(n−1/2), as long asWo is rich enough to include a good density ratio estimator ω∗ that satisfies IQ(ω
∗; D̂n) = O(n −1/2) and has a bounded norm ‖ω∗‖Wo .
We can expect to achieve IQ(ω∗; D̂n) = O(n−1/2), when (1) Q has an O(n−1/2) sequential Rademacher complexity (Rakhlin et al., 2015) (e.g., a finite ball in RKHS); and (2) D̂n is collected following a Markov chain with strong mixing condition and weakly converges to some limit distribution D∞ whose support is X , and therefore we can define ω∗ as the density ratio between Dπ and D∞. Refer to Appendix C 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.
Choice ofW and Q To ensure the concentration inequality in Theorem 4.1 is valid, the choice of Wo cannot depend on the data D̂n. Therefore, we should use a separate holdout data to construct a data-dependentWo. In contrast, the choice of Q can depend on the data D̂n arbitrarily, since it is a part of the optimization bound Eq. (8) but not in the tail bound Eq. (9). In this light, one can construct the best possible Q by exploiting the data information in the most favourable way. For example, we can construct an estimator of q̂ ≈ qπ based on any state-of-the-art method (e.g., Q-learning or model-based methods), and set Q to be a ball centering around q̂ such that qπ − q̂ ∈ Q. This enables post-hoc analysis based on prior information on qπ , as suggested in Feng et al. (2020).
Mis-specification of Q and Oracle Upper/Lower Estimates Our result relies on the assumption that qπ ∈ Q. However, as with other statistical estimation problems, there exists no provably way to empirically verify the correctness of model assumptions such as qπ ∈ Q. Because empirical data only reveals the information of the unknown function (in our case qπ) on a finite number data points, but no conclusion can be made on the unseeing data points without imposing certain smoothness assumption. Typically, what we can do is the opposite: reject qπ ∈ Q when the Bellman loss LW(q; D̂n) of all q in Q is larger than the threshold εn. We highlight that, even without verifying qπ ∈ Q, our method can still be viewed as a confidence interval of a best possible (oracle) upper and lower estimation given the data D̂n plus the assumption that qπ ∈ Q, defined as
Ĵ+Q,∗ = sup q∈Q
{ EDπ,0 [q] s.t. R̂q(xi, yi) = R̂qπ(xi, yi), ∀i = 1, . . . , n } . (18)
In fact, it is impossible to derive empirical upper bounds lower than Ĵ+Q,∗, as there is no way to distinguish q and qπ if R̂q(xi, yi) = R̂qπ(xi, yi) for all i. But our interval [ĴQ,K, Ĵ+Q,K] provides a 1− δ confidence outer bound of [Ĵ−Q,∗, Ĵ+Q,∗] once Eq. (9) holds, regardless if qπ ∈ Q holds or not. Hence, it is of independent interest to further explore the dual form of Eq. (18), which is another starting point for deriving our bound. We have more discussion in Appendix G.
Lastly, we argue that it is important to include the Q in the bound. Proposition G.1 in Appendix shows that removing the q ∈ Q constraint in Eq. (18) would lead to an infinite upper bound,
unless the {si, s′i}ni=1 from D̂n almost surely covers the whole state space S in the sense that Prs∼D0(s ∈ {si, s′i}ni=1) = 1.

heading: 5 EXPERIMENTS
text: We compare our method with a variety of existing algorithms for obtaining asymptotic and nonasymptotic 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: InvertedPendulum 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 our paper. Throughout the experiment, we always setW = K, the unit ball of the RKHS with positive definite kernel k, and set Q = rQK̃, the ball of radius rQ in the RKHS with another kernel k̃. We take both kernels to be Gaussian RBF kernel and choose rQ and the bandwidths of k and k̃ 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 D. Once ω is found, we evaluate the bound in Eq. (16) exactly to ensure that the theoretical guarantee holds.
Baseline Algorithms We compare our method with four existing baselines, including the IS-based non-asymptotic 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 & Nachum (2020), and the bound in Feng et al. (2020) which is equivalent to the primal bound in (8) 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
1 https://github.com/jxx123/simglucose.
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 rQ and coefficient cqπ,k in Eq. (13) (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.

heading: 6 CONCLUSION
text: We develop a dual approach to construct high confidence bounds for off-policy evaluation with an improved rate over Feng et al. (2020). Our method can handle dependent data, and does not require a global optimization to get a valid bound. Empirical results demonstrate that our bounds is tight and valid compared with a range of existing baseline. Future directions include leveraging our bounds for policy optimization and safe exploration.

heading: A PROOF OF THE DUAL BOUND IN THEOREM 4.2
text: Proof. Introducing a Lagrange multiplier, the bound in (8) is equivalent to
Ĵ+Q,W = max q∈Q min λ≥0
{ EDπ,0 [q] − λ ( max h∈W 1 n n∑ i=1 h(xi)R̂q(xi, yi)− εn )}
= max q∈Q min λ≥0 min h∈W
{ EDπ,0 [q] − λ ( 1
n n∑ i=1
h(xi)R̂q(xi, yi)− εn )}
= max q∈Q min ω∈Wo
{ EDπ,0 [q] − 1
n n∑ i=1 ω(xi)R̂q(xi, yi) + εn ‖ω‖Wo
} ,
where we use ω = λh(x), such that λ is replaced by ‖ω‖Wo . Define
M(q, ω; D̂n) = EDπ,0 [q] − 1
n n∑ i=1 ω(xi)R̂q(xi, yi) + εn ‖ω‖Wo
= ED̂ωn [r] + ∆(D̂ ω n , q) + εn ‖ω‖Wo .
Then we have
max q∈Q M(q, ω; D̂n) = ED̂ωn [r] + maxq∈Q ∆(D̂ ω n , q) + εn ‖ω‖Wo
= ED̂ωn [r] + IQ(ω; D̂n) + εn ‖ω‖Wo = F̂+Q (ω).
Therefore,
Ĵ+Q,W = max q∈Q min ω∈Wo M(q, ω; D̂n)
≤ min ω∈Wo max q∈Q M(q, ω; D̂n)
= min ω∈Wo
F̂+Q (ω).
The lower bound follows analogously. The strong duality holds when the Slater’s condition is satisfied (Nesterov, 2013), which amounts to saying that the primal problem in (8) is convex and strictly feasible; this requires that Q is convex and there exists at least one solution q ∈ Q that satisfy that constraint strictly, that is, LW(q; D̂n) < εn; note that the objective function Q is linear on q and the constraint function LW(q; D̂n) is always convex on q (since it is the sup a set of linear functions on q following (3)).

heading: B PROOF OF CONCENTRATION BOUND IN THEOREM 4.1
text: Our proof require the following Hoeffding inequality on Hilbert spaces by Pinelis (Theorem 3, 1992); see also Section 2.4 of Rosasco et al. (2010).
Lemma B.1. (Theorem 3, Pinelis, 1992) Let H be a Hilbert space and {fi}ni=1 is a Martingale sequence inH that satisfies supi ‖fi‖H ≤ σ almost surely. We have for any > 0,
Pr (∥∥∥∥∥ 1n n∑ i=1 fi ∥∥∥∥∥ H ≥ ) ≤ 2 exp ( −n 2 2σ2 ) .
Therefore, with probability at least 1− δ, we have ∥∥ 1 n ∑n i=1 fi ∥∥ H ≤ √ 2σ2 log(2/δ) n .
Lemma B.2. Let k(x, x′) be a positive definite kernel whose RKHS isHk. Define fi(·) = R̂q(xi, yi)k(xi, ·)−Rπq(xi)k(xi, ·).
Assume Assumption 2.1 holds, then {fi}ni=1 is a Martingale difference sequence inHk w.r.t. T<i := (xj , yj)j<i ∪ (xi). That is, E [fi+1(·) | T<i] = 0. In addition,∥∥∥∥∥ 1n n∑ i=1 fi ∥∥∥∥∥ 2
Hk
= 1
n2 n∑ ij=1 ( R̂q(xi, yi)−Rπq(xi) ) k(xi, xj) ( R̂q(xj , yj)−Rπq(xj) ) ,
and ‖fi‖2Hk ≤ cq,k for ∀i = 1, . . . , n.
Proof of Theorem 4.1. Following Lemma B.1 and Lemma B.2, since {fi}ni=1 is a Martingale difference sequence inHk with ‖fi‖Hk ≤ cq,k almost surely, we have with probability at least 1− δ,
1
n2 n∑ ij=1 ( R̂q(xi, yi)−Rπq(xi) ) k(xi, xj) ( R̂q(xj , yj)−Rπq(xj) ) = ∥∥∥∥∥ 1n n∑ i=1 fi ∥∥∥∥∥ 2
Hk
≤ 2cq,k log(2/δ) n .
Using Lemma B.3 below, we have∣∣∣LK(q; D̂n)− L∗K(q; D̂n)∣∣∣ ≤ ∥∥∥∥∥ 1n n∑ i=1 fi ∥∥∥∥∥ Hk ≤ √ 2cq,k log(2/δ) n .
This completes the proof.
Lemma B.3. Assume k(x, x′) is a positive definite kernel. We have∣∣∣LK(q; D̂n)− L∗K(q; D̂n)∣∣∣2 ≤ 1n2 n∑
ij=1
( R̂q(xi, yi)−Rπq(xi) ) k(xi, xj) ( R̂q(xj , yj)−Rπq(xj) ) .
Proof. Define
ĝ(·) = 1 n n∑ i=1 R̂q(xi, yi)k(xi, ·) g(·) = 1 n n∑ i=1 Rπq(xi)k(xi, ·).
Then we have
‖ĝ‖2Hk = 1
n2 n∑ ij=1 R̂q(xi, yi)k(xi, xj)R̂q(xj , yj) = L̂K(q; D̂n),
‖g‖2Hk = 1
n2 n∑ ij=1 Rπq(xi)k(xi, xj)Rπq(xj) = L ∗ K(q; D̂n),
‖ĝ − g‖2Hk = 1
n2 n∑ ij=1 ( R̂q(xi, yi)−Rπq(xi) ) k(xi, xj) ( R̂q(xj , yj)−Rπq(xj) ) .
The result then follows the triangle inequality ∣∣‖ĝ‖Hk − ‖g‖Hk ∣∣ ≤ ‖ĝ − g‖Hk .
B.1 CALCULATION OF cqπ,k
The practical calculation of the coefficient cqπ,k in the concentration inequality was discussed in Feng et al. (2020), which we include here for completeness.
Lemma B.4. (Feng et al. (2020) Lemma 3.1) Assume the reward function and kernel function is bounded with supx |r(x)| ≤ rmax and supx,x′ |k(x, x′)| ≤ Kmax, we have:
cqπ,k := sup x∈X ,y∈Y
(R̂qπ(x, y)) 2k(x, x) ≤ 4Kmaxr
2 max
(1− γ)2 .
In practice, we get access to Kmax from the kernel function that we choose (e.g., Kmax = 1 for RBF kernels), and rmax from the knowledge on the environment.

heading: C MORE ON THE TIGHTNESS OF THE CONFIDENCE INTERVAL
text: 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 ifWo is chosen to be sufficiently rich so that it includes a ω ∈ Wo such that D̂ωn ≈ Dπ .
Theorem C.1. AssumeWo is sufficiently rich to include a “good” ω∗ inWo with D̂ω ∗ n ≈ Dπ in that
sup q∈Q ∣∣∣ED̂ω∗n [R̂q(x;x′, r)]− EDπ [R̂q(x;x′, r)]∣∣∣ ≤ cnα , (19) where c and α are two positive coefficients. Then we have
max { Ĵ+Q,W − Jπ, Jπ − Ĵ−Q,W } ≤ c nα + εn ‖ω∗‖Wo .
Assumption (19) holds if D̂n is collected following a Markov chain with certain strong mixing condition and weakly converges to some limit discussion D̂∞ whose support is X , for which we can define ω∗(x) = Dπ(x)/D∞(x). In this case, if Q is a finite ball in RKHS, then we can achieve (19) 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 function set R̂Q = {Rq(x, y) : q ∈ Q} Rakhlin et al. (2015). In our empirical reults, we observe that the gap of the practically constructed bounds tend to follow the O(n−1/2) rate.
Proof. Note that Jπ = EDπ [r] = EDπ [r],
and IQ(ω; D̂n) = sup
q∈Q
{ ED̂ωn [γq(x ′)− q(x)]− EDπ [γq(x′)− q(x)] } .
Because ω∗ ∈ W , we have Ĵ+W,Q − Jπ ≤ F̂+Q (ω∗)− Jπ
= ED̂ωn [r]− EDπ [r] + IQ(ωπ; D̂n) + εn ‖ω ∗‖Wo
= sup q∈Q
{ ED̂ωn [ R̂q(x, y) ] − EDπ [ R̂q(x, y) ]} + εn ‖ω∗‖Wo
≤ c nα + εn ‖ω∗‖Wo .
The case of lower bound follows similarly.
D OPTIMIZATION ONWo Consider the optimization of ω inWo
F̂+Q (ω) := 1
n n∑ i=1 riω(xi) + IQ(ω; D̂n) + ‖ω‖Wo
√ 2cqπ,k log(2/δ)
n (20)
AssumeWo is the RKHS of kernel k(x, x̄), that is,Wo = Hk. By the finite representer theorem of RKHS (Smola et al., 2007). the optimization of ω in RKHSHk can be reduced to a finite dimensional optimization problem. Specifically, the optimal solution of (20) can be written into a form of ω(x) = ∑n i=1 k(x, xi)αi with ‖ω‖ 2 Hk = ∑n i,j=1 k(xi, xj)αiαj for some vector α := [αi] n i=1 ∈ Rn. WriteK = [k(xi, xj)]ni,j=1 and r = [ri] n i=1. The optimization of ω reduces to a finite dimensional optimization on α:
min α∈Rn
1 n r>Kα+ IQ(Kα; D̂n) +
√ αKα
√ 2cqπ,k log(2/δ)
n ,
where
IQ(Kα; D̂n) = max q∈Q
{ EDπ,0 [q] + 1
n (T̂q)>Kα
} ,
and T̂q = [γq(x′i) − q(xi)]ni=1. When Q is RKHS, we can calculate IQ(Kα; D̂n) using (22) in section F.
This computation can be still expensive when n is large. Fortunately, our confidence bound holds for any ω; better ω only gives tighter bounds, but it is not necessary to find the 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(xi, θ), where g is any parametric function (such as a neural network) with parameter θ which can be much lower dimensional than α. We can then optimize θ with stochastic gradient descent, by approximating all the data averaging 1n ∑n i=1(·) 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.
2) Replacing kernel k Assume the kernel k yields a random feature expansion: k(x, x̄) = Eβ∼π[φ(x, β)φ(x̄, β)], where φ(x, β) is a feature map with parameter β and π is a distribution of β. We draw {βi}mi=1 i.i.d. from π, where m is taken to be much smaller than n. We replace k with k̂(x, x̄) = 1m ∑m i=1 φ(x, βi)φ(x̄, βi) andHk withHk̂, That is, we consider to solve
Ĵ+Q,W = min ω∈Hk̂ F̂+Q (ω) := 1n n∑ i=1 riω(xi) + IQ(ω; D̂n) + ‖ω‖Hk̂ √ 2cqπ,k̂ log(2/δ) n  . It is known that any function ω in Hk̂ can be represented as ω(x) = 1m ∑m i=1 wiφ(x, βi), for some w = [wi]mi=1 ∈ Rm and satisfies ‖ω‖2Hk̂ = 1 m ∑m i=1 w 2 i . In this way, the problem reduces to optimizing an m-dimensional vector w, which can be solved by standard convex optimization techniques.

heading: E CONCENTRATION INEQUALITIES OF GENERAL FUNCTIONAL BELLMAN LOSSES
text: When K is a general function set, one can still obtain a general concentration bound using Radermacher complexity. Define R̂q ◦ W := {h(x, y) = R̂q(x, y)ω(x) : ω ∈ W}. Using the standard derivation in Radermacher complexity theory in conjunction with Martingale theory (Rakhlin et al., 2015), we have
sup ω∈W
{ 1
n n∑ i=1
(R̂q(xi, yi)−Rπq(xi))ω(xi) } ≤ 2Rad(R̂q ◦W) + √ 2cq log(2/δ)
n ,
where Rad(R̂q ◦ K) is the sequential Radermacher complexity as defined in (Rakhlin et al., 2015). A triangle inequality yields
| Lk(q; D̂n)− Lk(q; D̂n) | ≤ sup ω∈W
{ 1
n n∑ i=1
(R̂q(xi, yi)−Rπq(xi))ω(xi) }
Therefore,
| LW(q; D̂n)− LW(q; D̂n) | ≤ 2Rad(R̂q ◦W) + √ 2cq log(2/δ)
n , (21)
where cq,W = supω∈W supx,y(R̂q(x, y)−Rπq(x))2ω(x)2. WhenW equals the unit ball K of the RKHS related to kernel k, we have cq,k = cq,W , and hence this bound is strictly worse than the bound in Theorem 4.1.
F CLOSED FORM OF IQ(ω; D̂n) WHEN Q IS RKHS
Similar to LK(q; D̂n), whenQ is taken to be the unit ball K̃ of the RKHS of a positive definite kernel k̃(x, x̄), (7) can be expressed into a bilinear closed form shown in Mousavi et al. (2020):
IQ(ω; D̂n) 2 = A− 2B + C, (22)
where A = E(x,x̄)∼Dπ,0×Dπ,0 [k(x, x̄)]
B = E(x,x̄)∼D̂ωn×Dπ,0 [ T̂xπk(x, x̄) ] C = E(x,x̄)∼D̂ωn×D̂ωn [ T̂xπT̂ x̄ πk(x, x̄) ] ,
were T̂πf(x) = γf(x′) − f(x); in T̂xπT̂x̄πk(x, x̄), we apply T̂x̄π and T̂xπ in a sequential order by treating k as a function of x̄ and then of x.

heading: G MORE ON THE ORACLE BOUND AND ITS DUAL FORM
text: The oracle bound (18) provides another starting point for deriving optimization-based confidence bounds. We derive its due form here. Using Lagrangian multiplier, the optimization in (18) can be rewritten into
Ĵ+Q,∗ = max q∈Q min ω M(q, ω; D̂n), (23)
where M∗(q, ω; D̂n) = EDπ,0 [q]− 1
n n∑ i=1 ω(xi) ( R̂q(xi, yi)− R̂qπ(xi, yi) ) ,
where ω now serves as the Lagrangian multiplier. By the weak duality, we have J∗Q,+ ≤ F̂+Q,∗(ω) := ED̂ωn [r] + IQ(ω; D̂n)︸ ︷︷ ︸
known
+R(ω, qπ)︸ ︷︷ ︸ unknown , ∀ω.
and
R(ω, qπ) = 1
n n∑ i=1 ω(xi)R̂qπ(xi).
The derivation follows similarly for the lower bound. So for any ω ∈ Wo, we have [Ĵ−Q,∗, Ĵ+Q,∗] ⊆ [F̂−Q,∗(ω), F̂ + Q,∗(ω)].
Here the first two terms of F̂+Q,∗(ω) can be empirically estimated (it is the same as the first two terms of (16)), but the third term R(ω, qπ) depends on the unknown qπ and hence need to be further upper bounded.
Our method can be viewed as constraining ω inW , which is assumed to be the unit ball ofWo, and applying a worst case bound:
F̂+Q,∗(ω) := ED̂ωn [r] + IQ(ω; D̂n) +R(ω, qπ), ∀ω ∈ Wo ≤ ED̂ωn [r] + IQ(ω; D̂n) + ‖w‖Wo suph∈W R(h, qπ), ∀ω ∈ Wo
≤ ED̂ωn [r] + IQ(ω; D̂n) + ‖w‖Wo LW(qπ, D̂n), ∀ω ∈ Wo w.p.1−δ ≤ ED̂ωn [r] + IQ(ω; D̂n) + ‖w‖Wo , ∀ω ∈ Wo
= F̂+Q (ω).
where the last step applies the high probability bound that Pr(LW(qπ, D̂n) ≤ ε) ≥ 1− δ. Similar derivation on the lower bound counterpart gives
Pr ([ F̂−Q,∗(ω), F̂ + Q,∗(ω) ] ⊆ [ F̂−Q (ω), F̂ + Q (ω) ]) ≥ 1− δ.
Therefore, our confidence bound [F̂−Q (ω), F̂ + Q (ω)] is a 1− δ confidence outer bound of the oracle bound [Ĵ−Q,∗, Ĵ + Q,∗] ⊆ [F̂−Q,∗(ω), F̂+Q,∗(ω)].
Introducing Q is necessarily Our method does not require any independence assumption between the transition pairs, the trade-off is that that we have to assume that qπ falls into a function set Q that imposes certain smoothness assumption. This is necessary because the data only provide information regarding qπ on a finite number of points, and qπ can be arbitrarily non-smooth outside of the data points, and hence no reasonable upper/lower bound can be obtained without any smoothness condition that extend the information on the data points to other points in the domain.
Proposition G.1. Unless Prs∼Dπ,0(s /∈ {si, s′i}ni=1) = 0, for any u ∈ R, there exists a function q : S ×A → R, such that
EDπ,0 [q] = u, R̂q(xi, yi) = R̂qπ(xi, yi), ∀i = 1, . . . , n.
Proof. Let Qnull be the set of functions that are zero on {si, s′i}ni=1, that is, Qnull = {g : S ×A → R : g(s, a) = 0, ∀s ∈ {si, s′i}ni=1, a ∈ A}.
Then we have R̂π(qπ + g)(xi, yi) = R̂πqπ(xi, yi), ∀i = 1, . . . , n.
and EDπ,0 [qπ + g] = EDπ,0[qπ] + EDπ,0 [g] = Jπ + EDπ,0 [g]. Taking g(s, a) = zI(s /∈ {si, s′i}ni=1), where z is any real number. Then we have EDπ,0 [qπ + g] = Jπ + zPrs∼Dπ,0(s /∈ {si, s′i}ni=1).
Because Prs∼Dπ,0(s /∈ {si, s′i}ni=1). 6= 0, we can take z to be arbitrary value to make EDπ,0 [qπ + g] to take arbitrary value.

heading: H ABLATION STUDY AND EXPERIMENTAL DETAILS
text: H.1 EXPERIMENTAL DETAILS
Environments and Dataset Construction We test our method on three environments: InvertedPendulum 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.(16)-(17). Throughout the experiment, we always set W = K, a unit ball of RKHS with kernel k(·, ·). We set Q = rQK̃, the zero-centered ball of radius rQ in an RKHS with kernel k̃(·, ·). We take both k and k̃ to be Gaussian RBF kernel. The bandwidth of k and k̃ 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 a q̂ so that its functional Bellman loss is small than n. Then we set rQ = 10 ∗ ‖q̂‖K̃. We optimize ω using the random feature approximation method described in Appendix D. Once ω+ and ω− are found, we evaluate the bound in Eq. (16) 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 andW .
Radius of Q Recall that we choose Q to be a ball in RKHS with radius rQ, that is,
Q = rQK̃ = {rQf : f ∈ K̃},
where K̃ is the unit ball of the RKHS with kernel k̃. Ideally, we want to ensure that rQ ≥ ‖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∗ inside K̃ and then infer the corresponding deterministic reward function r(x) by inverting the Bellman equation:
r(x) := q∗(x)− γEx′∼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∗:
ρ∗ = ‖q∗‖K̃ . 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 rQ 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× 103 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 and k̃ forW 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 of k̃ 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 dataset D̂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.
2. The dataset D̂n is the combination of the data collected from multiple behavior policies of form πα defined as above, with α ∈ {0.0, 0.2, 0.4, 0.6, 0.8}.
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 in D̂n As we collect D̂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 collect D̂n, while fixing the total number n of transition pairs. In this way, the number of trajectories in each D̂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).

heading: I MORE RELATED WORKS
text: We give a more detailed 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 et al., 2000; 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 := {sit, ait, rit}T−1t=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:
Ĵ IS = 1
m m∑ i=1 ω(τi)J(τi) , where ω(τi) = T−1∏ t=0 π(ait|sit) π0(ait|sit) , J(τi) = T−1∑ t=0 γtrit . (24)
One can construct non-asymptotic confidence bounds based on Ĵ 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 & Li, 2016; Thomas & 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 nonasymptotic 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.
Infinite-Horizon, Behavior-Agnostic OPE Our work is closely related to the recent advances in infinite-horizon and behavior-agnostic OPE, including, for example, Liu et al. (2018a); Feng et al. (2019); Tang et al. (2020a); Mousavi et al. (2020); Liu et al. (2020); Yang et al. (2020b); Xie et al. (2019); Yin & Wang (2020), as well as the DICE-family (e.g., Nachum et al., 2019a;b; Zhang et al., 2020a; Wen et al., 2020; Zhang et al., 2020b). 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; Jiang & Huang, 2020) that we elaborate in depth in Section 3.
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; Jiang & Huang, 2020; Duan et al., 2020; Dai et al., 2020; Feng et al., 2020; Tang et al., 2020b; Yin et al., 2020; Lazic et al., 2020). For example, Dai et al. (2020) 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 & 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 estimator P̂ of P using maximum likelihood estimation or other methods, and plug it into Eq. (1) to obtain a plug-in estimator Ĵ = Jπ,P̂. This yields the model-based approach to OPE (e.g., Jiang & Li, 2016; Liu et al., 2018b). One can also estimate the uncertainty in Jπ,P̂ by propagating the
uncertatinty in P̂ (e.g., Asadi et al., 2018; Duan et al., 2020), but it is hard to obtain non-asymptotic
and computationally efficient bounds unless P̂ 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 & White, 2010; Hanna et al., 2017; Kostrikov & 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 & 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., 2016b; 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.

Review summaries:
The paper studies an off-policy evaluation (OPE) problem for Markov decision processes (MDPs). It suggests an optimization-based method that can construct a non-asymptotic confidence interval, for a given confidence level, for the value function of a policy starting from a fixed initial distribution. The paper builds on the works of Feng et al. (2019, 2020); the main advantages of the current work with respect to the previous methods are that the suggested approach guarantees a faster convergence rate, it does not require full independence between transition pairs, and it does not need the global optimal solution of the underlying optimization problem, in order to construct guaranteed confidence intervals. The authors present some theoretical results about the construction, including a discussion on the special case of using RKHS approaches, and also present numerical experiments on benchmark problems, such as the inverted-pendulum, cartpole and type-1 diabetes.

This work constructs non-asymptotic confidence intervals for off-policy evaluation. This is achieved by assuming that the reward at any given time only depends on the state action pair, leveraging that assumed structure to define the difference between the empirical and estimated bellman residual operators as a Martingale difference sequence. This, in turn, then allows the authors to apply a Hoeffding-like concentration inequality which applies to Hilbert spaces. The authors then provide a derivation of the confidence bounds by considering the divergence between policies. The work improves on the rate of prior work from $O(n^{-\frac{1}{4}})$ to $O(n^{-\frac{1}{2}})$ and allows for estimation without the need of global optimality via the dual formulation, both of which are very nice additions to the literature. Experimental evaluation backs up the authors’ claims, showing very strong performance with respect to prior art.