Title: Efficient Inference and Exploration for Reinforcement Learning

Abstract: Despite an ever growing literature on reinforcement learning algorithms and applications, much less is known about their statistical inference. In this paper, we investigate the large-sample behaviors of the Q-value estimates with closed-form characterizations of the asymptotic variances. This allows us to efficiently construct confidence regions for Q-value and optimal value functions, and to develop policies to minimize their estimation errors. This also leads to a policy exploration strategy that relies on estimating the relative discrepancies among the Q estimates. Numerical experiments show superior performances of our exploration strategy than other benchmark approaches.

1. Introduction:

We consider the classical reinforcement learning (RL) problem where the agent interacts with a random environment and aims to maximize the accumulated discounted reward over time. The environment is formulated as a Markov decision process (MDP) and the agent is uncertain about the true dynamics to start with. As the agent interacts with the environment, data about the system dynamics are collected and the agent becomes increasingly confident about her decision. With finite data, however, the potential reward from each decision is estimated with errors and the agent may be led to a suboptimal decision. Our focus in this paper is on statistically efficient methodologies to quantify these errors and uncertainties, and to demonstrate their use in obtaining better policies.

More precisely, we investigate the large-sample behaviors of estimated Q-value, optimal value function, and their associated policies. Our results are in the form of asymptotic convergence to an explicitly identified and computable Gaussian (or other) distribution, as the collected data sizes increase. The motivation of our investigation is three-fold. First, these precise asymptotic statements allow us to construct accurate confidence regions for quantities related to the optimal policy, and, like classical statistical inference, they can assess the reliability of the current estimates with respect to the data noises. Second, our results complement some finite-sample error bounds developed in the literature (Kearns & Singh, 1998; Kakade, 2003; Munos & Szepesvári, 2008) , by supplementing a closed-form asymptotic variance that often shows up in the first-order terms in these bounds.

Our third and most important motivation is to design good exploration policies by directly using our tight error estimates. Motivated by recent autonomous-driving and other applications (e.g., Kalashnikov et al. (2018) ), we consider the pure exploration setting where an agent is first assigned an initial period to collect as much experience as possible, and then, with the optimal policy trained offline, starts deployment to gain reward. We propose an efficient strategy to explore by optimizing the worst-case estimated relative discrepancy among the Q-values (ratio of mean squared difference to variance), which provides a proxy for the probability of selecting the best policy. Similar criteria have appeared in the so-called optimal computing budget allocation (OCBA) procedure in simulation-based optimization (Chen & Lee, 2011) (a problem closely related to best-arm identification (Audibert & Bubeck, 2010) in online learning). In this approach, one divides computation (or observation) budget into stages in which one sequentially updates mean and variance estimates, and optimizes next-stage budget allocations according to the worst-case relative discrepancy criterion. Our proposed procedure, which we term Q-OCBA, follows this idea with a crucial use of our Q-value estimates and randomized policies to achieve the optimal allocation. We demonstrate how this idea consistently outperforms other benchmark exploration policies, both in terms of the probability in selecting the best policy and generating the tightest confidence bounds for value estimates at the end of the exploration period.

Regarding the problem of constructing tight error estimates in RL, the closest work to ours is Mannor et al. (2004; 2007) , which studies the bias and variance in value function estimates with a fixed policy. Our technique resolves a main technical challenge in Mannor et al. (2004; 2007) , which allows us to substantially generalize their variance results to Q-values, optimal value functions and asymptotic distributional statements. The derivation in Mannor et al. (2004; 2007) hinges on an expansion of the value function in terms of the perturbation of the transition matrix, which (as pointed out by the authors) is not easily extendable from a fixed-policy to the optimal value function. In contrast, our results utilize an implicit function theorem applied to the Bellman equation that can be verified to be sufficiently smooth. This idea turns out to allow us to obtain gradients for Q-values, translate to the optimal value function, and furthermore generalize to similar results for constrained MDP and approximate value iterations. We also relate our work to the line of studies on dynamic treatment regimes (DTR) (Laber et al., 2014) applied commonly in medical decision-making, which focuses on the statistical properties of polices on finite horizon (such as two-period). Our infinite-horizon results on the optimal value and Q-value distinguishes our developments from the DTR literature. Moreover, our result on the non-unique policy case can be demonstrated to correspond to the "non-regularity" concept in DTR, where the true parameters are very close to the decision "boundaries" that switch the optimal policy (motivated by situations of small treatment effects), thus making the obtained policy highly sensitive to estimation noises.

In the rest of this paper, we first describe our MDP setup and notations (Section 2). Then we present our results on large-sample behaviors (Section 3), demonstrate their use in exploration strategies (Section 4), and finally substantiate our findings with experimental results (Section 5). In the Appendix, we first present generalizations of our theoretical results to constrained MDP (A.1) and problems using approximate value iteration (A.2). Then we include more numerical experiments (B), followed by all the proofs (C).

2. Problem Setup:

Consider an infinite horizon discounted reward MDP, M = (S, A, R, P, γ, ρ), where S is the state space, A is the action space, R(s, a) denotes the random reward when the agent is in state s ∈ S and selects action a ∈ A, P (s |s, a) is the probability of transitioning to state s in the next epoch given current state s and taken action a, γ is the discount factor, and ρ is the initial state distribution. The distribution of the reward R and the transition probability P are unknown to the agent. We assume both S and A are finite sets. Without loss of generality, we denote S = {1, 2, . . . , m s } and A = {1, 2, . . . , m a }. Finally, we make the following stochasticity assumption: Assumption 1. R(s, a) has finite mean µ R (s, a) and finite variance σ 2 R (s, a) ∀ s ∈ S, a ∈ A. For any given s ∈ S and a ∈ A, R(s, a) and S ∼ P (•|s, a) are all independent random variables.

A policy π is a mapping from each state s ∈ S to a probability measure over actions a ∈ A. Specifically, we write π(a|s) as the probability of taking action a when the agent is in state s and π(•|s) as the m adimensional vector of action probabilities at state s. For convenience, we sometimes write π(s) as the realized action given the current state is s. The value function associated with a policy π is defined as

V π (s) = E[ ∞ t=0 γ t R(s t , π(s t ))|s 0 = s] with s t+1 ∼ P (.|s t , π(s t )).

The expected value function, under the initial distribution ρ, is denoted by χ π = s ρ(s)V π (s). A policy π * is said to be optimal if V π * (s) = max π V π (s) for all s ∈ S. For convenience, we denote

EQUATION

for any (s, a) ∈ S × A. Denoting the Bellman operator as

T µ R ,P (•), Q is a fixed point associated with T µ R ,P , i.e. Q = T µ R ,P (Q).

For the most part of this paper we make the following assumption about Q: Assumption 2. For any state s ∈ S, arg max a∈A Q(s, a) is unique.

Under Assumption 2, the optimal policy π * is unique and deterministic. Let a * (s) = arg max a∈A Q(s, a). Then π * (a|s) = 1 (a = a * (s)), where 1(•) denotes the indicator function.

We next introduce some statistical quantities arising from data. Suppose we have n observations (whose collection mechanism will be made precise later), which we denote as {(s t , a t , r t (s t , a t ), s t (s t , a t )) : 1 ≤ t ≤ n}, where r t (s t , a t ) is the realized reward at time t and s t (s t , a t ) = s t+1 . We define the sample mean μR,n and the sample variance σ2 R,n of the reward as

EQUATION

Similarly, we define the empirical transition matrix Pn as

EQUATION

and its m s × m s sampling covariance matrix Σ Ps,a (with one sample point of 1(s t = s, a t = a)) as

Σ Ps,a (k 1 , k 2 ) = P (k 1 |s, a)(1 -P (k 1 |s, a)) k 1 = k 2 -P (k 1 |s, a)P (k 2 |s, a) k 1 = k 2 . , for 1 ≤ k 1 ≤ m s , 1 ≤ k 2 ≤ m s .

