Efficient Inference and Exploration for Rein- forcement Learning
Year: 2019

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.

heading: 1 Introduction
text: 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).

heading: 2 Problem Setup
text: 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, . . . ,ms} and A = {1, 2, . . . ,ma}. Finally, we make the following stochasticity assumption: Assumption 1. R(s, a) has finite mean µR(s, a) and finite variance σ2R(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 madimensional 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 γ
tR(st, π(st))|s0 = s] with st+1 ∼ P (.|st, π(st)). 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 V ∗ = V π∗ and χ∗ = ∑ s ρ(s)V
∗(s). The Q-value, denoted by Q(s, a), is defined as Q(s, a) = µR(s, a) +γE[V ∗(S′)|s, a]. Correspondingly, V ∗(s) = maxaQ(s, a) and the Bellman equation for Q takes the form
Q(s, a) = µR(s, a) + γE [ max a′ Q(s′, a′)|s, a ] , (1)
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 maxa∈AQ(s, a) is unique.
Under Assumption 2, the optimal policy π∗ is unique and deterministic. Let a∗(s) = arg maxa∈AQ(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 {(st, at, rt(st, at), s′t(st, at)) : 1 ≤ t ≤ n}, where rt(st, at) is the realized reward at time t and s′t(st, at) = st+1. We define the sample mean µ̂R,n and the sample variance σ̂2R,n of the reward as
µ̂R,n(s = i, a = j) =
∑ 1≤t≤n rt(st, at)1(st = i, at = j)∑
1≤t≤n 1(st = i, at = j) , (2)
σ̂2R,n(s = i, a = j) =
∑ 1≤t≤n rt(st, at)
21(st = i, at = j)∑ 1≤t≤n 1(st = i, at = j) − µ̂R,n(i, j)2. (3)
Similarly, we define the empirical transition matrix P̂n as
P̂n(s ′ = k|s = i, a = j) =
∑ 1≤t≤n 1(st = i, at = j, s
′ t(st, at) = k)∑
1≤t≤n 1(st = i, at = j) (4)
and its ms ×ms sampling covariance matrix ΣPs,a (with one sample point of 1(st = s, at = a)) as
ΣPs,a(k1, k2) = { P (k1|s, a)(1− P (k1|s, a)) k1 = k2 −P (k1|s, a)P (k2|s, a) k1 6= k2. , for 1 ≤ k1 ≤ ms, 1 ≤ k2 ≤ ms.
With the data, we construct our estimate of Q, called Q̂n, which is the empirical fixed point of Tµ̂R,n,P̂n , i.e. Q̂n = Tµ̂R,n,P̂n(Q̂n). Correspondingly, we also write V̂ ∗ n (s) = maxa∈A Q̂n(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 Q̂n can be solved exactly. .

heading: 3 Quantifying Asymptotic Estimation Errors
text: We present an array of results regarding the asymptotic behaviors of Q̂n 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) = limn→∞ 1n ∑ 1≤t≤n 1(st = i, at = 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 ei as the i-th unit vector. The dimension of N (µ,Σ), I and ei should be clear from the context. When not specified, all the vectors are column vectors. Let N = msma. 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)ma + j, e.g. µR(i, j) = µR((i− 1)ma + 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)ma + j, (i′ − 1)ma + j′).

heading: 3.1 Limit Theorems under Sufficient Exploration
text: We first establish the asymptotic normality of Q̂n under exploration policy π: Theorem 1. Under Assumptions 1 and 2, if the data is collected according to π satisfying Assumption 3, then Q̂n is a strongly consistent estimator of Q, i.e. Q̂n → Q almost surely as n→∞. Moreover,
√ n(Q̂n −Q)⇒ N (0,Σ) as n→∞,
where Σ = (I − γP̃π ∗ )−1W−1(DR +DQ)((I − γP̃π ∗ )−1)T , (5)

heading: W , DR and DQ are N ×N diagonal matrices with
text: W ((i− 1)ma + j, (i− 1)ma + j) = w(i, j), DR((i− 1)ma + j, (i− 1)ma + j) = σ2R(i, j) and DQ((i− 1)ma + j, (i− 1)ma + j) = (V ∗)TΣPi,jV ∗ 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̃π∗)−1V ∗. The variances of these parameter estimates (i.e., (2) and (4)) involve σ2R(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) = maxa∈AQ(s, a) and V̂ ∗n (s) = maxa∈A Q̂n(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 ms ×ms transition matrix with Pπ ∗ (i, j) = P (j|s = i, a = a∗(s)), Wπ∗ , Dπ ∗
R and D π∗ V are ms ×ms 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ΣPi,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, a1)−Q(s, a2), i.e. the difference between action a1 and a2 when the agent is in state s. Define σ2∆Q as
σ2∆Q(s, a1, a2) = (e(s−1)ma+a1 − e(s−1)ma+a2) TΣ(e(s−1)ma+a1 − e(s−1)ma+a2) (6)
and its estimator σ̂2∆Q,n by replacing Q, V ∗, σ2R,n, w, P with Q̂n,V̂ ∗ n σ̂ 2 R,n, ŵn, P̂n in Σ, where ŵn is the empirical frequency of visiting each state-action pair, i.e. ŵn(i, j) = 1n ∑ 1≤t≤n 1(st =
i, at = j). Then the 100(1 − α)% confidence interval (CI) for Q(s, a1) − Q(s, a2) takes the form( Q̂n(s, a1)− Q̂n(s, a2) ) ± zασ̂2∆Q,n(s, a1, a2), where zα is the (1− α/2)-quantile of N (0, 1).