With the data, we construct our estimate of Q, called Qn , which is the empirical fixed point of T μR,n , Pn , i.e. Qn = T μR,n , Pn ( Qn ). Correspondingly, we also write V * n (s) = max a∈A Qn (s, a) and χ * n = s∈S ρ(s) V * n (s). We shall focus on the empirical errors due to noises of the collected data, and assume the MDP or Q-value evaluation can be done off-line so that the fixed point equation for Qn can be solved exactly. .

3. Quantifying Asymptotic Estimation Errors:

We present an array of results regarding the asymptotic behaviors of Qn and V * n . To prepare, we first make an assumption on our exploration policy π to gather data. Define the extended transition probability P π as P π (s , a |s, a) = P (s |s, a)π(a |s ). We make the assumption: Assumption 3. The Markov chain with transition probability P π is positive recurrent.

Under Assumption 3, P π has a unique stationary distribution, denoted w, equal to the long run frequency in visiting each state-action pair, i.e. w(s, a) = lim n→∞ 1 n 1≤t≤n 1(s t = i, a t = j), where all w(s, a)'s are positive. Note that Assumption 3 is satisfied if for any two states s, s , there exists a sequence of actions such that s is attainable from s under P , and, moreover, if π is sufficiently mixed, e.g., π satisfies π(a |s ) > 0 for all s , a .

Our results in the sequel uses the following further notations. We denote "⇒" as "convergence in distribution", and N (µ, Σ) as a multivariate Gaussian distribution with mean vector µ and covariance matrix Σ. We write I as the identity matrix, and e i as the i-th unit vector. The dimension of N (µ, Σ), I and e i should be clear from the context. When not specified, all the vectors are column vectors. Let N = m s m a . In our algebraic derivations, we need to re-arrange µ R , Q and w as N -dimensional vectors. We thus define the following indexing rule: (s = i, a = j) is re-indexed as (i -1)m a + j, e.g. µ R (i, j) = µ R ((i -1)m a + j). We also need to re-arrange P π as an N × N matrix following the same indexing rule, i.e. P π (i , j |i, j) = P π ((i -1)m a + j, (i -1)m a + j ).

3.1. Limit Theorems under Sufficient Exploration:

We first establish the asymptotic normality of Qn under exploration policy π: Theorem 1. Under Assumptions 1 and 2, if the data is collected according to π satisfying Assumption 3, then Qn is a strongly consistent estimator of

Q, i.e. Qn → Q almost surely as n → ∞. Moreover, √ n( Qn -Q) ⇒ N (0, Σ) as n → ∞, where Σ = (I -γ P π * ) -1 W -1 (D R + D Q )((I -γ P π * ) -1 ) T , (5) W , D R and D Q are N × N diagonal matrices with W ((i -1)m a + j, (i -1)m a + j) = w(i, j), D R ((i -1)m a + j, (i -1)m a + j) = σ 2 R (i, j) and D Q ((i -1)m a + j, (i -1)m a + j) = (V * ) T Σ Pi,j V * respectively .

In addition to the asymptotic Gaussian behavior, a key element of Theorem 1 is the explicit form of the asymptotic variance Σ. This is derived from the delta method (Serfling, 2009) and, intuitively, is the product of the sensitivities (i.e., gradient) of Q with respect to its parameters and the variances of the parameter estimates. Here the parameters are µ R and P , with corresponding gradients (I -γ P π * ) -1 and (I -γ P π * ) -1 V * . The variances of these parameter estimates (i.e., (2) and ( 4)) involve σ 2 R (i, j) and Σ Pi,j , and the sample size allocated to estimate each parameter, which is proportional to w(i, j).

Using the relations that V * n (s) = max a∈A Q(s, a) and V * n (s) = max a∈A Qn (s, a), we can leverage Theorem 1 to further establish the asymptotic normality of V * n and χ * n : Corollary 1. Under Assumptions 1, 2 and 3,

√ n( V * n -V * ) ⇒ N (0, Σ V ) and √ n( χ * n -χ * ) ⇒ N (0, σ 2 χ ) as n → ∞ where Σ V = (I -γP π * ) -1 (W π * ) -1 [D π * R + D π * V ]((I -γP π * ) -1 ) T , σ 2 χ = ρ T Σ V ρ, P π * is an m s × m s transition matrix with P π * (i, j) = P (j|s = i, a = a * (s)), W π * , D π * R and D π * V are m s × m s diagonal matrices with W π * (i, i) = w(i, a * (i)), D π * R (i, i) = σ 2 R (i, a * (i)) and D π * V (i, i) = (V * ) T Σ P i,a * (i) V * respectively.

In the Appendix we also prove, using the same technique as above, a result on the large-sample behavior of the value function for a fixed policy (Corollary 2), which essentially recovers Corollary 4.1 in Mannor et al. (2007) . Different from Mannor et al. (2007) , we derive our results by using an implicit function theorem on the corresponding Bellman equation to obtain the gradient of Q, viewing the latter as the solution to the equation and as a function of µ R , P . This approach is able to generalize the results for fixed policies in Mannor et al. (2007) to the optimal value functions, and also provide distributional statements as Theorem 1 and Corollary 1 above. We also note that another potential route to obtain our results is to conduct perturbation analysis on the linear program (LP) representation of the MDP, which would also give gradient information of V * (and hence also Q), but using the implicit function theorem here seems sufficient.

Theorem 1 and Corollary 1 can be used immediately for statistical inference. In particular, we can construct confidence regions for subsets of the Q-value jointly, or for linear combinations of the Q-values. A quantity of interest that we will later utilize in designing good exploration policies is Q(s, a 1 ) -Q(s, a 2 ), i.e. the difference between action a 1 and a 2 when the agent is in state s. Define σ 2 ∆Q as

EQUATION

and its estimator σ2 ∆Q,n by replacing

Q, V * , σ 2 R,n , w, P with Qn , V * n σ2 R,n , ŵn , Pn in Σ,

where ŵn is the empirical frequency of visiting each state-action pair, i.e. ŵn (i, j)

= 1 n 1≤t≤n 1(s t = i, a t = j). Then the 100(1 -α)% confidence interval (CI) for Q(s, a 1 ) -Q(s, a 2 ) takes the form Qn (s, a 1 ) -Qn (s, a 2 ) ± z α σ2 ∆Q,n (s, a 1 , a 2 )

, where z α is the (1 -α/2)-quantile of N (0, 1).

3.2. Non-Unique Optimal Policy:

Suppose the optimal policy for the MDP M is not unique, i.e., Assumption 2 does not hold. In this situation, the estimated Qn and V * n may "jump" around different optimal actions, leading to a more complicated large-sample behavior as described below:

Theorem 2. Suppose Assumptions 1 and 3 hold but there is no unique optimal policy. Then there exists

K ≥ 1 distinct m s ×(N m s +N ) matrices {G k } 1≤k≤K and a deterministic partition of U = {u ∈ R msN +ms : ||u|| = 1} = ∪ 1≤k≤K U k such that √ n( V * n -V * ) ⇒ K k=1 G k 1 (Z/||Z|| ∈ U k ) Z, where Z = N (0, Σ R,P ), Σ R,P = Diag(W -1 D R , D P ) and D P = Diag(Σ P1,1 /w(0m a + 1), . . . , Σ Pi,j /w((i - 1)m a + j), . . . , Σ Pm s,ma /w((m s -1)m a + m a )).

In the case that K > 1 in Theorem 2, the limit distribution becomes non-Gaussian. This arises because the sensitivity to P or µ R can be very different depending on the perturbation direction, which is a consequence of solution non-uniqueness that can be formalized as a non-degeneracy in the LP representation of the MDP. We note that this phenomenon is analogous to the "non-regularity" concept in DTR that arises because the "true" parameters in these problems are very close to the decision "boundaries", which makes the obtained policy highly sensitive to estimation noises and incurs a 1/ √ n-order bias behavior. Our case of non-unique optimal policy here captures precisely this same behavior, where we see in Theorem 2 that when K > 1 the asymptotic limit no longer has mean zero and consequently a 1/ √ n-order bias arises.

We also develop two other generalizations of large-sample results, for constrained MDP and approximate value iteration respectively (see Appendices A.1 and A.2).

4. Efficient Exploration Policy:

We utilize our results in Section 3 to design exploration policies. We focus on the setting where an agent is assigned a period to collect data by running the state transition with an exploration policy.

The goal is to obtain the best policy at the end of the period in a probabilistic sense, i.e., minimize the probability of selecting a suboptimal policy for the accumulated reward.

We propose a strategy that maximizes the worst-case relative discrepancy among all Q-value estimates. More precisely, we define, for i ∈ S, j ∈ A and j = a * (i), the relative discrepancy as

h ij = (Q(i, a * (i)) -Q(i, j)) 2 /σ 2 ∆Q (i, a * (i), j), where σ 2 ∆Q (i, a * (i), j) is defined in (6).

Our procedure attempts to maximize the minimum of h ij 's, max

EQUATION

where w denotes the proportions of visits on the state-action pairs, within some allocation set W η (which we will explain). Intuitively, h ij captures the relative "difficulty" in obtaining the optimal policy given the estimation errors of Q's. If the Q-values are far apart, or if the estimation variance is small, then h ij is large which signifies an "easy" problem, and vice versa. Criterion ( 7) thus aims to make the problem the "easiest". Alternatively, one can also interpret ( 7) from a large deviations view (Glynn & Juneja, 2004; Dong & Zhu, 2016) . Suppose the Q-values for state i between two different actions a * (i) and j are very close. Then, one can show that the probability of suboptimal selection between the two has roughly an exponential decay rate controlled by h ij . Obviously, there can be many more comparisons to consider, but the exponential form dictates that the smallest decay rate dominates the calculation, thus leading to the inner min's in (7). Criterion like ( 7) is motivated from the OCBA procedure in simulation optimization (which historically has considered simple mean-value alternatives (Chen & Lee, 2011) ). Here, we consider the Q-values. For convenience, we call our procedure Q-OCBA.

Implementing criterion (7) requires two additional considerations. First, solving (7) needs the model primitives Q, P and σ 2 R that appear in the expression of h ij . These quantities are unknown a priori, but as we collect data they can be sequentially estimated. This leads to a multi-stage optimization plus parameter update scheme. Second, since data are collected through running a Markov chain on the exploration actions, not all allocation w is admissible, i.e., realizable as the stationary distribution of the MDP. To resolve this latter issue, we will derive a convenient characterization for admissibility.

Call π(•|s) admissible if the Markov Chain with transition probability P π , defined for Assumption 3, is positive recurrent, and denote w π as its stationary distribution. Define the set W = w > 0 : 1≤j≤ma w((i -1)m a + j) = 1≤k≤ms 1≤l≤ma w((k -1)m a + l)P (i|k, l)

∀1 ≤ i ≤ m s , 1≤i≤ms 1≤j≤ma w((i -1)m a + j) = 1 .

The following provides a characterization of the set of admissible π: Lemma 1. For any admission policy π, w π ∈ W. For any w ∈ W, π w with π w (a = j|s

= i) = w((i -1)m a + j)/ ( ma k=1 w((i -1)m a + k)) is an admissible policy.

In other words, optimizing over the set of admissible policies is equivalent to optimizing over the set of stationary distributions. The latter is much more tractable thanks to the linear structure of W.

In practice, we will use W η = W ∩ {w ≥ η} for some small η > 0 to ensure closedness of the set (our experiments use η = 10 -6 ).

Algorithm 1 describes Q-OCBA. In our experiments shown next, we simply use two stages, i.e., K = 2. Finally, we also note that criterion like (7) can be modified according to the decision goal.

For example, if one is interested in obtaining the best estimate of χ * , then it would be more beneficial to consider min w∈Wη σ 2 χ . We showcase this with additional experiments in the Appendix.

Input: Number of iterations K, length of each batch {B k } 1≤k≤K , initial exploration policy π 0 ; 

Initialization: k = 0; while k ≤ K do Run π k for B k steps and set k = k + 1; Calculate PB k , μR,B k ,σ 2 R,

; Set π k (a = j|s = i) = w k ((i -1)m a + j)/ ma l=1 w k ((i -1)m a + l); end

Algorithm 1: Q-OCBA sequential updating rule for exploration Note that ( 7) is equivalent to min w max i∈S max j∈A,j =a * (i) s,a c ij (s, a)/w s,a subject to w ∈ W η , where c ij (s, a)'s are non-negative coefficients. Based on the closed-form characterization of Σ in Theorem 1, c ij (s, a)'s can be estimated with plug-in estimators using data collected in earlier stages.

5. Numerical Experiments:

We conduct several numerical experiments to support our large-sample results in Sections 3 and demonstrate the performance of Q-OCBA against some benchmark methods. We use the RiverSwim problem in (Osband et al., 2013) with m s states and two actions at each state: swim left (0) or swim right (1) (see Figure 1 ). The triplet above each arc represents i) the action, 0 or 1, ii) the transition probability to the next state given the current state and action, iii) the reward under the current state and action. Note that, in this problem, rewards are given only at the left and right boundary states (where the value of r L will be varied). We consider the infinite horizon setting with γ = 0.95 and ρ = [1/m s , . . . , 1/m s ] T .

Figure 1 : RiverSwim Problem

We first demonstrate the validity of our large-sample results. We use a policy that swims right with probability 0.8 at each state, i.e. π(1|s) = 0.8. Tables 1 and 2 show the coverage rates of the constructed 95% CIs, for a small m s = 6 (using Theorem 1 and Corollary 1) and a large m s = 31 (using Theorem 4 in the Appendix) respectively. The latter case uses a linear interpolation with S 0 = {1, 4, . . . , 28, 31}. All coverage rates are estimated using 10 3 independent experimental repetitions (the bracketed numbers in the tables show the half-widths of 95% CI for the coverage estimates). For the Q-values, we report the average coverage rate over all (s, a) pairs. When the number of observations n is large enough (≥ 3 × 10 4 for exact update and ≥ 10 5 for interpolation), we see highly accurate CI coverages, i.e., close to 95%. Next we investigate the efficiency of our exploration policy. We compare Q-OCBA with K = 2 to four benchmark policies: i) -greedy with different values of , ii) random exploration (RE) with different values of π(1|s), iii) UCRL2 (a variant of UCRL) with δ = 0.05 (Jaksch et al., 2010) , iv) PSRL with different posterior updating frequencies (Osband et al., 2013) , i.e., PSRL(x) means PSRL is implemented with x episodes. We use m s = 6 and vary r L from 1 to 3. To ensure fairness, we use a two-stage implementation for all policies, with 30% of iterations first dedicated to RE (with π(1|s) = 0.6) as a warm start, i.e., the data are used to estimate the parameters needed for the second stage. To give enough benefit of the doubt, we notice the probabilities of correct selection for both UCRL2 and PSRL are much worse without the warm start.