heading: 3.2 Non-Unique Optimal Policy
text: Suppose the optimal policy for the MDPM is not unique, i.e., Assumption 2 does not hold. In this situation, the estimated Q̂n 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 ms×(Nms+N) matrices {Gk}1≤k≤K and a deterministic partition of U = {u ∈ RmsN+ms : ||u|| = 1} = ∪1≤k≤KUk such that √ n(V̂ ∗n − V ∗)⇒ ∑K k=1Gk1 (Z/||Z|| ∈ Uk)Z, where Z = N (0,ΣR,P ), ΣR,P = Diag(W−1DR, DP ) and DP = Diag(ΣP1,1/w(0ma + 1), . . . ,ΣPi,j/w((i− 1)ma + j), . . . ,ΣPms,ma /w((ms − 1)ma +ma)).
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).

heading: 4 Efficient Exploration Policy
text: 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 6= a∗(i), the relative discrepancy as
hij = (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 hij ’s,
max w∈Wη min i∈S min j∈A,j 6=a∗(i) hij , (7)
where w denotes the proportions of visits on the state-action pairs, within some allocation set Wη (which we will explain). Intuitively, hij 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 hij 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 hij . 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 σ2R that appear in the expression of hij . 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)ma + j) = ∑ 1≤k≤ms ∑ 1≤l≤ma w((k − 1)ma + l)P (i|k, l)
∀1 ≤ i ≤ ms, ∑ 1≤i≤ms ∑ 1≤j≤ma w((i− 1)ma + 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)ma + j)/ ( ∑ma k=1 w((i− 1)ma + 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 minw∈Wη σ2χ. We showcase this with additional experiments in the Appendix.
Input: Number of iterations K, length of each batch {Bk}1≤k≤K , initial exploration policy π0; Initialization: k = 0; while k ≤ K do
Run πk for Bk steps and set k = k + 1; Calculate P̂Bk , µ̂R,Bk ,σ̂2R,Bk and ŵBk based on the Bk data points collected ; Apply value-iteration using P̂Bk and µ̂2R,Bk to obtain Q̂Bk ; Plug the estimates P̂Bk , σ̂2R,Bk and Q̂Bk into (7) to solve for the optimal wk ; Set πk(a = j|s = i) = wk((i− 1)ma + j)/ ∑ma l=1 wk((i− 1)ma + l);
end Algorithm 1: Q-OCBA sequential updating rule for exploration Note that (7) is equivalent to minw maxi∈S maxj∈A,j 6=a∗(i) ∑ s,a cij(s, a)/ws,a subject to w ∈ Wη, where cij(s, a)’s are non-negative coefficients. Based on the closed-form characterization of Σ in Theorem 1, cij(s, a)’s can be estimated with plug-in estimators using data collected in earlier stages.

heading: 5 Numerical Experiments
text: 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 ms 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 rL will be varied). We consider the infinite horizon setting with γ = 0.95 and ρ = [1/ms, . . . , 1/ms]T .
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 ms = 6 (using Theorem 1 and Corollary 1) and a large ms = 31 (using Theorem 4 in the Appendix) respectively. The latter case uses a linear interpolation with S0 = {1, 4, . . . , 28, 31}. All coverage rates are estimated using 103 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× 104 for exact update and ≥ 105 for interpolation), we see highly accurate CI coverages, i.e., close to 95%.
Table 1: Exact tabular update
n 104 3× 104 5× 104 Q 0.77(0.03) 0.93(0.02) 0.96(0.01)
χπ ∗ 0.77(0.03) 0.93(0.02) 0.96(0.01)
Table 2: Approximate value iteration
n 104 105 106
Q 0.53(0.02) 0.95(0.01) 0.95(0.01)
χπ ∗ 0.80(0.03) 0.94(0.02) 0.95(0.01)
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 ms = 6 and vary rL 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 Q̂n’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 103 replications of the procedure. We observe that Q-OCBA substantially outperforms the other methods, both with a small data size (n = 103 in Table 3) and a larger one (n = 104 in Table 4). Generally, these benchmark policies perform worse for larger values of rL. This is because for small rL, the (s, a) pairs that need to be explored more also tend to have larger Q-values. However, as rL 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.