Tables 3 and 4 compare the probabilities of obtaining the optimal policy (based on the estimated Qn 's). For -greedy, RE, and PSRL, we report the results with the parameters that give the best performances in our numerical experiments. The probability of correct selection is estimated using 10 3 replications of the procedure. We observe that Q-OCBA substantially outperforms the other methods, both with a small data size (n = 10 3 in Table 3 ) and a larger one (n = 10 4 in Table 4 ). Generally, these benchmark policies perform worse for larger values of r L . This is because for small r L , the (s, a) pairs that need to be explored more also tend to have larger Q-values. However, as r L increase, there is a misalignment between the Q-values and the (s, a) pairs that need more exploration.

The superiority of our Q-OCBA in these experiments come as no surprise to us. The benchmark methods like UCRL2 and PSRL are designed to minimize regret which involves balancing the exploration-exploitation trade-off. On the other hand, Q-OCBA focuses on efficient exploration only, i.e., our goal is to minimize the probability of incorrect policy selection, and this is achieved by carefully utilizing the variance information gathered from the first stage that is made possible by our derived asymptotic formulas. We provide additional numerical results in Appendix B. 

A Additional Theoretical Results:

In this section, we present additional results on large-sample behaviors for constrained MDPs and also estimations based on approximation value iteration.

A.1 Constrained Problems:

We consider the constrained MDP setting for budgeted decision-making (Boutilier & Lu, 2016 ) and more recently safety-critical applications (Achiam et al., 2017; Chow et al., 2017) . Suppose now we aim to maximize the long-run accumulated discounted reward,

V π (s) = E[ ∞ t=0 γ t R(s t , π(s t ))|s 0 = s],

while at the same time want to ensure that a long-run accumulated discounted cost, denoted as L π (s) = E[ ∞ t=0 γ t C(s t , π(s t ))|s 0 = s] which we call the loss function, is constrained by some given value η, i.e., max

EQUATION

We assume data coming in like before and, in addition, that we have observations on the incurred cost at each sample of (s, a). Call the empirical estimate of the cost μC,n . We follow our paradigm to solve the empirical counterpart of the problem, namely to find a policy π * n that solves (8) by using V π n (s) and Lπ n (s) instead of V π (s) and L π (s), where V π n (s)'s and Lπ n (s)'s are the value functions and loss functions evaluated using the empirical estimates μR,n , μC,n , Pn . We focus on the estimation error of the optimal value (instead of the feasibility, which could also be important but not pursued here).

To understand the error, we first utilize an optimality characterization of constrained MDPs. In general, an optimal policy for (8) is a "split" policy (Feinberg & Rothblum, 2012) , namely, a policy that is deterministic except that at one particular state a randomization between two different actions is allowed. This characterization can be deduced from the associated LP using occupancy measures (Altman, 1999) . We call the randomization probability the mixing parameter α * , i.e., whenever this particular state, say s r , is visited, action a * 1 (s r ) is chosen with probability α * and action a * 2 (s r ) is chosen with probability 1 -α * . We then have the following result: Theorem 3. Suppose Assumptions 1 and 3 hold and there is a unique optimal policy. Moreover, assume that there is no deterministic policy π that satisfies s ρ(s)L π (s) = η. Then we have √ n( V * n -V * ) ⇒ N (0, Σ) as n → ∞, where one of the following cases hold:

1. The optimal policy is deterministic. We then have Σ = Σ V where Σ V is defined in Theorem 1.

2. The optimal policy is deterministic, except at one state where a randomization between two actions occurs, with the mixing parameter α * . Denote the state where the randomization occurs by s r and two possible actions for s r by a * 1 (s r ) and a * 2 (s r ), We have

Σ = (I -γP π * ) -1 [G π * , 0, H π * V ] - (I -γP π * ) -1 q V ρ T (I -γP π * ) -1 [0, G π * , H π * L ] ρ T (I -γP π * ) -1 q L Σ R,C,P (I -γP π * ) -1 [G π * , 0, H π * V ] - (I -γP π * ) -1 q V ρ T (I -γP π * ) -1 [0, G π * , H π * L ] ρ T (I -γP π * ) -1 q L T , where Σ R,C,P = Diag(W -1 D R , W -1 D C , D P ), q V and q L are m s -dimensional vectors with q V (s) = (µ R (s, a * 1 (s))-µ R (s, a * 2 (s))+ ms j=1 γV π * (j)(P (j|s, a * 1 (s))-P (j|s, a * 2 (s))) and q L (s) = (µ C (s, a * 1 (s))- µ C (s, a * 2 (s)) + ms j=1 γL π * (j)(P (j|s, a * 1 (s)) -P (j|s, a * 2 (s))

) when s = s r , and q v (s) = q L (s) = 0 when s = s r ,

G π * =    π * (•|1) T . . . π * (•|m s ) T    , H π * V =    q π * V (1) T . . . q π * V (m s ) T    where q π * V (i) T = γ[π * (1|i)(V π * ) T , . . . , π * (j|i)(V π * ) T , . . . , π * (m a |i)(V π * ) T ],

which is an Ndimensional row vector, and H π * L and q π * L (i) are defined similarly by substituting V with L.

Case 1 in Theorem 3 corresponds to the case where the constraint in ( 8) is non-binding. This effectively reduces to the unconstrained scenario in Theorem 1 since a small perturbation of µ R , µ C , P does not affect feasibility. Case 2 is when the constraint is binding. In this case, α * must be chosen such that the split policy ensures equality in the constraint, and when µ R , µ C , P is perturbed the estimated α * n would adjust accordingly. Thus, in this case the estimation of V * incurs two sources of noises, one from the uncertainty of R, P that appears also in unconstrained problems, and one from the uncertainty in calibrating α * n that is in turn affected by C, P , thus leading to the extra terms in the variance expression.

A.2 Approximate Value Iteration:

When the state space S is large, updating an m s × m a look-up table via T µ R ,P (.) can be computationally infeasible. Approximate value iteration operates by applying a mapping M over T µ R ,P . In many cases, M = M g • M S0 I where M S0 I is a dimension-reducing "inherit" mapping R msma → R ms 0 ma , and M g is the "generalization" mapping R ms 0 ma → R msma that lifts back to the full dimension. By selecting a "representative" subset S 0 ⊂ S with cardinality

|S 0 | = m s0 m s , M S0 I is defined as M S0 I (x) = [x(i, j)] i∈S0,1≤j≤ma

where [x i ] i∈I denotes the set of entries of x whose index i ∈ I. In this setup, we define Q M as a fixed point of the operator M • T µ R ,P (•), and V M (s) = max a Q M (s, a). We also define

Q M S0 = M S0 I • T µ R ,P (Q M )

as the dimension-reduced Q-value. We derive large-sample error estimates in this case. For this, we first assume there is a welldefined metric on S. To guarantee the existence of Q M , we make the following assumption on the generalization map M g : Assumption 4. M g is a max-norm non-expansion mapping in S, i.e., ||M g

(x) -M g (y)|| ∞ < ||x -y|| ∞ ∀x, y ∈ S.

We also need the following analogs of Assumptions 2 and 3 to Q M and S 0 : Tables 6, 7 and 8 compare the CI coverage rates when the state space is large, i.e. m s = 31, using RE with different values of π(1|s), i.e., π(1|s) = 0.8, 0.85, and 0.9. Compared to exact update, the coverage convergence for approximate update appears generally slower. Specifically, comparing Tables 5 and 6 that use the same RE with π(1|s) = 0.8, we see that the coverages on the averages of Q's and V * 's for approximate update are only around 23% -25% when n = 10 4 , whereas they are 77% -78% for exact update. Also, while the nominal coverage, i.e., 95%, is obtained when n = 5 × 10 4 in the exact update for all studied quantities, this sample size is not enough for approximate update, where it appears we need n to be of order 10 7 to obtain the nominal coverage.

Furthermore, Tables 6, 7 and 8 show that the rates of convergence to the nominal coverage are quite different for different values of π(1|s)'s. The convergence rate when π(1|s) = 0.85 seems to be the fastest, with the coverage close to 95% already when n = 10 5 . On the other hand, when π(1|s) = 0.8, the coverage is close to 95% only when n is as large as 10 7 , and when π(1|s) = 0.9, even n = 10 7 is not enough for convergence to kick in. We also see that, when the coverage is very far from the nominal rate, discrepancies can show up among the estimates of Q, V * and χ * . For example, when π(1|s) = 0.8 and n = 10 4 , the coverages of Q and V * are around 23% -25% but the coverage of χ * is as low as 1%, and when π(1|s) = 0.9 and n = 10 7 , the coverages of Q and V * are around 75% -77% but that of χ * is only 29%. However, in settings where the coverage is close to 95%, all these quantities appear to attain this accuracy simultaneously in all considered cases. These caution that coverage accuracy can be quite sensitive to the specifications of the exploration policy. Nonetheless, the convergence behaviors predicted by Theorem 1, Corollary 1 and Theorem 4 are all observed to hold. 

B.2 Efficiency of Exploration Policies:

In Q-OCBA, our second-stage exploration policy is derived by maximizing the worst-case relative discrepancy among all Q-value estimates. If one is interested in obtaining the best estimate of χ * (i.e., the optimal value function initialized at a distribution ρ), then it would be more beneficial to consider solving min

w∈Wη σ 2 χ (9)

to derive the optimal second-stage exploration policy π w (recall Lemma 1). The motivation is that by doing so we would obtain a CI for χ * as short as possible.

Table 9 compares the 95% CI lengths and coverages for this exploration policy with other benchmark strategies, for r L ranging from 1 to 3. For each r L , we show the averages of the coverages and CI lengths of Q estimates among all (s, a) pairs, and also the coverage and CI length of χ * estimates. Note that our strategy intends to shorten the CI lengths of χ * estimates. Like our experiment in Section 5 in the main paper, we use a total observation budget n = 10 4 , and devote 30% to the initial stage where RE with π(1|s) = 0.8 is used to estimate the parameters used to plug in the criterion to be optimized in the second stage. For convenience and the consistency of terminology, we continue to call our procedure to attain criterion (9) Q-OCBA. We compare this with pure RE and -greedy, with ranging from 0.01 to 0.2.

Table 9 shows that our budget is enough to achieve the nominal 95% coverages for both the Q-values and χ * using all strategies, which is consistent with the conclusion from Theorem 1 and Corollary 1. However, Q-OCBA leads to desirably much shorter CI's generally, with the shortest CI lengths in all settings and sometimes by a big margin. For example, when r L = 2, the CI length derived by Q-OCBA is at least 80% less than those derived by all the other methods. We also observe that Q-OCBA performs much more stably than RE and -greedy, the latter varying quite significantly for different values of r L . When r L = 1, -greedy with = 0.01 can perform almost as well as Q-OCBA, with both CI lengths for Q being 2.45 -2.46 and for χ * being 2.41 -2.42. But when r L = 2, -greedy with the same = 0.01 cannot even explore all (s, a) pairs. The situation worsens when r L = 3, where none of the considered values of can explore all (s, a) pairs. This observation on -greedy is consistent with Table 3 in the main paper where we consider the criterion using the probability of correct selection. Regardless of using that or the current criteria, the performances of -greedy depend fundamentally on whether the (s, a) pairs that need to be explored more also tend to have larger Q-values. Note that when changing r L , the corresponding changes in the Q-values would change the exploration "preference" for -greedy. However, as the underlying stochasticity of the system does not change with r L , the states that need more exploration remain unchanged. This misalignment leads us to observe quite different performances for -greedy when r L varies. Lastly, again consistent with the results on the probability of correct selection shown in Table 3 of the main paper, we observe that Q-OCBA outperforms pure RE in all cases in Table 5 , with at least 40% shorter CI lengths for the χ * estimates. This is attributed to the efficient use of variance information in the second stage of Q-OCBA.

Table 9 : Length of CI comparison for different exploration policies a NA means that some (s, a) pair has never been visited.

= 0.2 = 0.1 = 0.01 π(1|s) = 0.8 Q-OCBA r L = 1 Q Coverage 0.

C Proofs of Main Results:

In this section, we present the proofs of the main results. In the proofs, we shall treat P as an N m s -dimensional vector following the index rule:

P ((i -1)N + (j -1)m s + k) = P (k|i, j) Proof of Theorem 1. Define F (Q , r , P ) as a mapping from R N × R N × R N ms to R N . Specifically, F (Q , r , P )((i -1)m a + j) = Q ((i -1)m a + j) -r ((i -1)m a + j) -γ 1≤k≤ms P ((i -1)N + (j -1)m s + k) g k (Q ) for 1 ≤ i ≤ m s and 1 ≤ j ≤ m a , where g k (Q ) = max l Q ((k -1)m a + l), for 1 ≤ k ≤ m s .

By Assumption 2, there exists an open neighborhood of Q, which we denote as Ω, such that ∀Q ∈ Ω, arg max j Q ((i -1)m a + j) is still unique for each 1 ≤ i ≤ m s . Then, for each 1 ≤ k ≤ m s , g k (Q ) has all its partial derivatives exist and continuous. This implies that

F (Q , r , P ) is continuously differentiable in Ω × R N × R N ms .

Denote the partial derivatives of F as ∂F ∂(Q , r , P ) = ∂F ∂Q ∂F ∂r

∂F:

∂P .

Note that ∂F ∂Q is an N ×N matrix. Denote its element at the ((i-1)m a +j)-th row, ((k -1)m a +l)-th column by

∂F (i-1)ma +j ∂Q (k-1)ma +l

. Then we have

∂F (i-1)ma+j ∂Q (k-1)ma+l = 1 (k = i, j = l) -γP ((i -1)N + (j -1)m s + k) ×1 Q ((k -1)m a + l) = max u Q ((k -1)m a + u) .

Putting all the elements together, we have ∂F ∂Q = I -γ P ,

where P is an N × N matrix with P ((i -1)m a + j, (k -1)m a + l) = P ((i -1)N + (j -1)m s + k)

×1(Q ((k -1)m a + l) = max u Q ((k -1)m a + u)), for 1 ≤ i ≤ m s , 1 ≤ j ≤ m a ,1 ≤ k ≤ m s ,1 ≤ l ≤ m a .

Since all rows of P sum up to one, P can be interpreted as the transition matrix of a Markov Chain with state space {(i, j) : 1 ≤ i ≤ m s , 1 ≤ j ≤ m a }. Note that ∂F ∂Q is invertible for any Q ∈ Ω. We can then apply the implicit function theorem to the equation F (Q, µ R , P ) = 0. In particular, there exists an open set U around µ R × P ∈ R N × R N ms , and a unique continuously differentiable function φ: U → R N , such that for any r × P ∈ U φ(µ R , P ) = Q F (φ(r , P ), r , P ) = 0. In addition, the partial derivatives of φ satisfy

∇φ(µ R , P ) := ∂φ ∂(r , P ) r =µ R ,P =P = - ∂F ∂Q -1 ∂F ∂r , ∂F ∂P Q =Q,r =µ R ,P =P

It is also easy to verify that ∂F ∂r Q =Q,r =µ R ,P =P = I N ×N .

We also note that

∂F (i-1)ma+j ∂P (k-1)N +(l-1)ms+v Q =Q,r =µ R ,P =P = γ max u Q((v -1)m a + u)1 (k = i, j = l) = γV * (v)1 (k = i, j = l) , for 1 ≤ i ≤ m s , 1 ≤ j ≤ m a ,1 ≤ k ≤ m s ,1 ≤ l ≤ m a , and 1 ≤ v ≤ m s . Then C π * := ∂F ∂P Q =Q,r =µ R ,P =P =         (V * ) T . . . (V * ) T . . . (V * ) T         , which is an N × N m s matrix. Next, for μR,n ((i -1)m a + j) = 1≤t≤n r t (s t , a t )1(s t = i, a t = j) 1≤t≤n 1(s t = i, a t = j) = 1≤t≤n r t (s t , a t )1(s t = i, a t = j) w((i -1)m a + j)n w((i -1)m a + j)n 1≤t≤n 1(s t = i, a t = j) and Pn ((i -1)N + (j -1)m s + k) = 1≤t≤n 1(s t = i, a t = j, s t (s t , a t ) = k) 1≤t≤n 1(s t = i, a t = j) = 1≤t≤n 1(s t = i, a t = j, s t (s t , a t ) = k) w((i -1)m a + j)n

w((i -1)m a + j)n 1≤t≤n 1(s t = i, a t = j) by Assumption 1, 3 and Slutsky's theorem, we have

EQUATION

where Σ R,P = W -1 D R 0 0 D P , and

D P =           Σ P 1,1 w(0ma+1)

. . .

Σ P i,j w((i-1)ma+j) . . . Σ Pm s,ma w((ms-1)ma+ma))          

, which is an N m s × N m s matrix. By the continuous mapping theorem, we have φ(μ R,n , Pn ) -φ(µ R , P ) → 0 a.s. as n → ∞, which implies Qn → Q a.s.. In addition, using the delta method, we have

√ n( Qn -Q) = √ n(φ(μ R,n , Pn ) -φ(µ R , P ))

⇒ N (0, ∇φ(µ R , P )Σ R,P ∇φ(µ R , P ) T ) as n → ∞.

We also have

∇φ(µ R , P )Σ R,P ∇φ(µ R , P ) T = (I -γ P π * ) -1 [I, C π * ] W -1 D R 0 0 D P [I, C π * ] T ((I -γ P π * ) -1 ) T = (I -γ P π * ) -1 (W -1 D R + C π * D P (C π * ) T )((I -γ P π * ) -1 ) T = (I -γ P π * ) -1 W -1 [D R + D Q ]((I -γ P π * ) -1 ) T . Proof of Corollary 1. Define g V (Q): R N → R ms as g V (Q) = (g 1 (Q), . . . , g ms (Q)) = (V π * (1), . . . , V π * (m s )),

which is continuously differentiable in an open neighborhood of Q. Then we can apply the delta method to get

√ n( V * n -V π * ) = √ n(g V ( Qn ) -g V (Q)) ⇒ N (0, ∇g V (Q)Σ(∇g V (Q)) T ) as n → ∞. Note that ∇g V (Q) is a m s × N matrix with ∇g V (Q)(i, (j -1)m a + k) = 1 (i = j, k = a * (i))

. By rearranging the index such that P π * = P π * 0 * 0 (where " * " denotes a placeholder of some quan- tities), we have ∇g

V (Q) = [I, 0], D R = Diag(D π * R , * ), W = Diag(W π * , * ) and D Q = Diag(D π * V , * ). We also note that ∇g V (Q)(I -γ P π * ) -1 = [I, 0] ∞ i=0 γ i ( P π * ) i = [I, 0] ∞ i=0 γ i (P π * ) i 0 * 0 = ∞ i=0 γ i (P π * ) i , 0 = (I -γP π * ) -1 , 0 . Thus, ∇g V (Q)Σ(∇g V (Q)) T = (I -γP π * ) -1 , 0 W -1 [D R + D Q ] (I -γP π * ) -1 , 0 T = (I -γP π * ) -1 (W π * ) -1 [D π * R + D π * V ]((I -γP π * ) -1 ) T .

Lastly, the asymptotic normality of χπ * n follows from the continuous mapping theorem.

We next establish the asymptotic normality for the estimated value function under a given policy π.

In this case, the value function V π satisfies a Bellman equation

V π (s) = a µ R (s, a)π(a|s) + γ a π(a|s) s P (s |s, a)V π (s ).

Denote the the estimator of V π by V π n . In particular, V π n is the fixed point of the corresponding empirical Bellman equation that replaces (µ R , P ) by (μ R,n , Pn ). We have:

Corollary 2. Under Assumptions 1 and 3, √ n( V π n -V π ) ⇒ N (0, Σ π V ) where Σ π V = X W X T , X = (I -γP π ) -1 , P π is an m s × m s transition matrix with P π (i, j) = a P (j|s = i, a)π(a|s = i), W is an m s × m s diagonal matrix with W (i, i) = j π(j|i) 2 w(i, j) [(γV π ) T Σ Pi,j (γV π ) + σ 2 R (i, j)]

Proof of Corollary 2. Similar to the proof of Theorem 1, define F π as a mapping from Note that F π (V π , µ R , P ) = 0, F π is continuously differentiable and I -γP π is invertible. We can thus apply the implicit function theorem. In particular, there exists an open set U π around µ R × P ∈ R N × R N ms , and a unique continuously differentiable function φ

EQUATION

π : U π → R N , such that φ π (µ R , P ) = V π F π (φ π (r , P ), r , P ) = 0

for any r × P ∈ U π . For the partial derivatives of φ π , we have

EQUATION

H π V := ∂F π ∂P V =V π ,r =µ R ,P =P =         (q π 1 ) T . . . (q π i ) T . . . (q π ms ) T         where (q π i ) T = γ[π(1|i)(V π ) T , . . . π(j|i)(V π ) T , . . . π(m a |i)(V π ) T ], which is an N -dimensional row vector. ûn = ( √ n( Pn -P ), √ n(μ R,n -µ R ))/ || √ n( Pn -P )|| 2 + || √ n(μ R,n -µ R )|| 2 is a continuous mapping of ( √ n( Pn -P ), √ n(μ R,n -µ R ))

. By taking n → ∞, we get the result by the continuous mapping theorem.

Proof of Theorem 3. We use the LP representation of the constrained MDP. Define x s,a as the occupancy measure

x s,a = ∞ t=0 γ t P ρ (S t = s)

where P ρ denotes the distribution of S t 's with initial distribution ρ. Then, x s,a satisfies the LP max s,a µ R (s, a)x s,a subject to s,a µ C (s, a)x s,a ≤ η a x s,a -γ s ,a P (s|s , a)x s ,a = ρ(s), ∀s x s,a ≥ 0, ∀s, a

(This is the dual formulation in the proof of Theorem 2 with an extra constraint.) The objective and the first constraint correspond to the objective and constraint in the constrained MDP formulation. The second constraint can be deduced by a one-step analysis on the definition of occupancy measure. Once ( 11) is solved to obtain an optimal solution (x * s,a ) s,a , it can be translated into an optimal policy π * (a|s) =

x * s,a a x * s,a

Note that the number of structural constraints in (11) is m s + 1, and a corner-point optimal solution has m s + 1 basic variables. Moreover, by our assumptions, the optimal solution is unique and the LP is non-degenerate, so that perturbing the parameters µ R (s, a), µ C (s, a), P (s|s , a) does not immediately imply an overshoot to negativity for the reduced costs of the non-basic variables. Now consider two cases depending on whether the first constraint is non-binding or binding. The first case corresponds to a deterministic optimal policy, i.e., for any s, x s,a > 0 for only one a. In this case a small perturbation of the parameters still retains the same basic and non-basic variables, and the derived perturbed policy still retains the non-binding first constraint. In this case, the analysis reduces back to Corollary 2.

In the second case, x s,a > 0 for only one a, for all s except one state s r , where we can have x sr,a * 1 (sr) > 0 and x sr,a * 2 (sr) > 0 for two distinct actions a * 1 (s r ), a * 2 (s r ). Again, perturbing the parameters retains these basic and non-basic variables. In particular, the first constraint remains binding in the perturbation, so that the perturbed optimal policy π * is still split at the same state and satisfies s ρ(s)L π * (s) = η. Now denote the mixing parameter by α * := π * (a * 1 (s r )|s r ), and so π * (a * 2 (s r )|s r ) = 1 -α * . By applying the implicit function theorem to the Bellman equation, there exists a continuously differentiable function φ L such that L π * (s) = φ L (µ C , P, α * ). By applying the implicit function theorem again to the equation s ρ(s)φ L (µ C , P, α * ) = η, we know α * is a continuously differentiable function of µ C and P . Thus V * can be viewed as a function of µ R , P and α * , the latter in turn depending on µ C and P . It can also be viewed as a function of µ R , µ C and P directly. We use ∇ µ R ,µ C ,P V * (µ R , µ C , P ) to denote the Jacobian of V * with respect to µ R , µ C , P when viewing V * as a function of these variables. We also use ∇ µ R ,µ C ,P,α V * (µ R , P, α * ) to denote the Jacobian of V * with respect to µ R , µ C , P, α * , this time viewing V * as a function of µ R , P, α * .

We use similar notations throughout, and their meanings should be clear from the context. To facilitate derivations, we also distinguish the notation ∇ x f , which denotes the multi-dimensional Jacobian matrix, from ∂ x f , which is used to denote the Jacobian when x is a scalar (1-dimensional). 

L π *               q 0 L = (I -γP π * ) -1 q L ,

where q 0 L is a vector with q 0 L (s r ) = µ R (s r , a * 1 (s r )) -µ R (s r , a * 2 (s r )), q 0 L (m s + (s r -1)m a + j) = P (j|s r , a * 1 (s r )) -P (j|s r , a * 2 (s r )) for 1 ≤ j ≤ m a and q 0 L (i) = 0 for any other index i. Similarly, we can define r R by r R (s) = ma j=1 µ R (s, j)π * (j|s), and we have

∂ α * V * (µ R , P, α * ) = (I -γP π * ) -1 q V

The derivation of ∇ µ R ,µ C ,P V * (µ R , P, α * ) and ∇ µ R ,µ C ,P L π * (µ C , P, α * ) follows exactly the same line of analysis as how we derive G π and H π V in the proof of Corollary 2.

Proof of Theorem 4. Denote Using the delta method, we have √ n( QM n -Q M ) ⇒ N (0, ∇φ M (µ R , P )Σ R,P,S (∇φ M (µ R , P )) T ) = N (0, Σ M S0 ).

Proof of lemma 1. For any given policy π, by the balance equation for Markov Chains, its induced stationary distribution w π satisfies k,l w π ((k -1)m a + l)P (i|s = k, a = l)π(a = j|s = i) = w π ((i -1)m a + j)

:



Figures and Tables:

s, a)π(a|s) -γ a π(a|s) s P (s |s, a)V (s ).

µ R ,µ C ,P V * (µ R , µ C , P ) = ∇ µ R ,µ C ,P,α * V * (µ R , P, α * )[I, ∇ µ R ,µ C ,P α * (µ C , P ) T ] T = ∇ µ R ,µ C ,P V * (µ R , P, α * ) + ∂ α * V * (µ R , P, α * )∇ µ R ,µ C ,P α * (µ C , P ) (12) Differentiating ρ T L π * (µ R , µ C , P ) = η, we have ρ T ∇ µ R ,µ C ,P L π * (µ C , P, α * ) + ρ T ∂ α * L π * (µ C , P, α * )∇ µ R ,µ C ,P α * (µ C , P ) = 0.By rearranging the equation, we have∇ µ R ,µ C ,P α * (µ C , P ) = -1 ρ T ∂ α * L π * (µ C , P, α * ) ρ T ∇ µ R ,µ C ,P L π * (µ C , P, α * )Substituting this into (12), we get∇ µ R ,µ C ,P V * (µ R , µ C , P ) = ∇ µ R ,µ C ,P V (µ R , P, α * ) -ρ T ∇ µ R ,µ C ,P L π * (µ C , P, α * )∂ α * V * (µ R , P, α * ) ρ T ∂ α * L π * (µ C , P, α * ) Next, define an m s -dimensional vector r C by r C (s) = ma j=1 µ C (s, j)π * (j|s). Then ∂ α * L π * (µ C , P, α * ) = ∇ r C ,P π * L π * (r C , P π * )[(∂ α * r C (α * )) T , (∂ α * P π * (α * )) T ] T Note that (I -γP π * )L π * =r C . By applying the implicit function theorem, we have ∂ α * L π * (µ C , P, α * ) = (I -γP π * )

[μ R,n , Pn ] S0 = [μ R,n ((i -1)m a + j), Pn ((i -1)N + (j -1)m s + k)] i∈S0,1≤j≤ma,1≤k≤ms , and [µ R , P ] S0 = [µ R ((i -1)m a + j), P ((i -1)N + (j -1)m s + k)] i∈S0,1≤j≤ma,1≤k≤ms . By Assumption 6, we have√ n([μ R,n , Pn ] S0 -[µ R , P ] S0 ) ⇒ N (0, ΣR,P,S0 ) where Σ R,P,S0 = Notice M g • M S0 I • T μR,n , Pn only involves random variables [μ R,n , Pn ] S0 . Changing the distribution of [μ R,n , Pn ] S\S0 will not change the distribution of QM n . We can thus assign auxiliary random variables to μR,n and Pn for alli / ∈ S 0 , 1 ≤ j ≤ m a , 1 ≤ k ≤ m s .In particular, we define independent random variables for each i / ∈ S 0 by lettingμR,n ((i -1)m a + j) R ((i -1)m a + j), 1) Pn ((i -1)N + (j -1)m s + k) ((i -1)N + (j -1)m s + k), 1).By doing so, we extend them s0 m a -dimensional random variable [μ R,n , Pn ] S0 to an m s m a -dimensional random variable [μ R,n , Pn ] S and √ n([μ R,n , Pn ] S -[µ R , P ] S ) ⇒ N (0, Σ R,P,S )where Σ R,P,S = Σ R,P,S0 0 0I .Similar to the proof of Theorem 1, defineF M (Q , r , P ) = Q -M g • M S0 I • T r ,P (Q ). By Assumption 4, M g is max-norm non-expansion. Then, M g • M S0I is also max-norm non-expansion, implying that ∇(M g • M S0I ) has all its eigenvalues less than or equal to 1. Thus,∂F M ∂Q = ∇M (T r ,P (Q ))(I -γ P )is invertible. By Assumption 4, we have F M (Q M , µ R , P ) = 0. By Assumption 5, there exists a neighborhood Ω M around (Q M , µ R , P ), such that F M is continuously differentiable on Ω M . Then, applying the implicit function theorem, we have that there exists an open set E M ⊂ Ω M and a continuously differentiable mapping φ M on E M , such that φ M (µ R , P ) = Q M and ∇φ M (µ R , P ) = -=Q M ,r =µ R ,P =P .

B k and ŵB k based on the B k data points collected ; Apply value-iteration using PB k and μ2 R,B k to obtain QB k ; Plug the estimates PB k , σ2 R,B k and QB k into (7) to solve for the optimal w k

Table content: <table/>

Exact tabular update

Table content: <table><tr><td>n</td><td>10 4</td><td>3 × 10 4</td><td>5 × 10 4</td></tr><tr><td>Q</td><td colspan="3">0.77(0.03) 0.93(0.02) 0.96(0.01)</td></tr><tr><td colspan="4">χ π  *  0.77(0.03) 0.93(0.02) 0.96(0.01)</td></tr></table>

Approximate value iteration

Table content: <table><tr><td>n</td><td>10 4</td><td>10 5</td><td>10 6</td></tr><tr><td>Q</td><td colspan="3">0.53(0.02) 0.95(0.01) 0.95(0.01)</td></tr><tr><td colspan="4">χ π  *  0.80(0.03) 0.94(0.02) 0.95(0.01)</td></tr></table>

Probability of correct selection for different exploration policies, n = 10 3

Table content: <table><tr><td colspan="2">r L 0.2-greedy</td><td>RE(0.6)</td><td>UCRL2</td><td>PSRL(100) Q-OCBA</td></tr><tr><td>1</td><td colspan="4">0.95(0.01) 0.70(0.03) 0.44(0.03) 0.53(0.03) 0.87(0.02)</td></tr><tr><td>2</td><td colspan="4">0.15(0.02) 0.29(0.03) 0.11(0.02) 0.33(0.03) 0.55(0.03)</td></tr><tr><td>3</td><td colspan="4">0.00(0.00) 0.45(0.03) 0.21(0.02) 0.41(0.03) 0.84(0.02)</td></tr></table>