heading: A Additional Theoretical Results
text: 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 γ
tR(st, π(st))|s0 = s], while at the same time want to ensure that a long-run accumulated discounted cost, denoted as Lπ(s) = E[ ∑∞ t=0 γ
tC(st, π(st))|s0 = s] which we call the loss function, is constrained by some given value η, i.e.,
max π ∑ s ρ(s)V π(s) subject to ∑ s ρ(s)Lπ(s) ≤ η (8)
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, P̂n. 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 sr, is visited, action a∗1(sr) is chosen with probability α∗ and action a∗2(sr) 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.

heading: 2. The optimal policy is deterministic, except at one state where a randomization between two
text: actions occurs, with the mixing parameter α∗. Denote the state where the randomization occurs by sr and two possible actions for sr by a∗1(sr) and a∗2(sr), We have
Σ = ( (I − γPπ ∗ )−1[Gπ ∗ , 0, Hπ ∗ V ]− (I − γPπ∗)−1qV ρT (I − γPπ ∗ )−1[0, Gπ ∗ , Hπ ∗ L ]
ρT (I − γPπ∗)−1qL ) ΣR,C,P ( (I − γPπ ∗ )−1[Gπ ∗ , 0, Hπ ∗ V ]− (I − γPπ∗)−1qV ρT (I − γPπ ∗ )−1[0, Gπ ∗ , Hπ ∗ L ]
ρT (I − γPπ∗)−1qL
)T ,
where ΣR,C,P = Diag(W−1DR,W−1DC , DP ), qV and qL are ms-dimensional vectors with qV (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 qL(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 = sr, and qv(s) = qL(s) = 0 when s 6= sr,
Gπ∗ = π ∗(·|1)T
. . . π∗(·|ms)T
 , Hπ∗V = q π∗ V (1) T
. . . qπ ∗
V (ms) T  where qπ ∗
V (i) T = γ[π∗(1|i)(V π∗)T , . . . , π∗(j|i)(V π∗)T , . . . , π∗(ma|i)(V π ∗ )T ], which is an N-
dimensional 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 ms×ma 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 = Mg ◦MS0I where M S0 I is a dimension-reducing “inherit” mapping Rmsma → Rms0ma , and Mg is the “generalization” mapping Rms0ma → Rmsma that lifts back to the full dimension. By selecting a “representative” subset S0 ⊂ S with cardinality |S0| = ms0 ms, M S0 I is defined as MS0I (x) = [x(i, j)]i∈S0,1≤j≤ma where [xi]i∈I denotes the set of entries of x whose index i ∈ I. In this setup, we define QM as a fixed point of the operator M ◦ TµR,P (·), and VM (s) = maxaQM (s, a). We also define QMS0 = M S0 I ◦ TµR,P (QM ) 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 QM , we make the following assumption on the generalization map Mg: Assumption 4. Mg is a max-norm non-expansion mapping in S, i.e., ||Mg(x) −Mg(y)||∞ < ||x− y||∞ ∀x, y ∈ S.
We also need the following analogs of Assumptions 2 and 3 to QM and S0:
Assumption 5. For any state s ∈ S, arg maxa∈AQM (s, a) is unique.
Assumption 6. For the Markov Chain with transition probability P̃π, the set of states {(s, a) : s ∈ S0, a ∈ A} is in the same communication class and this class is positive recurrent.
Let N0 = ms0ma and IS0 = {(i − 1)ma + j : i ∈ S0, a ∈ A}. With Assumption 6, we denote P̃MS0 as a sub-matrix of the matrix P̃π that only contains rows with indexes in IS0 . We also denote S0(i) as the i-th element (state) in S0. We define Q̂Mn as the empirical estimator of QM built on n observations. Then we have:
Theorem 4. Under Assumption 4, 5 and 6, if Mg is continuously differentiable, then √ n(Q̂Mn −

heading: QM )⇒ N (0,ΣMS0) as n→∞, where
text: ΣMS0 = (I − γ∇Mg(Q M S0)P̃ M S0 ) −1∇Mg(QMS0)(W S0)−1[DS0R +D S0 Q ]
∇Mg(QMS0) T ((I − γ∇Mg(QMS0)P̃ M S0 ) −1)T ,
∇Mg is the Jacobian of mapping Mg, WS0 , DS0R , D S0 Q are N0×N0 diagonal matrices with WS0((i− 1)ma + j, (i− 1)ma + j) = w(S0(i), j), DS0R ((i− 1)ma + j, (i− 1)ma + j) = σ2R(S0(i), j), D S0 Q ((i− 1)ma + j, (i− 1)ma + j) = (VM )TΣPS0(i),jV M .
Assumption 4 is generally satisfied by “local” approximation methods such as linear interpolation, k-nearest neighbors and local weighted average (Gordon, 1995). In all these cases, ∇Mg in Theorem 4 is actually a constant matrix.

heading: B Additional Numerical Results
text: This section reports additional numerical experiments. Sections B.1 present further results on the estimation quality of Q-values, V ∗and χ∗. Section B.2 provides additional results to demonstrate the efficiency of our proposed exploration strategy.
B.1 Statistical Quality of Interval Estimators
In this section, we provide additional numerical results about Tables 1 and 2 in the main paper. For Q-values in Table 1 of main paper, we only report the average coverage rate over all (s, a) pairs. Table 5 presents the coverage rates for each individual Q-values. It also provides the coverage rates for the value functions and χ∗, where χ∗ is the uniformly initialized value function, i.e. ρ = [1/ms, . . . , 1/ms]
T . We use RE with π(1|s) = 0.5 as our exploration policy. We see that the behaviors of these individual estimates are largely consistent. The coverage rates all converge to the nominal 95% as the number of observations n increases. Moreover, the coverages for the individual Q’s, V ∗’s, and the averages of these quantities are similar at any given n. Specifically, when n = 104, the coverages are all around 77% − 78%, when n = 3 × 104 they are all around 93%, and when n = 5 × 104 they are all very close to 95%. These suggest a sample size of 5 × 104 (or lower) is enough to elicit our asymptotic results in Theorem 1 and Corollary 1 in this problem.
Tables 6, 7 and 8 compare the CI coverage rates when the state space is large, i.e. ms = 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 = 104, whereas they are 77% − 78% for exact update. Also, while the nominal coverage, i.e., 95%, is obtained when n = 5 × 104 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 107 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 = 105. On the other hand, when π(1|s) = 0.8, the coverage is close to 95% only when n is as large as 107, and when π(1|s) = 0.9, even n = 107 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 = 104, 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 = 107, 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 rL ranging from 1 to 3. For each rL, 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 = 104, 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 rL = 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 rL. When rL = 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 rL = 2, -greedy with the same = 0.01 cannot even explore all (s, a) pairs. The situation worsens when rL = 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 rL, 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 rL, the states that need more exploration remain unchanged. This misalignment leads us to observe quite different performances for -greedy when rL 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.

heading: C Proofs of Main Results
text: In this section, we present the proofs of the main results. In the proofs, we shall treat P as an Nms-dimensional vector following the index rule: P ((i− 1)N + (j − 1)ms + k) = P (k|i, j)
Proof of Theorem 1. Define F (Q′, r′, P ′) as a mapping from RN × RN × RNms to RN . Specifically, F (Q′, r′, P ′)((i− 1)ma + j) = Q′((i− 1)ma + j)− r′((i− 1)ma + j) −γ
∑ 1≤k≤ms P ′((i− 1)N + (j − 1)ms + k) gk(Q′)
for 1 ≤ i ≤ ms and 1 ≤ j ≤ ma, where gk(Q′) = maxlQ′((k − 1)ma + l), for 1 ≤ k ≤ ms.
By Assumption 2, there exists an open neighborhood of Q, which we denote as Ω, such that ∀Q′ ∈ Ω, arg maxj Q
′((i− 1)ma + j) is still unique for each 1 ≤ i ≤ ms. Then, for each 1 ≤ k ≤ ms, gk(Q′) has all its partial derivatives exist and continuous. This implies that F (Q′, r′, P ′) is continuously differentiable in Ω× RN × RNms . 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)ma+j)-th row, ((k−1)ma+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)ms + k)
×1 ( Q′((k − 1)ma + l) = max
u Q′((k − 1)ma + u)
) .
Putting all the elements together, we have ∂F
∂Q′ = I − γP̃ ′,
where P̃ ′ is an N ×N matrix with P̃ ′((i− 1)ma + j, (k − 1)ma + l) = P ′((i− 1)N + (j − 1)ms + k)
×1(Q′((k − 1)ma + l) = max u Q′((k − 1)ma + u)),
for 1 ≤ i ≤ ms, 1 ≤ j ≤ ma,1 ≤ k ≤ ms,1 ≤ l ≤ ma.
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 ≤ ms, 1 ≤ j ≤ ma}. 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 ∈ RN × RNms , and a unique continuously differentiable function φ: U → RN , 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 = IN×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)ma + u)1 (k = i, j = l)
= γV ∗(v)1 (k = i, j = l) ,
for 1 ≤ i ≤ ms, 1 ≤ j ≤ ma,1 ≤ k ≤ ms,1 ≤ l ≤ ma, and 1 ≤ v ≤ ms. Then
Cπ ∗ := ∂F
∂P ′ ∣∣∣∣ Q′=Q,r′=µR,P ′=P =  (V ∗)T . . . (V ∗)T
. . . (V ∗)T  , which is an N ×Nms matrix. Next, for
µ̂R,n((i− 1)ma + j) = ∑
1≤t≤n rt(st, at)1(st = i, at = j)∑ 1≤t≤n 1(st = i, at = j)
=
∑ 1≤t≤n rt(st, at)1(st = i, at = j)
w((i− 1)ma + j)n w((i− 1)ma + j)n∑ 1≤t≤n 1(st = i, at = j)
and
P̂n((i− 1)N + (j − 1)ms + k)
=
∑ 1≤t≤n 1(st = i, at = j, s
′ t(st, at) = k)∑
1≤t≤n 1(st = i, at = j)
=
∑ 1≤t≤n 1(st = i, at = j, s ′ t(st, at) = k)
w((i− 1)ma + j)n w((i− 1)ma + j)n∑ 1≤t≤n 1(st = i, at = j)
by Assumption 1, 3 and Slutsky’s theorem, we have
[µ̂R,n, P̂n]− [µR, P ] −→ 0 a.s and
√ n([µ̂R,n, P̂n]− [µR, P ])⇒ N (0,ΣR,P ) (10)
where ΣR,P = ( W−1DR 0
0 DP
) , and
DP =  ΣP1,1 w(0ma+1) . . . ΣPi,j w((i−1)ma+j) . . .
ΣPms,ma w((ms−1)ma+ma))
 ,
which is an Nms ×Nms matrix. By the continuous mapping theorem, we have
φ(µ̂R,n, P̂n)− φ(µR, P )→ 0 a.s. as n→∞,
which implies Q̂n → Q a.s.. In addition, using the delta method, we have √ n(Q̂n −Q) = √ n(φ(µ̂R,n, P̂n)− φ(µ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−1DR 0
0 DP
) [I, Cπ ∗ ]T ((I − γP̃π ∗ )−1)T
= (I − γP̃π ∗ )−1(W−1DR + C π∗DP (C π∗)T )((I − γP̃π ∗ )−1)T = (I − γP̃π ∗ )−1W−1[DR +DQ]((I − γP̃π ∗ )−1)T .
Proof of Corollary 1. Define gV (Q): RN → Rms as
gV (Q) = (g1(Q), . . . , gms(Q)) = (V π∗(1), . . . , V π ∗ (ms)),
which is continuously differentiable in an open neighborhood of Q. Then we can apply the delta method to get
√ n(V̂ ∗n − V π ∗ ) = √ n(gV (Q̂n)− gV (Q))⇒ N (0,∇gV (Q)Σ(∇gV (Q))T ) as n→∞.
Note that ∇gV (Q) is a ms ×N matrix with ∇gV (Q)(i, (j − 1)ma + k) = 1 (i = j, k = a∗(i)). By
rearranging the index such that P̃π ∗ = ( Pπ ∗ 0 ∗ 0 ) (where “∗” denotes a placeholder of some quantities), we have ∇gV (Q) = [I, 0], DR = Diag(Dπ ∗ R , ∗), W = Diag(Wπ ∗ , ∗) and DQ = Diag(Dπ ∗
V , ∗). We also note that
∇gV (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, ∇gV (Q)Σ(∇gV (Q))T = [ (I − γPπ ∗ )−1,0 ] W−1[DR +DQ] [ (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, P̂n). 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 ms ×ms transition matrix with P π̃(i, j) = ∑ a P (j|s = i, a)π̃(a|s = i),

heading: W ′ is an ms ×ms diagonal matrix with
text: W ′(i, i) = ∑ j π̃(j|i)2 w(i, j) [(γV π̃)TΣPi,j (γV π̃) + σ2R(i, j)]
Proof of Corollary 2. Similar to the proof of Theorem 1, define F π̃ as a mapping from Rms ×RN × RNms → Rms :
F π̃(V ′, r′, P ′)(s) = V ′(s)− ∑ a r′(s, a)π̃(a|s)− γ ∑ a π̃(a|s) ∑ s′ P ′(s′|s, a)V ′(s′).
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 ∈ RN × RNms , and a unique continuously differentiable function φπ̃: U π̃ → RN , such that
φπ̃(µR, P ) = V π̃
F π̃(φπ̃(r′, P ′), r′, P ′) = 0
for any r′ × P ′ ∈ U π̃. For the partial derivatives of φπ̃, we have
∂φπ̃
∂(r′, P ′) ∣∣∣∣ r′=µR,P ′=P = − [ ∂F π̃ ∂V ′ ]−1 [ ∂F π̃ ∂r′ , ∂F π̃ ∂P ′ ]∣∣∣∣∣ V ′=V π̃,r′=µR,P ′=P
where ∂F π̃
∂V ′ ∣∣∣∣ V ′=V π̃,r′=µR,P ′=P = I − γP π̃,
Gπ̃ := ∂F π̃
∂r′ ∣∣∣∣ V ′=V π̃,r′=µR,P ′=P =  π̃(·|1)T . . . π̃(·|i)T
. . . π̃(·|ms)T  , and
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 , . . . π̃(ma|i)(V π̃)T ], which is an N -dimensional row vector.
Applying the delta method, we have √ n(V̂ π̃n − V π̃) = √ n(φπ̃(µ̂R,n, P̂n)− φπ̃(µR, P ))
⇒ N (0,∇φπ̃(µR, P )ΣR,P∇φπ̃(µR, P )T ) as n→∞, where
∇φπ̃(µR, P )ΣR,P∇φπ̃(µR, P )T = (I − γP π̃)−1[Gπ̃, H π̃V ] ( W−1DR 0
0 DP
) [Gπ̃, H π̃V ] T ((I − γP π̃)−1)T
= (I − γP π̃)−1(Gπ̃W−1DR(Gπ̃)T ) +H π̃VDP (H π̃V )T )((I − γP π̃)−1)T
and the conclusion follows.
Proof of Theorem 2. Write the MDP problem in its LP representation
max ∑ s ρ(s)V (s)
subject to V (s) ≥ r(s, a) + γ ∑ s′∈S P (s ′|s, a)V (s′), ∀s, a
with the dual problem max
∑ s,a µR(s, a)xs,a
subject to ∑ a xs,a − γ ∑ s′,a P (s|s′, a)xs′,a = ρ(s), ∀s
xs,a ≥ 0, ∀s, a
The decision variables in the dual problem, xs,a’s, in particular represent the occupancy measures of the MDP. If the MDP has non-unique optimal policies, the dual problem also has non-unique optimal solutions, which implies a degeneration of the primal problem. Degeneration here means that some constraints are redundant at the primal optimal solution (i.e., the corner-point solution is at the intersection of more than ms hyperplanes). Since the rows of the primal LP are linearly independent, we know that in this case, there are multiple (K ′ > 1) choices for the set of basic variables (vBk )1≤k≤K′ at the optimal solution. When the coefficients in the intersecting hyperlanes perturb slightly along a given direction, the objective value will change by a perturbation of the objective coefficients along a chosen set of basic variables vBk . In other words, we can partition the set of directions U into subsets {Uk}1≤k≤K′ such that, if the direction of perturbation of (P, µR), say u, lies in Uk, then the LP optimal value perturbs by fixing the basic variables as vBk . Denote Gρk as the gradient vector corresponding to this direction. If some of the G ρ k’s are equal, we merge the corresponding Uk’s into one partition set. Thus, we have Kρ ≥ 1 distinct Gρk’s and a partition of U = ∪1≤k≤KρUρk , where U ρ k = {u : Dρu(P, µR) = G ρ k}, and Dρu(P, µR) denotes the gradient of∑
s ρ(s)V (s) with respect to (P, µR) along the direction u. Note that the argument so far focuses on the LP with objective value ∑ s ρ(s)V (s). However, we can repeat the same argument for each V (s) by setting ρ(s) = es. For any u ∈ U and s ∈ S, denote the directional gradient of V (s) with respect to P, µR by DuV (P, µR)(s), thus we can define the directional Jacobian of V with respect to P, µR as DuV (P, µR) := [DuV (P, µR)(1), . . . , DuV (P, µR)(s), . . . , DuV (P, µR)(ms)]
T , which leads to K ≥ 1 (potentially larger than Kρ) distinct Gk’s, where Gk ∈ Rms×(Nms+N) and U = ∪1≤k≤KUk is partitioned with Uk = {u : DuV (P, µR) = Gk}. Define ûn = (P̂n−P, µ̂R,n−µR)/ √ ||P̂n − P ||2 + ||µ̂R,n − µR||2. We
have
V̂ ∗n − V ∗ = K∑ k=1 Gk1 (ûn ∈ Uk) (P̂n − P, µ̂R,n − µR) + oP (||(P̂n − P, µ̂R,n − µR)||)
Multiply √ n on both sides and notice that
ûn = ( √ n(P̂n − P ), √ n(µ̂R,n − µR))/ √ || √ n(P̂n − P )||2 + || √ n(µ̂R,n − µR)||2
is a continuous mapping of ( √ n(P̂n − 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 xs,a as the occupancy measure
xs,a = ∞∑ t=0 γtPρ(St = s)
where Pρ denotes the distribution of St’s with initial distribution ρ. Then, xs,a satisfies the LP max ∑ s,a µR(s, a)xs,a
subject to ∑ s,a µC(s, a)xs,a ≤ η∑ a xs,a − γ ∑ s′,a P (s|s′, a)xs′,a = ρ(s), ∀s
xs,a ≥ 0, ∀s, a
(11)
(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 ms + 1, and a corner-point optimal solution has ms + 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, xs,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, xs,a > 0 for only one a, for all s except one state sr, where we can have xsr,a∗1(sr) > 0 and xsr,a∗2(sr) > 0 for two distinct actions a ∗ 1(sr), a ∗ 2(sr). 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(sr)|sr), and so π∗(a∗2(sr)|sr) = 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 ,PV ∗(µ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 ∇xf , which denotes the multi-dimensional Jacobian matrix, from ∂xf , which is used to denote the Jacobian when x is a scalar (1-dimensional).
∇µR,µC ,PV ∗(µR, µC , P ) = ∇µR,µC ,P,α∗V ∗(µR, P, α∗)[I,∇µR,µC ,Pα∗(µC , P )T ]T
= ∇µR,µC ,PV ∗(µR, P, α∗) + ∂α∗V ∗(µR, P, α∗)∇µR,µC ,Pα∗(µC , P ) (12)
Differentiating ρTLπ ∗ (µR, µC , P ) = η,
we have
ρT∇µR,µC ,PLπ ∗ (µ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 ,PLπ ∗ (µC , P, α ∗)
Substituting this into (12), we get
∇µR,µC ,PV ∗(µR, µC , P ) = ∇µR,µC ,PV (µR, P, α∗)
− ρ T∇µR,µC ,PLπ
∗ (µC , P, α ∗)∂α∗V ∗(µR, P, α ∗)
ρT∂α∗Lπ ∗(µC , P, α∗) Next, define an ms-dimensional vector rC by rC(s) = ∑ma j=1 µC(s, j)π ∗(j|s). Then
∂α∗L π∗(µC , P, α ∗) = ∇rC ,Pπ∗L π∗(rC , P π∗)[(∂α∗rC(α ∗))T , (∂α∗P π∗(α∗))T ]T
Note that (I − γPπ∗)Lπ∗ = rC . By applying the implicit function theorem, we have
∂α∗L π∗(µC , P, α ∗) = (I − γPπ ∗ )−1 I,  Lπ∗ . . . Lπ∗
. . . Lπ∗
  q0L = (I − γPπ ∗ )−1qL,
where q0L is a vector with q 0 L(sr) = µR(sr, a ∗ 1(sr)) − µR(sr, a∗2(sr)), q0L(ms + (sr − 1)ma + j) = P (j|sr, a∗1(sr))− P (j|sr, a∗2(sr)) for 1 ≤ j ≤ ma and q0L(i) = 0 for any other index i. Similarly, we can define rR by rR(s) = ∑ma j=1 µR(s, j)π ∗(j|s), and we have
∂α∗V ∗(µR, P, α ∗) = (I − γPπ ∗ )−1qV
The derivation of ∇µR,µC ,PV ∗(µR, P, α∗) and ∇µR,µC ,PLπ ∗ (µ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
[µ̂R,n, P̂n]S0 = [µ̂R,n((i− 1)ma + j), P̂n((i− 1)N + (j − 1)ms + k)]i∈S0,1≤j≤ma,1≤k≤ms ,
and [µR, P ]S0 = [µR((i− 1)ma + j), P ((i− 1)N + (j − 1)ms + k)]i∈S0,1≤j≤ma,1≤k≤ms .
By Assumption 6, we have
√ n([µ̂R,n, P̂n]S0 − [µR, P ]S0)⇒ N (0,ΣR,P,S0) where ΣR,P,S0 =
( W−1S0 D S0 R 0
0 DS0Q
) .
NoticeMg ◦MS0I ◦Tµ̂R,n,P̂n only involves random variables [µ̂R,n, P̂n]S0 . Changing the distribution of [µ̂R,n, P̂n]S\S0 will not change the distribution of Q̂ M n . We can thus assign auxiliary random variables to µ̂R,n and P̂n for all i /∈ S0, 1 ≤ j ≤ ma, 1 ≤ k ≤ ms. In particular, we define independent random variables for each i /∈ S0 by letting
µ̂R,n((i− 1)ma + j) D = 1√ n N (µR((i− 1)ma + j), 1)
P̂n((i− 1)N + (j − 1)ms + k) D = 1√ n N (P ((i− 1)N + (j − 1)ms + k), 1).
By doing so, we extend thems0ma-dimensional random variable [µ̂R,n, P̂n]S0 to anmsma-dimensional random variable [µ̂R,n, P̂n]S and
√ n([µ̂R,n, P̂n]S − [µR, P ]S)⇒ N (0,ΣR,P,S) where ΣR,P,S =
( ΣR,P,S0 0 0 I ) .
Similar to the proof of Theorem 1, define
FM (Q ′, r′, P ′) = Q′ −Mg ◦MS0I ◦ Tr′,P ′(Q ′).
By Assumption 4,Mg is max-norm non-expansion. Then,Mg ◦MS0I is also max-norm non-expansion, implying that ∇(Mg ◦MS0I ) has all its eigenvalues less than or equal to 1. Thus,
∂FM ∂Q′ = ∇M(Tr′,P ′(Q′))(I − γP̃ ′)
is invertible. By Assumption 4, we have FM (QM , µR, P ) = 0. By Assumption 5, there exists a neighborhood ΩM around (QM , µR, P ), such that FM is continuously differentiable on ΩM . Then, applying the implicit function theorem, we have that there exists an open set EM ⊂ ΩM and a continuously differentiable mapping φM on EM , such that φM (µR, P ) = QM and
∇φM (µR, P ) = − [ ∂FM ∂Q′ ]−1 [ ∂FM ∂r′ , ∂FM ∂P ′ ]∣∣∣∣∣ Q′=QM ,r′=µR,P ′=P .
Using the delta method, we have √ n(Q̂Mn −QM )⇒ N (0,∇φM (µR, P )ΣR,P,S(∇φM (µR, P ))T ) = N (0,ΣMS0).
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)ma + l)P (i|s = k, a = l)π(a = j|s = i) = wπ((i− 1)ma + j)
for any i ∈ S, j ∈ A. Summing up across j’s for each i, we have∑ j wπ((i− 1)ma + j)
= ∑ j ∑ k,l wπ((k − 1)ma + l)P (i|s = k, a = l)π(a = j|s = i)
= ∑ k,l wπ((k − 1)ma + l)P (i|s = k, a = l)
On the other hand, for any w in W, πw satisfies∑ k,l w((k − 1)ma + l)P (i|s = k, a = l)πw(a = j|s = i)
= ∑ k,l w((k − 1)ma + l)P (i|s = k, a = l)w((i− 1)ma + j)/ ∑ u w((i− 1)ma + u)
= ∑ u w((i− 1)ma + u)w((i− 1)ma + j)/ ∑ u w((i− 1)ma + u) = w((i− 1)ma + j)
for all i ∈ S. Thus, w is the stationary distribution of transition matrix P̃πw .

Review summaries:
This paper studies the asymptotic properties of action value function $Q$ and value function $V$. Specifically, the authors assume that we can collect $n$ data points $\{(s_t,a_t,r_t(s_t,a_t),s_{t+1})\}_{t=1}^{n}$. Based on the collected data, the authors calculated the sample mean of the unknown reward function $\widehat r_n$ and transition probability $\widehat{P}(\cdot|s,a)$. They further defined an estimator $\widehat Q_n$, which is the fixed point of the empirical Bellman operator $\widehat{\mathcal{T}}_n$ derived from $\widehat r_n$ and $\widehat P_n$. The authors proved under certain assumptions that $\widehat Q_n\rightarrow Q^*$ almost surely. Based on this argument, they also derived a similar convergence of value function $\widehat V_n$ in distribution. Confidence intervals can also be established based on these results.

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.