Probability of correct selection for different exploration policies, n = 10 4 Shie Mannor, Duncan Simester, Peng Sun, and John N Tsitsiklis. Bias and variance in value function estimation. In Proceedings of the twenty-first international conference on Machine learning, pp. 72. ACM, 2004. Shie Mannor, Duncan Simester, Peng Sun, and John N Tsitsiklis. Bias and variance approximation in value function estimates. Management Science, 53(2):308-322, 2007. Rémi Munos and Csaba Szepesvári. Finite-time bounds for fitted value iteration. Journal of Machine Learning Research, 9(May):815-857, 2008. Ian Osband, Daniel Russo, and Benjamin Van Roy. (More) efficient reinforcement learning via posterior sampling. In Advances in Neural Information Processing Systems, pp. 3003-3011, 2013. Robert J Serfling. Approximation Theorems of Mathematical Statistics, volume 162. John Wiley & Sons, 2009.

Table content: <table><tr><td colspan="2">r L 0.2-greedy</td><td>RE(0.6)</td><td>UCRL2</td><td>PSRL(100) Q-OCBA</td></tr><tr><td>1</td><td colspan="4">1.00(0.00) 0.95(0.01) 0.82(0.02) 1.00(0.00) 1.00(0.00)</td></tr><tr><td>2</td><td colspan="4">0.55(0.03) 0.80(0.03) 0.52(0.03) 0.94(0.02) 1.00(0.00)</td></tr><tr><td>3</td><td colspan="4">0.21(0.03) 0.94(0.01) 0.75(0.03) 0.76(0.03) 1.00(0.00)</td></tr></table>

Coverage for Q(s, a), V * and χ * values using exact tabular update

Table content: <table><tr><td>n</td><td>10 4</td><td>3 × 10 4</td><td>5 × 10 4</td></tr><tr><td>Q(1, 0)</td><td colspan="3">0.773(0.026) 0.930(0.016) 0.956(0.013)</td></tr><tr><td>Q(1, 1)</td><td colspan="3">0.773(0.026) 0.930(0.016) 0.956(0.013)</td></tr><tr><td>Q(2, 0)</td><td colspan="3">0.773(0.026) 0.930(0.016) 0.956(0.013)</td></tr><tr><td>Q(2, 1)</td><td colspan="3">0.771(0.026) 0.930(0.016) 0.955(0.013)</td></tr><tr><td>Q(3, 0)</td><td colspan="3">0.771(0.026) 0.930(0.016) 0.955(0.013)</td></tr><tr><td>Q(3, 1)</td><td colspan="3">0.771(0.026) 0.929(0.016) 0.957(0.013)</td></tr><tr><td>Q(4, 0)</td><td colspan="3">0.771(0.026) 0.929(0.016) 0.957(0.012)</td></tr><tr><td>Q(4, 1)</td><td colspan="3">0.775(0.026) 0.934(0.015) 0.958(0.012)</td></tr><tr><td>Q(5, 0)</td><td colspan="3">0.775(0.026) 0.934(0.015) 0.958(0.012)</td></tr><tr><td>Q(5, 1)</td><td colspan="3">0.775(0.026) 0.935(0.015) 0.954(0.013)</td></tr><tr><td>Q(6, 0)</td><td colspan="3">0.775(0.026) 0.935(0.015) 0.954(0.013)</td></tr><tr><td>Q(6, 1)</td><td colspan="3">0.750(0.026) 0.920(0.017) 0.950(0.014)</td></tr><tr><td colspan="4">Average of Q 0.771(0.026) 0.931(0.016) 0.956(0.013)</td></tr><tr><td>V  *  (1)</td><td colspan="3">0.773(0.026) 0.930(0.016) 0.956(0.013)</td></tr><tr><td>V  *  (2)</td><td colspan="3">0.772(0.026) 0.930(0.016) 0.955(0.013)</td></tr><tr><td>V  *  (3)</td><td colspan="3">0.774(0.026) 0.929(0.016) 0.957(0.013)</td></tr><tr><td>V  *  (4)</td><td colspan="3">0.776(0.026) 0.934(0.015) 0.958(0.012)</td></tr><tr><td>V  *  (5)</td><td colspan="3">0.772(0.026) 0.935(0.015) 0.954(0.013)</td></tr><tr><td>V  *  (6)</td><td colspan="3">0.752(0.027) 0.920(0.017) 0.950(0.014)</td></tr><tr><td colspan="4">Average of V  *  0.769(0.026) 0.929(0.016) 0.955(0.012)</td></tr><tr><td>χ  *</td><td colspan="3">0.772(0.026) 0.933(0.015) 0.958(0.012)</td></tr></table>

Linear interpolation in approximate value iteration with π(1|s) = 0.8

Table content: <table><tr><td>n</td><td>10 4</td><td>10 5</td><td>10 6</td><td>10 7</td></tr><tr><td colspan="5">Average of Q 0.23(0.01) 0.55(0.03) 0.92(0.02) 0.95(0.01)</td></tr><tr><td colspan="5">Average of V  *  0.25(0.01) 0.56(0.02) 0.93(0.01) 0.95(0.01)</td></tr><tr><td>χ  *  coverage</td><td colspan="4">0.01(0.01) 0.35(0.03) 0.94(0.02) 0.94(0.02)</td></tr><tr><td colspan="5">Table 7: Linear interpolation in approximate value iteration with π(1|s) = 0.85</td></tr><tr><td>n</td><td>10 4</td><td>10 5</td><td>10 6</td><td>10 7</td></tr><tr><td colspan="5">Average of Q 0.53(0.02) 0.95(0.01) 0.95(0.01) 0.95(0.01)</td></tr><tr><td colspan="5">Average of V  *  0.59(0.02) 0.95(0.01) 0.95(0.01) 0.95(0.01)</td></tr><tr><td>χ  *  coverage</td><td colspan="4">0.80(0.03) 0.94(0.02) 0.95(0.01) 0.94(0.01)</td></tr></table>

Linear interpolation in approximate value iteration with π(1|s) = 0.9

Table content: <table><tr><td>n</td><td>10 4</td><td>10 5</td><td>10 6</td><td>10 7</td></tr><tr><td colspan="5">Average of Q 0.31(0.01) 0.49(0.01) 0.67(0.01) 0.75(0.02)</td></tr><tr><td colspan="5">Average of V  *  0.37(0.01) 0.54(0.01) 0.71(0.01) 0.77(0.02)</td></tr><tr><td>χ  *  coverage</td><td colspan="4">0.25(0.03) 0.30(0.03) 0.28(0.02) 0.29(0.03)</td></tr></table>

Peer review (iclr2020):
This paper studies the inference problem of reinforcement learning. With a given exploration policy that satisfies some strong property to collect n data, the paper studies the distribution of the estimated optimal value function and Q-function when n goes to infinity. Both unique and non-unique optimal policy cases are studied. The non-unique case has a very different behavior as it is no longer Gaussian. The paper then uses these estimations to design a method that better explores and proposes a method Q-OCBA. Experiments were performed to compare this method with previous algorithms, e.g., UCRL. The inference results for the Q-function is definitely important and will be useful for the community. But the authors should not ignore all the finite-sample bound results in their related work section, especially for those who achieving minimax rates. For instance: MG Azar, I Osband, R Munos, Minimax regret bounds for reinforcement learning, 2017 MG Azar, R Munos, H Kappen, Minimax PAC bounds on the sample complexity of reinforcement learning with a generative model, 2013 A Sidford, M Wang, X Wu, L Yang, Y Ye, Near-optimal time and sample complexities for solving Markov decision processes with a generative model, 2018 A Agarwal, S Kakade, L Yang, On the Optimality of Sparse Model-Based Planning for Markov Decision Processes Please also see references therein. You can derive an inference result for the generative model as well. More comments: * The assumption on exploration policy is quite restrictive. It essentially says the bound does not hold if the MDP is not communicating. However, the Q-function is still well defined. * The comparison to UCRL and other exploration-based method is not fair. In these results, the probability of picking *the* optimal policy is not important. Also, picking the correct action is not feasible if the problem is more complex. It is my understanding that your algorithm has a much larger regret due to that it requires to sample from every state. * What about finite horizon MDP?