Batch Policy Learning in Average Reward Markov Decision Processes

Introduction

Mobile health (mHealth) is a rapidly growing field due to the recent advances in mobile and sensing technologies. The mHealth intervention provides a unique opportunity to promote the healthy behaviors (e.g., regular physical activity and adherence to medications) and has been successfully applied in many health fields (e.g., smoking cessation, physical activity, drug abuse and diabetes). Just-in-time adaptive interventions (JITAI, Nahum-Shani et al. 2016) use a decision rule (i.e., a treatment policy or policy) that maps real-time information about the individual’s context to a particular treatment. In this work we study the problem of how to use data consisting of multiple trajectories to estimate a policy that leads to good long-term performance.

We model the sequential decision making process by a time-homogeneous Markov Decision Process (MDP) (Puterman 1994) over infinite time horizon. This framework is natural for mobile health applications in which the number of decision times is often large. For example, in HeartSteps, a physical activity mHealth study, there are five decision times per day, resulting in thousands of decision times over a year. Tremendous progress has been made in finite horizon setting; see the recent review by Kosorok and Laber 2019 for references therein. However when the number of time points is very large, methods that are based on the idea of backward iteration (e.g., Q-learning) or importance sampling (Precup 2000) may suffer a large variance in problems or even be unpractical (Voloshin et al. 2019; Laber et al. 2014).

We propose to estimate the policy that optimizes the long-term average outcomes (rewards) using data consisting of multiple trajectories of finite length. The majority of existing methods focuses on the alternative, the discounted sum of rewards (Sutton and Barto 2018); see the recent works in statistics (Luckett et al. 2019; Ertefaie and Strawderman 2018; Shi et al. 2020; Shi et al. 2021). The discounted formulation weighs immediate rewards more heavily than rewards further in the future, which is practical in some applications (e.g., finance). However, for mHealth applications, choosing an appropriate discount rate could be non-trivial. The rewards (i.e., the health outcomes) in the distant future are as important as the near-term ones, especially when considering maintenance of health behaviors as well as longer term treatment burden. This suggests using a large discount rate. However, it is well known that algorithms developed in the discounted setting can become increasingly unstable as the discount rate goes to one; see for example Naik et al. 2019. The long-term average reward framework provides a good approximation to the long-term performance of a desired treatment policy in mHealth. Indeed, it can be shown that under regularity conditions the finite average of the expected rewards converges sublinearly to the long-term average reward as time goes to infinity (Hernández-Lerma and Lasserre 1999). Therefore, a policy that optimizes the average reward would approximately maximize the sum of the rewards over a sufficiently long time horizon.

In this work, we present a novel algorithm that estimates the optimal policy in a prespecified, parametric policy class. Various methods have been proposed to estimate the global optimal policy by estimating the optimal Q-function;see for example Ormoneit and Sen 2003; Lagoudakis and Parr 2003; Ernst et al. 2005; Munos and Szepesvári 2008; Antos, Szepesvári and Munos 2008a; Antos, Szepesvári and Munos 2008b; Ertefaie and Strawderman 2018; Fujimoto, Meger and Precup 2019; Kumar et al. 2019; Agarwal, Schuurmans and Norouzi 2020. In practice, the optimal Q-function could be highly non-smooth and complex, thus requiring the use of a flexible function class. This usually results in a learned policy that is also complex. If interpretability is important, this is problematic. Furthermore, when the training data is limited, the flexible function class might overfit the data and thus the variance of the estimated value function and the corresponding policy could be high. Restricting to a pre-specified policy class was studied by Zhang et al. 2012; Zhang et al. 2013; Zhou et al. 2017; Zhao et al. 2015; Zhao et al. 2019; Athey and Wager 2017 in finite time horizon problems and by Luckett et al. 2019; Murphy et al. 2016; Liu et al. 2019 in infinite time horizon problems. The restriction to a simple policy class enhance the interpretability of the learned policy and reduces the variance of the learned policy, although this induces a bias when the optimal policy is not in the class (i.e., trading off the bias and variance).

To efficiently learn an optimal policy in a prespecified policy class, the main statistical challenge is to construct an estimator for the average reward of a policy that is both data-efficient and performs uniformly well when optimizing over the policy class. Our first contribution of this work is a novel doubly robust estimator (see Section 3); we show that this estimator achieves the semiparametric efficiency bound under certain conditions on the estimation error of nuisance functions (see Section 5). Doubly robust estimators have been developed in the finite time horizon problems (Robins, Rotnitzky and Zhao 1994; Murphy et al. 2001; Dudík et al. 2014; Jiang and Li 2016; Thomas and Brunskill 2016) and recently in the discounted reward infinite horizon setting (Kallus and Uehara 2019a; Tang et al. 2020). To the best of our knowledge, our doubly robust estimator for the long-term average reward is new. In the literature of the average MDP in the batch setting, only the non-doubly robust estimator proposed by Liao, Klasnja and Murphy 2019 for the long-term average reward can be shown to achieve the semi-parametric efficiency, although they did not explicitly derive it. Most of the previous works on the policy optimization/evaluation under this framework are focused the online setting or under parametric models (Mahadevan 1996; Abounadi, Bertsekas and Borkar 2001; Wan, Naik and Sutton 2021, e.g.,). Theoretical studies on the average reward MDP in the batch setting are very limited, especially under non-parametric models.

To establish the semiparametric efficiency of the doubly robust estimator and the regret bound, we derive finite-sample error bounds for two nuisance function estimators, a relative value estimator and a ratio estimator. The obtained error bounds are shown to hold uniformly over the prespecified class of policies. Both the relative value and ratio estimators are both derived from the same principle (i.e., coupled estimation; see Section 4). In the case of the ratio estimator, we use an iterative procedure to obtain a near-optimal error bound for the ratio estimator. To the best of our knowledge, this is the first theoretical result characterizing the ratio estimation error, which might be of independent interest.

The rest of the article is organized as follows. Section 2 formalizes the decision making problem and introduces the average reward MDP. Section 3 presents the proposed method of learning the in-class optimal policy, including the doubly robust estimator for average reward (Section 3.3). In Section 4, the coupled estimators of the policy-dependent nuisance functions are introduced. Section 5 provides a thorough theoretical analysis on the regret bound of our proposed method. In Section 6, we describe a practical optimization algorithm when Reproducing Kernel Hilbert Spaces (RKHSs) are used to model the nuisance functions. We further conduct several simulation studies to demonstrate the promising performance of our method in Section 7. All the technical proofs are postponed to the supplementary material.

Problem Setup

Suppose we observe a training dataset, $\mathcal{D}_{n}=\{D_{i}\}_{i=1}^{n}$ that consists of $n$ independent, identically distributed (i.i.d.) observations of $D$:

We use $t$ to index the decision time. The length of the trajectory, $T$, is a fixed constant. $S_{t}\in\mathcal{S}$ is the state at time $t$ and $A_{t}\in\mathcal{A}$ is the action (treatment) selected at time $t$. We assume the action space, $\mathcal{A}$, is finite. To eliminate unnecessary technical distractions, we assume that the state space, $\mathcal{S}$, is finite; this assumption imposes no practical limitations and can be extended to the general state space.

Consider a time-stationary, Markovian policy, $\pi$, that takes the state as input and outputs a probability distribution on the action space, $\mathcal{A}$, that is, $\pi(a|s)$ is the probability of selecting action, $a$, at state, $s$. The average reward of the policy, $\pi$, is defined as

The induced Markov chain, $P^{\pi}$, is irreducible for $\pi\in\Pi$.

The goal of this paper is to develop a method that can efficiently use the training data, $\mathcal{D}_{n}$, to learn a policy that maximizes the average reward over $\Pi$. We propose to construct $\hat{\eta}_{n}^{\pi}$, an efficient estimator for the average reward, $\eta^{\pi}$, for each policy $\pi\in\Pi$ and learn an optimal policy by solving

The performance of $\hat{\pi}_{n}$ is measured by its regret, defined as

Note that although the average reward of the learned policy, $\hat{\pi}_{n}$, is defined over an infinite horizon, the goal here is to characterize the regret based on using a finite number of trajectories, $n$, hence the finite sample regret bound is in terms of $n$. Indeed while the average reward, $\eta^{\pi}$ is defined as $t^{*}\rightarrow\infty$ (2.1), the Markovian and stationary assumptions allow us to estimate $\eta^{\pi}$ using fixed length trajectories.

Doubly Robust Estimator for Average Reward

In this section we present a doubly robust estimator for the average reward for a given policy. The estimator is derived from the efficient influence function (EIF). Below we first introduce two functions that occur in the EIF of the average reward. Throughout this section we fix a time-stationary Markovian policy, $\pi$, and focus on the setting where the induced Markov chain, $P^{\pi}$, is irreducible (Assumption 1).

First, we define the relative value function by

The relative value function, $Q^{\pi}$, and the average reward, $\eta^{\pi}$, are closely related via the Bellman equation:

We now introduce the ratio function. For $t=1,\dots,T$, let $d_{t}(s,a)$ be the probability mass of state-action pair at time $t$ in the trajectory $D$ generated by the behavior policy. Denote by $d_{D}(s,a):=(1/T)\sum_{t=1}^{T}d_{t}(s,a)$ the average probability mass across the $T$ decision times in $D$. Similarly, define $d_{t}(s)$ as the marginal distribution of $S_{t}$ and $d_{D}(s)=(1/T)\sum_{t=1}^{T}d_{t}(s)$ as the average distribution of states in the trajectory $D$. Recall that $T$ is the fixed length of the trajectory, $D$; $d_{D}$ describes the distribution of this finite length trajectory. Further recall that under Assumption 1, the stationary distribution of $P^{\pi}$ exists and is denoted by $d^{\pi}(s)$. We assume the following conditions on the data-generating process.

There exists some $p_{\min}>0$, such that $\pi_{b,t}(a|H_{t})\geq p_{\min}$ for all $a\in\mathcal{A}$ and $1\leq t\leq T$ almost surely.

The average distribution $d_{D}(s)>0$ for all $s\in\mathcal{S}$.

Under Assumption 2, it is easy to see that $d_{D}(s,a)\geq p_{\min}\cdot\left(\min_{s}d_{D}(s)\right)>0$ for all state-action pair, $(s,a)$. It essentially states that the data generating process ensures that every state-action pair $(s,a)\in\mathcal{S}\times\mathcal{A}$ has a positive probability of being visited, which is a standard assumption in the literature. See e.g., Theorem 7 of Kallus and Uehara 2019b and (A2) of Shi et al. 2020. In particular, Assumption (2-1) is often satisfied in randomized trials. See our mobile health application in Section 8. In addition, the batch data of our mobile health application consist of $37$ trajectories with $210$ decision points on each trajectory. In this application, as long as every state has a positive probability of being visited in at least one of $210$ decision points, Assumption (2-2) is also satisfied. Assumption (2-2) is imposed on the data generating process. We essentially require that across an infinite number of draws from this data generating process/trajectory, every state $s\in\mathcal{S}$ will be observed. Note that Assumption 2 does not require that the form of the behavior policy is known. Now we can define the ratio function:

The ratio function plays a similar role as the importance weight in finite horizon problems. While the classic importance weight only corrects the distribution of actions between behavior policy and target policy, the ratio here also involves the correction of the states’ distribution. The ratio function is connected with the average reward by

for any fixed trajectory length, $T$. An important property of $\omega^{\pi}$ is that for any state-action function $f(s,a)$ (not only $Q^{\pi}$),

This orthogonality is key to develop the estimator for $\omega^{\pi}$ (see Section 4.3).

2 Efficient influence function

In this subsection, we derive the EIF of $\eta^{\pi}$ for a fixed policy $\pi$ under time-homogeneous Markov Decision Process described in Section 2. Recall that the semiparametric efficiency bound is the supremum of the Cramèr-Rao bounds for all parametric submodels (Newey 1990). EIF is defined as the influence function of a regular estimator that achieves the semiparametric efficiency bound. For more details, refer to Bickel et al. 1993 and Van der Vaart 2000. The EIF of $\eta^{\pi}$ is given by the following theorem. The proof is provided in Appendix A.

Suppose the states in the trajectory, $D$, evolve according to the time-homogeneous Markov process and Assumption 2 holds. Consider a policy, $\pi$, such that Assumption 1 holds. Then the EIF of the average reward, $\eta^{\pi}$, is

Recall that we impose the Markovian and time-stationary assumptions on the data-generating process. Even though the parameter of interest here (i.e. the average reward $\eta^{\pi}$) is one-dimensional, there may exist multiple, non-efficient influence functions as a result of these assumptions on the multivariate distribution.

3 Doubly robust estimator

We have the following doubly robustness of this estimator (the proof is given in Appendix A).

Suppose $\hat{U}_{n}^{\pi}(s,a)$ and $\hat{\omega}_{n}^{\pi}(s,a)$ converge in probability to deterministic limits $\bar{U}^{\pi}(s,a)$ and $\bar{\omega}^{\pi}(s,a)$ uniformly over $\mathcal{S}\times\mathcal{A}$. If either $\bar{U}^{\pi}=U^{\pi}$ or $\bar{\omega}^{\pi}=\omega^{\pi}$, then $\hat{\eta}_{n}^{\pi}$ converges to $\eta^{\pi}$ in probability.

The uniform convergence in probability can be relaxed to $L_{2}$ convergence by using uniform laws of large numbers. The doubly robustness can protect against potential model mis-specifications since we only require one of two models is correct. Moreover, the doubly robust structure can be used to relax the required rate for each of the nuisance function estimation to achieve the semiparametric efficiency bound, especially if we use sample-splitting techniques (see Section 9), as discussed in Chernozhukov et al. 2018.

Estimators for the Nuisance Functions

Recall the doubly robust estimator (3.6) requires the estimation of two nuisance functions, $U^{\pi}$ and $\omega^{\pi}$. It turns out that although these two nuisance functions are defined from different perspectives, both nuisance functions can in fact be characterized in a similar way. Both estimators can be obtained by minimizing an objective function that involves a minimizer of another objective function (“coupled estimation”). This can be viewed as a generalization of the classical M-estimator with a “plug-in estimator” in the sense that the the second objective function also involves the unknown parameters to be estimated. The idea of coupled estimation was previously used by Antos, Szepesvári and Munos 2008a; Farahmand et al. 2016 to estimate the value function in the discounted reward setting and recently by Liao, Klasnja and Murphy 2019 in the average reward setting. In what follows we provide a general coupled estimation framework and discuss the motivation for using it. We then review the coupled estimator for relative value function and ratio function in Liao, Klasnja and Murphy 2019.

Consider a setting where the true parameter (or function), $\theta^{*}$, can be characterized as the minimizer of the following objective function:

2 Relative value function estimator

Let $Z_{t}=(S_{t},A_{t},S_{t+1})$ be the transition sample at time $t$. For a given $(\eta,Q)$ pair, let

where $\hat{g}_{n}^{\pi}(\cdot,\cdot;\eta,Q)$ is the projected Bellman error at $(\eta,Q)$:

Given the estimator of the (shifted) relative value function, $\hat{Q}_{n}^{\pi}$, we form the estimator of $U^{\pi}$ by $\hat{U}_{n}^{\pi}(s,a,s^{\prime})=\sum_{a^{\prime}}\pi(a^{\prime}|s^{\prime})\hat{Q}^{\pi}_{n}(s^{\prime},a^{\prime})-\hat{Q}_{n}^{\pi}(s,a).$

Throughout this paper, we use tuning parameters, $(\lambda_{n},\mu_{n})$, that do not depend on the policy. In the setting where the policy class is highly complex and the corresponding relative value functions are very different, it could be beneficial to select the tuning parameters locally at a cost of higher computation burden.

3 Ratio function estimator

Below we derive the estimator for the ratio function, $\omega^{\pi}$ using the coupled estimation framework. In particular we estimate a scaled version of the ratio function (denoted by $e^{\pi}$ below) and then convert this back to an estimator of $\omega^{\pi}$. To estimate $e^{\pi}$, we first construct a new MDP and estimate the relative value function for this new MDP (denoted by $H^{\pi}$) using the coupled estimation framework. The estimator of $e^{\pi}$ is then derived from the estimator of $H^{\pi}$.

By definition, $\sum_{s,a}e^{\pi}(s,a)d^{\pi}(s)\pi(a|s)=1$. If we were to replace the reward function in our MDP by $1-e^{\pi}(s,a)$, then the “average reward” of $\pi$ in this new MDP is constant and equal to zero under Assumption 1 (i.e., $\sum_{s,a}\left\{1-e^{\pi}(s,a)\right\}d^{\pi}(s)\pi(a|s)=0$). The “relative value function” of policy $\pi$ under the new MDP is,

Note that $H^{\pi}$ is well-defined under Assumption 1. Furthermore, consider the following Bellman equation for the new MDP:

where for any $H\in\mathcal{F}$, $\hat{g}_{n}^{\pi}(\cdot,\cdot;H)$ solves

Recall that $e^{\pi}$ can be written in terms of $H^{\pi}$ by (4.7); that is, re-arranging terms,

The above ratio function estimator was developed by Liao, Klasnja and Murphy 2019. In this paper we, for the first time, derive a finite-sample error bound for this ratio function estimator, uniformly over the policy class (Theorem B.2 in the appendix). This is the key element in establishing the finite-sample bound regret bound for the estimated optimal policy.

Our ratio function estimator is different from most in the existing literature, such as Liu et al. 2018; Uehara and Jiang 2019; Nachum et al. 2019; Zhang et al. 2020, which are obtained by min-max based estimating methods. For example, Liu et al. 2018 aimed to estimate the ratio between stationary distribution induced by a known, Markovian time-stationary behavior policy and target policy, which is then used to estimate the average reward of a given policy. This is not suitable for the setting where the behavior policy is history dependent. Uehara and Jiang 2019 estimated the ratio, $\omega^{\pi}(s,a)$, based on the observation that for every state-action function $f$,

where $\Delta$ is a simplex space and $\cal F^{\prime}$ is a set of discriminator functions. This method minimizes the upper bound of the bias of their average reward estimator if the state-action value function is contained in $\mathcal{F}^{\prime}$. They proved consistency of their ratio and average reward estimators in the parametric setting, that is, where $\omega^{\pi}(S_{t},A_{t})$ can be modelled parametrically and $\mathcal{F}^{\prime}$ is a finite dimensional space. Subsequently Zhang et al. 2020 developed a general min-max based estimator by considering variational $f$-divergence, which subsumes the case in Uehara and Jiang 2019. Unfortunately, there are no error bounds guarantee for ratio function estimators developed in Uehara and Jiang 2019 and Zhang et al. 2020. Our ratio estimator appears closely related to the estimation developed by Nachum et al. 2019 as they also formulated the ratio estimator as a minimizer of a loss function. However, relying on the Fenchel’s duality theorem, they still use the min-max based method to estimate the ratio. Furthermore, their method cannot be applied in average reward settings. Instead of using min-max based estimators, we use coupled estimation. This will facilitate the derivation of estimation error bounds as will be seen below. We will derive the estimation error of the ratio function, which will enable us to provide a strong theoretical guarantee, and finally demonstrate the efficiency of our average reward estimator without imposing restrictive parametric assumptions on the nuisance function estimations, see Section 5 below.

Theoretical Results

In this section, we provide a finite sample bound on the regret of $\hat{\pi}_{n}$ defined in (2.4), i.e., the difference between the optimal average reward in the policy class, $\Pi$, and the average reward of the estimated policy, $\hat{\pi}_{n}$.

We make use of the following assumption on $\Pi$.

There exists $L_{\Theta}>0$, such that for $\theta_{1},\theta_{2}\in\Theta$ and for all $(s,a)\in\mathcal{S}\times\mathcal{A}$, the following holds

There exists constants $C_{0}>0$ and $0\leq\beta<1$, such that for every $\pi\in\Pi$, the following hold for all $t\geq 1$:

The Lipschitz property of the policy class (3-2) is used to control the complexity of nuisance function induced by $\Pi$, that is, $\{U^{\pi}(\cdot,\cdot,\cdot):\pi\in\Pi\}$ and $\{\omega^{\pi}(\cdot,\cdot):\pi\in\Pi\}$. This is commonly assumed in the finite-time horizon problems (e.g., Zhou et al. 2017). Assumptions (3-1) and (3-2) can be easily satisfied by many policy classes such the one we used in Section 7. Our analysis can be extended to more general policy classes if a similar complexity property holds for these two nuisance function classes. Intuitively the constant $\beta$ in the assumption (3-3) relates to the “mixing time” of the Markov chain induced by $\pi\in\Pi$. A similar assumption was used by Van Roy 1998; Liao, Klasnja and Murphy 2019 in average reward setting. Specifically, Equation (5.1) in Assumption (3-3) is used to show that two nuisance functions $U^{\pi}$ and $\omega^{\pi}$ are Lipschitz continuous with respect to the policy parameter $\theta$ so that we can quantify their estimation error uniformly over the policy class. See Lemma C.1 of Supplementary Material for more details. Equation (5.2) in Assumption (3-3) basically requires an exponential convergence rate of the policy induced Markov chain to the stationary distribution in terms of the expectation under the $L_{2}$-norm with respect to the data generating process. This assumption, together with Assumption (5-3) stated below is used to guarantee the Bellman operator for $U^{\pi}$ based on Equation (3.2) (or a similar quantity related to the ratio function estimation defined in Lemma B.4 of Supplementary Material) is well-posed in the sense of $L_{2}$-norm with respect to the data generating process so as to derive their estimation errors. See Lemma B.5 of Liao, Klasnja and Murphy 2019 and Lemma B.4 of Supplementary Material for more details.

Recall that we use the same pair of function classes $(\mathcal{F},\mathcal{G})$ in the coupled estimation for both $U^{\pi}$ and $\omega^{\pi}$. We make the following assumptions on $(\mathcal{F},\mathcal{G})$.

The function classes, $(\mathcal{F},\mathcal{G})$, satisfy the following:

$\mathcal{F}\subset\mathcal{B}(\mathcal{S}\times\mathcal{A},F_{\max})$ and $\mathcal{G}\subset\mathcal{B}(\mathcal{S}\times\mathcal{A},G_{\max})$

The regularization functionals, $J_{1}$ and $J_{2}$, are pseudo norms and induced by the inner products $J_{1}(\cdot,\cdot)$ and $J_{2}(\cdot,\cdot)$, respectively.

Let $\mathcal{F}_{M}=\{f\in\mathcal{F}:J_{1}(f)\leq M\}$ and $\mathcal{G}_{M}=\{g\in\mathcal{G}:J_{2}(g)\leq M\}$. There exists $C_{1}$ and $\alpha\in(0,1)$ such that for any $\epsilon,M>0$,

We now introduce the assumption that is used to bound the estimation error of value function uniformly over the policy class. Define the projected Bellman error operator:

The triplet, $(\Pi,\mathcal{F},\mathcal{G})$, satisfies the following:

A similar set of conditions are employed to bound the estimation of ratio function. For $\pi\in\Pi$ and $H\in\mathcal{F}$, define the projected error:

where, as before, $\Delta^{\pi}(Z_{t};H)=1-H(S_{t},A_{t})+\sum_{a^{\prime}}\pi(a^{\prime}|S_{t+1})H(S_{t+1},a^{\prime})$.

The triplet, $(\Pi,\mathcal{F},\mathcal{G})$, satisfies the following:

$e^{\pi}(\cdot,\cdot)\in\mathcal{G}$, for $\pi\in\Pi$.

There exists two constants $C_{2}^{\prime},C_{3}^{\prime}$ such that $J_{2}\left\{g^{*}_{\pi}(\cdot,\cdot;H)\right\}\leq C_{2}^{\prime}+C_{3}^{\prime}J_{1}(H)$ holds for $H\in\mathcal{F}$ and $\pi\in\Pi$.

Suppose Assumptions 1 to 6 hold. Let $\hat{\pi}_{n}$ be the estimated policy (2.3) in which the nuisance functions are estimated with tuning parameters $\mu_{n}=\lambda_{n}=\mu_{n}^{\prime}=\lambda_{n}^{\prime}=Ln^{-1/(1+\alpha)}$, for some constant $L>0$. Define $\beta_{k}=\frac{1}{1+\alpha}\left\{1-(1-\alpha)2^{-k+1}\right\}$. Fix any integer $k\geq 2$, $\delta\in(0,1)$ and sufficiently large $n$. With probability at least $1-\delta$, we have

Recall that $p$ is the number of parameters in the policy, $\alpha$ is given in (4-4), and $n$ is the number of trajectories in the data. Theorem 5.1 shows that when the tuning parameters are of the order $O(n^{-1/(1+\alpha)})$, the regret of the estimated policy is $O(p^{1/2}n^{-1/2}+pn^{-\beta_{k}})$. The leading term (in terms of $n$), $O(\sqrt{p/n})$, corresponds to the regret of an estimated policy as if the nuisance functions are known beforehand. The second term is due to the estimation error of nuisance functions. In particular, we show in Theorem B.1 in Section B of the appendix that the uniform estimation error of the relative value function is of $O(pn^{-1/(1+\alpha)})$ and in Theorem B.2 in the same section that the uniform estimation error of ratio is of $O(pn^{-\beta_{k}})$ (see the remark after Theorem B.2 for why the rate depends on $k$). Note that the error of ratio is the dominant term as $\beta_{k}<1/(1+\alpha)$ and $\beta_{k}$ can be chosen arbitrarily close to $\frac{1}{1+\alpha}$ by choosing a sufficiently large $k$. Therefore the proposed ratio estimator can achieve the near-optimal nonparametric convergence rate. See the proof of Theorem B.2 in Section B.1 of the appendix for more details. To the best of our knowledge, this is the first result that characterizes the regret of the estimated optimal in-class policy in the infinite horizon setting.

2 Asymptotic results

In this section, we prove that the average reward for our estimator of the optimal policy converges to the optimal average reward at a parametric rate (i.e., $\sqrt{n}$). Recall $\phi^{\pi}(D)$ is the efficient influence function of $\eta^{\pi}$ given in Theorem 3.1.

Suppose Assumptions 1 to 6 hold. For each $n\geq 1$, let $\hat{\eta}_{n}^{\pi}$ be the doubly robust estimator defined in (3.6) and $\hat{\pi}_{n}$ be the estimated policy defined in (2.3) with tuning parameters $\mu_{n}=\lambda_{n}=\mu_{n}^{\prime}=\lambda_{n}^{\prime}=Ln^{-1/(1+\alpha)}$, for some constant $L>0$. Then as $n\rightarrow\infty$,

The first result shows that the estimated average reward by the doubly robust estimator reaches the semiparametric efficiency bound when we plug in the estimator for the two nuisance functions. The double robustness structure ensures that the estimation error of nuisance functions is only of lower order and does not impact the asymptotic variance of the estimated average reward. The second result shows the asymptotic of the estimated optimal value, $\hat{\eta}_{n}^{\hat{\pi}_{n}}$, converges to the maximum of the Gaussian process at the optimal policies. When there is a unique optimal policy $\pi^{*}=\operatorname{argmax}_{\pi\in\Pi}\eta^{\pi}$, we have $\sqrt{n}(\hat{\eta}_{n}^{\hat{\pi}_{n}}-\eta^{\pi^{*}})$ weakly converges to a Gaussian distribution. Estimating the limiting distribution could be challenging (especially when there exists non-unique policies) and is left for future work. Alternatively one can consider resampling-based method to construct confidence interval for $\sup_{\pi}\eta^{\pi}$ (see the recent work by Wu and Wang 2020 in single-stage problem).

Practical Implementation

In this section, we describe an algorithm to estimate an in-class optimal policy based on our efficient average reward estimator $\hat{\eta}_{n}^{\pi}$. Without loss of generality, we consider a binary-action setting, i.e., $\mathcal{A}=\{0,1\}$, and the following stochastic parametrized policy class $\Pi$ indexed by $\theta$:

for some pre-specified constant $c>0$. Note that other link functions such as the probit function might be used here instead. Here $\|\cdot\|_{\infty}$ refers to sup-norm in Euclidean space. We fix $c=10$ throughout our paper. In addition, we set $\mathcal{F}$ and $\mathcal{G}$ in the estimation of both value and ratio functions to be Reproducing Kernel Hilbert Spaces (RKHSs) associated with Gaussian kernels because of the representer theorem and the property of universal consistency.

The constraint on $\theta$, $\|\theta\|_{\infty}\leq c$, is used to maintain sufficient stochasticity in our learned policy. The stochasticity facilitates the use of $\hat{\pi}_{n}$ as a “warm start" policy for use by an online algorithm with future individuals. A nice side effect is that the restriction on $\theta$ provides a computational stability and can avoid degenerative cases in policy optimization similar to that when using logistic regression in classification problems (Friedman, Hastie and Tibshirani 2001). As discussed in the introduction, we consider the simple policy class $\Pi$ instead of nonparametric models such as neural networks or tree-based models mainly due to the concern of overfitting. In the batch setting, data are limited and often noisy. Using flexible function classes for modeling the policy may lead to overfitting and thus the variance of the resulting policy could be very large. The use of a simple policy class can reduce the variance while it may induce some possible bias. In addition, interpretability is critical in our batch policy learning problem. The interpretability of decision tree models are often not very stable, whereas neural networks are not very interpretable. Therefore we prefer using this simple policy class $\Pi$.

To obtain $\hat{\pi}_{n}\in\Pi$, we solve a multi-level optimization problem (6.1)-(6.5). Recall a multi-level optimization problem (Richardson 1995) is a optimization problems in which the feasible set is implicitly determined by a sequence of nested optimization problems. It typically consists of an upper level optimization task that represents the objective function, and a series of (possibly nested) lower level optimization tasks that represents the feasible set.

As a reminder, recall that in Section 4 we have defined

Also, the ratio estimator $\omega_{n}^{\pi}$ can be obtained from $\hat{H}^{\pi}_{n}(\cdot,\cdot)$ by using (4.10).

The upper optimization task (6.1) is used to search for $\hat{\pi}_{n}$ and the two parallel lower optimization tasks (6.2)-(6.3) and (6.4)-(6.5) are used to compute two nuisance function estimators for a given $\pi\in\Pi$, i.e., the feasible set, respectively. Note that each nuisance function estimation is itself a nested optimization sub-problem. Multi-level optimization problems in general cannot be computed by iteratively updating solutions to lower problems (6.2)-(6.3) and (6.4)-(6.5), and solutions to the upper problem (6.1), in a similar manner to coordinate descent. Hence, in order to solve this problem, one common approach is to replace the inner optimization problems (6.2)-(6.3) and (6.4)-(6.5) by their corresponding Karush-Kuhn-Tucker (KKT) conditions so that the overall problem can be equivalently formulated as a nonlinear constraint optimization problem. However, this approach can be computationally expensive and may not be suitable for large scale settings. Instead we overcome this computational obstacle by using the representer theorem and obtain the closed-form solutions for our inner optimization problems (6.2)-(6.3) and (6.4)-(6.5) respectively. After plugging these closed-form solutions into (6.1), we can use a gradient-based method to find $\hat{\pi}_{n}$.

In the following subsection, we briefly discuss how to simplify our multi-level optimization problem (6.1) using the representer theorem. The details of computation can be found in Appendix E. For the ease of illustration, we rewrite the training data $\mathcal{D}_{n}$ into tuples $Z_{h}=\{S_{h},A_{h},R_{h},S_{h}^{\prime}\}$ where $h=1,\dots,N=nT$ indexes the tuple of the transition sample in the training set $\mathcal{D}_{n}$, $S_{h}$ and $S_{h}^{\prime}$ are the current and next states and $R_{h}$ is the associated reward. Let $W_{h}=(S_{h},A_{h})$ be the state-action pair, and $W_{h}^{\prime}=(S_{h},A_{h},S_{h}^{\prime})$. Suppose the kernel function for the state is denoted by $k_{0}(s_{1},s_{2})$, where $s_{1},s_{2}\in\mathcal{S}$. In order to incorporate the action space, we can define $k((s_{1},a_{1}),(s_{2},a_{2}))=\mathds{1}_{\{a_{1}=a_{2}\}}k_{0}(s_{1},s_{2})$. Basically, we model each $Q(\cdot,a)$ separately for each arm in the RKHS with the same kernel $k_{0}$. Recall that we have to restrict the function space $\mathcal{F}$ such that $Q(s^{*},a^{*})=0$ for all $Q\in\mathcal{F}$ so as to avoid the identification issue. Thus for any given kernel function $k$ defined on $\mathcal{S}\times\mathcal{A}$, we make the following transformation by defining $k(W_{h},W_{j})=k_{0}(W_{h},W_{j})-k_{0}((s^{*},a^{*}),W_{h})-k_{0}((s^{*},a^{*}),W_{j})+k_{0}((s^{*},a^{*}),(s^{*},a^{*}))$ for any $1\leq h,j\leq N$. One can check that the induced RKHS by $k(\cdot,\cdot)$ satisfies the constraint in $\mathcal{F}$ automatically.

We denote kernel functions for $\mathcal{F}$ and $\mathcal{G}$ by $k(\cdot,\cdot),l(\cdot,\cdot)$ respectively. The corresponding inner products are defined as $\langle\cdot,\cdot\rangle_{\mathcal{F}}$ and $\langle\cdot,\cdot\rangle_{\mathcal{G}}$. We first discuss the inner minimization problem (6.2)-(6.3). Note that this is indeed a nested kernel ridge regression problem, different from the standard ridge regression. The closed form solution can be obtained as $\hat{g}_{n}^{\pi}(\cdot,\cdot;\eta,Q)=\sum_{h=1}^{N}l(W_{h},\cdot)\hat{\gamma}(\eta,Q)$. In particular, $\hat{\gamma}(\eta,Q)=(L+\mu I_{N})^{-1}\delta_{N}^{{\pi}}(\eta,Q)$, where $Q\in\mathcal{F}$ and $L$ is the kernel matrix induced by $l$, $\mu=\mu_{n}N$, and $\delta^{\pi}_{N}(\eta,Q)=(\delta^{\pi}(Z_{h};\eta,Q))_{h=1}^{N}$ is a vector of TD error. Each TD error can be further written as $\delta^{\pi}(Z_{h};\eta,Q)=R-\eta-\langle Q,f_{W^{\prime}_{h}}\rangle_{\mathcal{F}}$ where

Summarizing together and plugging all the intermediate results into (6.1), the multi-level optimization problem can be simplified as:

where $1_{N}$ is a length-$N$ vector of all ones.

2 Optimization

Note that problem (6.6) becomes a smooth nonlinear optimization with box constraints. We use limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with box constraints (L-BFGS-B) to compute the solution $\hat{\theta}$ (Liu and Nocedal 1989). The gradient computing is provided in appendix. The computational complexity/operations of our algorithm is of order $\Upsilon N^{3}p$, where $\Upsilon$ is the number of iterations in our optimization algorithm. The memory requirement is of order $N^{2}p$. One may implement some sub-sampling methods such as stochastic gradient decent to further improve both computation and memory complexity of our algorithm. We will leave it for future work. Although the overall optimization problem is non-convex and, thus an optimal solution may not be achievable, the performance of our numerical experiments in the following section are quite stable and promising. Recently, there is a growing interest in studying statistical properties of algorithm-type of nonconvex M-estimators, e.g., (Mei et al. 2018; Loh et al. 2017). For many practical applications, gradient decent methods with a random initialization have been demonstrated to converge to local minima (or even global minima) that are statistically good. While this is not the focus of our paper, it will be interesting to pursue toward this direction for future research such as studying the landscape of $\eta^{\pi}$ and its related properties.

3 Tuning parameters selection

In this subsection, we discuss the choice of tuning parameters in our method. The bandwidths in the Gaussian kernels are selected using median heuristic, e.g., median of pairwise distance (Fukumizu et al. 2009). The tuning parameters $(\lambda_{n},\mu_{n})$ and $(\lambda^{\prime}_{n},\mu^{\prime}_{n})$ are selected based on 3-fold cross-validation. Given assumptions in Theorems B.1 and B.2 of the appendix that these tuning parameters are independent of the policy $\pi$, we can select them for the ratio and value functions separately. Specifically, for the tuning parameters $(\lambda_{n},\mu_{n})$ in the estimation of value function, we focus on (6.2)-(6.3). For the tuning parameters $(\lambda^{\prime}_{n},\mu^{\prime}_{n})$ in the estimation of ratio function, we focus on (6.4)-(6.5). At the first glance, one may think the selection of tuning parameters will be the same as those in the standard supervised learning. However, this actually requires an additional step as we cannot observe responses when estimating these two coupled estimators (recalled that we need to first compute projected bellman errors), in contrast to the standard kernel regression setting. In the following, we discuss our selection procedure of $(\lambda_{n},\mu_{n})$ and $(\lambda^{\prime}_{n},\mu^{\prime}_{n})$ with more details.

We first randomly choose a set of candidate policies used to gauge our tuning parameters. For each candidate policy, $\pi$, in this set, we can firstly estimate $(\hat{\eta}_{n}^{\pi},\hat{\alpha}(\pi))$ by the proposed method using two folds of data. Then for the value function estimation, we calculate temporal difference errors $\delta^{\pi}(\cdot;\hat{\eta}_{n}^{\pi},\hat{\alpha}(\pi))$ for each transition sample in the validation set. Since we cannot observe/calculate the true bellman error, following the idea in (Farahmand and Szepesvári 2011), we estimate the Bellman error by projecting these temporal differences on the space of $\mathcal{S}\times A$ in the validation set using the standard Gaussian kernel regression. Thus for each policy $\pi$ and each pair of tuning parameters, we output the squared estimated Bellman error in the validation set as a criterion to evaluate the performance of our value function estimation. Since tuning parameters are assumed independent of policies, we then select the tuning parameters that minimize the worst case of estimated Bellman errors among the set of all candidate policies. We use the same strategy to select the tuning parameters for our ratio estimation. The details are given in the Algorithm 1. Without the independent assumptions of tuning parameters from the policies in $\Pi$, one may alternatively choose these tuning parameters jointly by maximizing $\hat{\eta}_{n}^{\pi}$ on the validation set, which requires large computational costs and we omit here. But it would be very interesting to study the theoretical properties of these two cross-validation procedures, or more generally, the selection of tuning parameters in the framework of couple estimation, which we leave it as future work.

Simulation Studies

In this section, we consider two scenarios to evaluate the proposed algorithm. For both scenarios, we consider $S_{t}=(S_{t,1},S_{t,2},S_{t,3})$ as a three-dimensional state at each decision point $t$, and the action space is binary, i.e., $\mathcal{A}=\{0,1\}$. The behavior policy used to generate actions follows Bernoulli distribution with equal probabilities. In addition, the initial state $S_{1}$ is sampled from standard multi-variate normal distribution, i.e., $S_{1}\sim MVN(0,I_{4})$

The first scenario we consider is a standard MDP setting. Let $\xi_{t}$ follows a standard multi-variate normal distribution. Then we generate the transition of states and reward functions via following models:

for $t=1,\cdots,T$. Here $S_{t,3}$ can be interpreted as the treatment burden or fatigue.

The second scenario we consider is a non-stationary environment. In particular, we consider the same transition models as above, but let the reward function to be time-dependent. More specifically, we consider

where the time-varying parameters $\beta_{t}=0.25\times\exp(-0.05(t-1))$ and $\tau_{t}=0.4\times\exp(-0.05(t-1))$. This generative model represents the scenario in which as the study progresses, the overall impact of intervention is decreasing. Note that since the reward function is non-stationary, we do not have a guarantee for our proposed algorithm to find an optimal policy.

We compare with four baseline methods, which were proposed in the setting of the discounted sum of rewards. The first two are recently proposed deep off-policy RL algorithms (Fujimoto, Meger and Precup 2019; Kumar et al. 2019) denoted by BCQ and BEAR respectively. The underlying idea behind these two state-of-art algorithms is to conservatively estimate the optimal $Q$-function on the less explored state-action pair and restrict the resulting policy close to the behavior one. The third method is the celebrated fitted-Q iteration (FQI) method proposed by Ernst et al. 2005. At each iteration, relying on the optimal Bellman equation, FQI algorithm updates the estimation of the optimal $Q$-function via solving a supervised learning problem. The last method is V-learning proposed by (Luckett et al. 2019), which also aims to learn an optimal in-class policy. Since our goal is to maximize the long-term average reward, we set the discount factor $\gamma$ in these four methods as $0.99$ to approximate the average reward for comparison. In addition, to draw a relatively fair comparison, we implement these four methods using the same policy class as ours. Specifically, the first three methods will output an estimation of the optimal $Q$-function (defined in the discounted setting), after which we implement a weighted logistic regression to estimate the optimal in-class policy. For V-learning, we keep the default setup and use the same policy class as ours. Finally, for BEAR and BCQ, we use two-hidden layers neural networks with 32 nodes for each and ReLU activation functions to model the optimal $Q$-function. The other hyper-parameters are either tuned for their best performance, or recommended in the official implementation as robust choices. For FQI, we implement a kernel ridge regression at each iteration with tuning parameters selected similar to our nuisance parameter estimation.

Results of the above two scenarios can be found in Table 1. As we can see, our algorithm performs well in finding optimal in-class stationary policies, compared with the other four baseline methods. Compared with the oracle one with the best average reward about $10$, the regret of our algorithm is almost the smallest among all these methods, which is expected as we aim to maximize the average reward while the other four methods are for maximizing the discounted sum of rewards. For BEAR and BCQ with neural network models, due to the relatively small sample size and a large discount factor $\gamma$, the performances seem unstable. FQI and V-learning methods overall show competitive performances. But it can be seen that V-learning may suffer some large variance. In addition, one possible reason for the high quality performance of our method in Scenario 2 is that the time-dependent effect in reward function is exponentially decaying by our design. Therefore we expect the performance of our algorithm may not be affected severely by the non-stationarity. In addition, it can be seen that as the sample size $n$ or the length of each trajectory $T$ increases, the average rewards of our estimated policies are also improved, demonstrating the appealing performance of the proposed method. Finally, we remark that the maximum running time of our method for one replication in our simulation studies is less than 40 minutes.

Application to mobile health

We apply the proposed method to HeartSteps. HeartSteps is mobile health application focusing on physical activity. Three studies were conducted to develop the intervention. In this work, we apply the proposed method to the data collected from the first study, which we will refer to as HS1 in the throughout. HS1 is a 42-day micro-randomized trial (Klasnja et al. 2015; Liao et al. 2016). Each participant was provided with a Jawbone wrist tracker to collect step count data and specified five decision times, roughly 2.5 hours apart during each day, that would be good times to potentially receive contextually tailored activity suggestion message. In HS1, the activity message was sent with a fixed probability 0.6 at each of the five decision times. Our goal is to use HS1 data to learn a treatment policy that determines at each decision time whether to send the activity message (i.e., binary action).

We construct the state variable using the previous step count (the 30-min step count prior to the decision time and from yesterday), location, temperature and past notifications. We set the reward to be the log transformation of the step count in 30-min window after each decision time. In this analysis, we include 37 participants’ data and exclude the decision times when participants were traveling abroad or experiencing technical issues or when the reward (i.e., post 30-min step count) was considered as missing (Klasnja et al. 2018).

Next, we construct the policy class. In this analysis, we include two state variables in the policy. The first variable is the location (home/work vs. other locations). Location is important because people, in a more structured environment (i.e. at home or work), may respond better to an activity suggestion as compared to when they are at other locations. As a proxy for participant burden, the second variable included in the policy is “dosage”, a discounted sum of the number of past activity messages sent with the discount rate chosen as 0.95. The rationale for using this variable is that receiving too many notifications in the recent past is likely to decrease the effectiveness of sending the activity message due to over-burdening participants. We consider the policy class of the form $\pi_{\theta}(1|s)=\operatorname{expit}(\theta^{\top}\phi(s)),\theta\in\Theta$, where the feature vector $\phi(s)=(1,\text{dosage},\text{location})$ and $\Theta$ is the box constraint within -10 and 10. Here dosage is standardized to be within 0 and 1.

We apply the proposed method with the tuning parameter selected by cross-validation in Algorithm 1. The estimated coefficients are $[10,-10,-4.788]$. Figure 1 shows the estimated policy at different combination of dosage and location. As one would expect, the learned policy tends to send fewer suggestions if the participant received many suggestions in the recent past. Also, the policy indicates that it is more effective to send the message when the user is at home/work location. The estimated average reward of this policy is 3.301. As a comparison, the estimated average reward of the simple location-based policy (i.e., send only when the user is at home/work) is 3.15 and the send-nothing policy is 2.96. Transforming to the scale of the raw step count as that in (Klasnja et al. 2018), the learned policy can result in $16\%$ (i.e., $\exp(3.301-3.15)-1=0.16$) improvement, which is equivalent to $40$ more steps (the mean step count across all decision times in the data is $248$) compared with the simple location-based policy, and $40\%$ (i.e., $\exp(3.301-2.96)-1=0.40$) improvement, or equivalently $101$ steps more, compared with the send-nothing policy. Lastly, we remark that the running time of our real data analysis is about $2$ hours, which is acceptable in the batch setting. This is ultimately different from online RL domains where the policy is usually updated upon the arrival of each observation.

Discussion

Double/Debiased machine learning An alternative way to construct the estimator for the average reward is based on the idea of double/debiased machine learning (a.k.a. cross-fitting, Bickel et al. 1993 and Chernozhukov et al. 2018). There is growing interest in using double machine learning in causal inference and in the policy learning literature (Zhao et al. 2019) in order to relax assumptions on the convergence rates of nuisance parameters. The basic idea is to split the data into $K$ folds. For each of the $K$ folds, construct the estimating equation by plugging in the estimated nuisance functions that are obtained using the remaining $(K-1)$ folds. The final estimator is obtained by solving the aggregated estimation equations. While cross-fitting requires weaker conditions on the nuisance function estimations, it indeed incurs additional computational cost, especially in our setting where nuisance functions are policy-dependent and we aim to search for the in-class optimal policy. Further, this sample splitting procedure may not be stable when the sample size is relatively small, e.g., in a typical mHealth clinical trial. A more efficient way of data splitting under the framework of MDP is needed, which we leave as future work

Computation and optimization Our current algorithm requires relatively large computation and memory because of the non-parametric estimation and the policy-dependent structure of nuisance functions. It is therefore desirable to develop a more efficient algorithm. One possible remedy is to consider a zero-order optimization method such as Bayesian optimization (Snoek, Larochelle and Adams 2012), which is suitable when the dimension of state variables is small. Another possible way to improve the computational efficiency is to first apply some simple algorithm to estimate a sub-optimal policy, based on which we can implement our method to estimate two nuisance parameters. Then one can develop the performance difference lemma in terms of the average reward MDP, similar to that in the discounted setting (Kakade and Langford 2002), to construct a lower bound for $\mathcal{V}(\pi)$ using two estimated nuisance parameters. The last step is to optimize this lower bound for obtaining a better policy. This method may require less computational cost.

Tuning parameters/Model selection In our proposed algorithm, we assume tuning parameters are independent of policies, based on which we develop a min-max cross-validation procedure for the selection of tuning parameters. Model selection in the offline RL setting, which is necessary for improving generalization of RL techniques, is often considered as a challenging task as there is no ground truth available for performance demonstration, in contrast to the online setting with simulated environment. Therefore, it will be interesting to systematically investigate how to perform model selection in offline RL and to provide theoretical guarantees.

Acknowledgements

Peng Liao was supported by NIH grants P50DA039838, R01AA023187, and U01 CA229437. Susan Murphy was supported by NIH grants P50DA039838, R01AA023187, P50DA054039, P41EB028242, U01 CA229437, UG3DE028723, and UH3DE028723. The authors would also like to thank two reviewers, the Associate Editor and the Editor for helpful comments and suggestions that led to substantial improvement in the presentation.

A Semi-parametric efficiency bound and doubly robustness

In this section, we calculate the semi-parametric efficient influence function and proves the doubly robustness. Denote by $L(D;\zeta)$ the likelihood of a parametric sub-model of the data collected over $T$ decision times:

Then the score function of the above parametric sub-model is given by

where $S_{\zeta}(S_{1})$ is $\frac{d\log p_{\zeta}(S_{1})}{d\zeta}$, $S_{\zeta}(S_{t+1}|S_{t},A_{t})=\frac{d\log S_{\zeta}(S_{t+1}|S_{t},A_{t})}{d\zeta}$, and $S_{\zeta}(A_{t}|H_{t})=\frac{d\log S_{\zeta}(A_{t}|H_{t})}{d\zeta}$ for $t=1,\cdots,T$.

If $\phi_{\text{eff}}$ is an influence function, then for any parametric submodel that contains the true parameter $\zeta_{0}$,

Plugging into RHS and using the definition of score function gives

where $(S,A,S^{\prime})\sim d^{\pi}(s,a)P_{\zeta_{0}}(s^{\prime}|s,a)$ follows the stationary distribution under the true model.

For the LHS, we start with the Bellman equation: for any $\zeta$, we have

where we write $U^{\pi}$ and $\eta^{\pi}$ as $U^{\pi}_{\zeta}$ and $\eta_{\zeta}(\pi)$ to explicitly indicate its dependency on $\zeta$. Taking the derivative implies

where in the second last line we use $\int\frac{d}{d\zeta}\mathcal{P}_{\zeta}(s^{\prime}|s,a)d\mu(s^{\prime})=0$. Now averaging over the stationary distribution $d^{\pi}_{\zeta}(s,a)=d_{\zeta}^{\pi}(s)\pi(a|s)$ of the state-action pair gives

For the first term, using the definition of stationary distribution we have

(i) The tangent space $\mathscr{T}$ is given by

(ii) The orthogonal complement of the tangent space $\mathscr{T}$ is

Given the expression of score function $S_{\zeta}$, we can obtain statement (i). In particular, $\mathcal{F}_{1}$ is induced by $S_{\zeta_{0}}(S_{1})$, $\mathcal{F}_{t}$ is induced by $S_{\zeta_{0}}(S_{t+1}|S_{t},A_{t})$ for $t=2,\cdots,T+1$, and $\mathcal{G}_{t}$ is induced by $S_{\zeta_{0}}(A_{t}|H_{t})$ for $t=1,\cdots,T$. See Theorem 1 of (Kallus and Uehara 2020). For any $2\leq t_{1},t_{2}\leq T+1$ and $t_{1}\neq t_{2}$, we can show that $\mathcal{F}_{t_{1}}$ is orthogonal to $\mathcal{F}_{t_{2}}$. Without loss of generality, suppose $t_{1}<t_{2}$. Then for any $q_{t_{1}}\in\mathcal{F}_{t_{1}}$ and $q_{t_{2}}\in\mathcal{F}_{t_{2}}$,

where the second equality is by Markov property. By the similar argument, we can also show that $\mathcal{G}_{t_{1}}$ is orthogonal to $\mathcal{G}_{t_{2}}$ for $1\leq t_{1},t_{2},\leq T$ and $t_{1}\neq t_{2}$. In addition, for $2\leq t\leq T$, we can show $\mathcal{G}_{t}$ is orthogonal to $\mathcal{F}_{t}$ by again similar argument.

In order to derive the orthogonal complement of tangent space $\mathscr{T}$, we first note that

is the space of all random functions with mean zero and finite variance, where

for $t=2,\cdots,(T+1)$, are orthogonal to each other. Then it is enough to project each elements in $\mathcal{F}^{\prime\prime}_{t}$ onto the orthogonal complement of $\left(\mathcal{F}_{t}\bigoplus\mathcal{G}_{t}\right)$ for $2\leq t\leq T$ and $\mathcal{F}_{T+1}$ respectively.

First of all, we can see that $\mathcal{F}^{\prime\prime}_{t}$ is orthogonal to $\mathcal{G}_{t}$ for $2\leq t\leq T$ by the definition of $\mathcal{G}_{t}$. Secondly, it is straightforward to show that $\mathcal{F}^{\prime\prime}_{t}$ is orthogonal to

which is indeed also orthogonal to $\mathcal{F}_{t}$ for $2\leq t\leq(T+1)$. Then projecting each element in $\mathcal{F}^{\prime\prime}_{t}$ onto the orthogonal complement of $\mathcal{F}_{t}$ is equivalent to projecting onto the orthogonal space of

which gives us exactly $\mathcal{F}^{\prime}_{t}$. This concludes statement (ii).

The first term is zero since $\omega^{\pi}(s_{k-1},a_{k-1})\delta^{\pi}(s_{k-1},a_{k-1},s_{k})\in\mathcal{F}_{k}$ by definition of TD error and thus orthogonal to $f(h_{k})\in\mathcal{F}_{k}^{\perp}$. For the second term, for any $1\leq t\leq k-2$, we have $(S_{t},A_{t},S_{t+1})\in\sigma(H_{k-1})$

If $\bar{w}^{\pi}=w^{\pi}$, then by the definition of stationary distribution, we have

If $Q^{\pi}=\bar{Q}^{\pi}$, which implies $U^{\pi}=\bar{U}^{\pi}$, then

converges to $0$ in probability. This concludes our statement. ∎

B Theoretical Results on Nuisance Function Estimation

In this section, we present two finite sample upper bounds for the estimation error of the nuisance functions that holds uniformly over the policy class. These results are needed to prove Theorem 5.1 and 5.2.

We first present the uniform bound for the estimation error of $U^{\pi}$ over $\pi\in\Pi$. This is a generalization of Theorem 1 in Liao, Klasnja and Murphy 2019 in which they focused only on a single policy.

where $\iota=4(\kappa)^{-2}(1+p_{\min}^{-1}(1+(1/T)\|\frac{d_{T+1}}{d_{D}}\|_{\infty}))(1+C_{0}\beta/(1-\beta))^{2}$

Since the proof is similar to that in (Liao, Klasnja and Murphy 2019) with additional efforts on controlling complexity of the policy class $\Pi$, we omit here. Next we present the uniform finite sample bound for the ratio estimator.

where $\omega_{k}=1-2^{-k+1}$ and $\iota^{\prime}=4(\kappa^{\prime})^{-2}(1+p_{\min}^{-1}(1+(1/T)\|\frac{d_{T+1}}{d_{D}}\|_{\infty}))(1+C_{0}\beta/(1-\beta))^{2}$.

Recall that the optimal convergence rate for the classical nonparametric regression problem under the entropy condition (6-4) is $n^{-1/(1+\alpha)}$. Theorem B.2 shows that the the ratio estimator achieves the near-optimal convergence rate. As we have seen in Theorem 5.2, the achieved error rate is enough to guarantee the asymptotic efficiency of the doubly robust estimator; in fact we only need to ensure the error decays faster than $n^{-1/4}$.

Suppose $\mu_{n}=\lambda_{n}=Ln^{-1/(1+\alpha)}$. Fix some $k\in\mathds{N}^{+}$. Under Assumptions (2-1), (2-2), (3-3), (4-1), (4-2), (6-1), (6-2), (6-3), (6-4) and (4-4), the followings hold with probability $1-(3+k)\delta$: for all $\pi\in\Pi$:

Denote the leading constant in Lemma B.1 by ${K_{0}}$. For the choice of tuning parameters $(\lambda_{n},\mu_{n})$, Lemma B.1 and Assumption (6-4) imply that w.p. $1-\delta$, for all $\pi\in\Pi$, the first term in (B.1) can be bounded by

where ${K_{1}}={K_{0}}L(1+2C_{2}^{2})$ and $C_{1}(\delta)={K_{0}}\left(L(1+2C_{1}^{2})+L^{-\alpha}+1+\log(1/\delta)\right)$.

Now consider the second term. Denote the leading constant in Lemma B.2 by ${K_{2}}$. Applying the Decomposition Lemma B.2 implies that w.p. $1-2\delta$, for all $\pi\in\Pi$

With the choice of $(\lambda_{n},\mu_{n})$, $\gamma_{2}(\delta,n,p,\mu_{n},\lambda_{n})$ can be bounded by

As a result, we obtain that with probability at least $1-3\delta$, for all $\pi\in\Pi$,

Initial Rate We derive an initial rate by bounding $\operatorname{Rem}(\pi)$ uniformly over $\pi\in\Pi$. Let

where the function class $\mathcal{F}_{0}$ is given by

and thus we have $\mathcal{F}_{0}\subset\mathcal{F}_{1}$, where

Applying Lemma B.3 with $M=1$ and $\sigma=f_{\max}:=8G_{\max}F_{\max}$ implies that the following holds with probability at least $1-\delta$:

Dividing $\lambda_{n}$ on both sides gives

Let $x=J_{1}(\hat{H}^{\pi}_{n})$ and the above inequality becomes $x^{2}\leq a+bx$ for some $a,b>0$. When $a\leq bx$, we have $x^{2}\leq 2bx$, or $x^{2}\leq 4b^{2}$. When $a>bx$, we have $x^{2}\leq a+bx\leq 2a$. Thus $x^{2}\leq\max(4b^{2},2a)\leq 2a+4b^{2}$. Now we have

Now using (B.3), w.p. $1-4\delta$ for all $\pi\in\Pi$:

Let $C^{}(\delta)=\max(C_{4}(\delta)+C_{5}(\delta),1)$, $\beta_{1}=\frac{\alpha}{1+\alpha}$ and $\omega_{1}=0$. We have shown that with probability at least $1-4\delta$, the inequalities (B.2), (B.3) and the followings hold:

Rate Improvement Let $c=4,\beta=\beta_{1},\omega=\omega_{1}$ and $C(\delta)=C^{}(\delta)$. Denote by $E_{n}$ the event that the inequalities (B.3), (B.2) holds and

We have shown that $\Pr(E_{n})\geq 1-c\delta$. Below we improve the rate by refining the bound of the remainder term, $\operatorname{Rem}(\pi)$. First, we note that for the constant, $\iota$, specified in the condition, under the event $E_{n}$,

see Lemma B.4 for the deviation of the second inequality and similarly,

Combing with (B.2), which holds under the event $E_{n}$, we have

Thus using the same argument as in the proof of Lemma B.2 gives

We now obtain that the following holds w.p. $1-(c+1)\delta$, for all $\pi\in\Pi$

Thus the convergence rate is improved to $\beta_{2}=\frac{\beta_{1}+1/(1+\alpha)}{2}$ and $\omega_{2}=(1+\omega_{1})/2$. The same procedure can be applied $k$ times. It is easy to verify that for any $k\geq 1$, $\beta_{k+1}=\frac{\beta_{k}+1/(1+\alpha)}{2}=\frac{1}{1+\alpha}-\frac{(1-\alpha)2^{-k}}{1+\alpha}$ and $\omega_{k+1}=(1+\omega)/2=1-2^{-k}$ and thus the desired result. ∎

Recall the ratio estimator, $\hat{\omega}^{\pi}_{n}$ in (4.10). From Theorem B.3 and Lemma B.1, w.p. $1-(3+k)\delta$ for all $\pi\in\Pi$ we have

where, as in the beginning of the proof of Theorem B.3, ${K_{1}}={K_{0}}L(1+2C_{2}^{2})$ and ${K_{0}}$ is the leading constant in Lemma B.1. As such, we have

For simplicity, let $f(g)(D)=(1/T)\sum_{t=1}^{T}g(S_{t},A_{t})$. On the other hand,

where $\mathcal{G}_{1}=\{f(g):g\in\mathcal{G},J_{2}(g)\leq 1\}$. Using Lemma B.1, w.p $1-\delta$ for all $\pi$,

where ${K_{0}}$ is the leading constant in Lemma B.1. Combining with the bound on $J_{1}(\hat{H}^{\pi}_{n})$ in Theorem B.3 gives that

Together with the bound on $\|\hat{e}^{\pi}_{n}-e^{\pi}\|$, we have

Recall that $\|\omega^{\pi}\|^{2}=\int\omega^{\pi}(s,a)d^{\pi}(s,a)>1$. As a result, $Pf(e^{\pi})=\frac{1}{\|\omega^{\pi}\|^{2}}=\|e^{\pi}\|^{2}$. Finally we have

where the leading constant only depends on $F_{\max},G_{\max},L_{\Theta},\text{diam}(\Theta),C_{3},\alpha$

For $g_{1},g_{2}\in\mathcal{G},q\in\mathcal{F},\pi\in\Pi$, we introduce

Thus we have $\|\hat{g}_{n}^{\pi}(q)-g^{*}_{\pi}(H)\|^{2}+\|\hat{g}_{n}^{\pi}(q)-g^{*}_{\pi}(H)\|_{n}^{2}+\mu_{n}J_{2}^{2}(\hat{g}_{n}^{\pi}(q))=I_{1}(q,\pi)+I_{2}(q,\pi)$, where

For the first term, the optimizing property of $\hat{g}_{n}^{\pi}(q)$ implies that

Thus, $I_{1}(q,\pi)\leq 5\mu_{n}J_{2}^{2}(g^{*}_{\pi}(H))+2\mu_{n}J_{1}^{2}(H)$ holds for all $(q,\pi)$.

where we introduce $f=f_{1}-f_{2}$. Fix some $t>0$.

Next we verify the conditions (A1 - A4) in Theorem 19.3 in (Györfi et al. 2006) with $\mathcal{F}=\mathcal{F}_{l}$, $\epsilon=1/2$ and $\eta=2^{l}t$ to get an exponential inequality for each term in the summation, similar to the proof of Lemma B.2 below.

The (A1) and (A2) conditions are easy to verify using Assumptions (4-1), (6-1) and (6-2). For (A1), it is easy to see that

and thus $K_{1}=6G_{\max}(G_{\max}+1+2F_{\max})$. For (A2), note that

To ensure the condition (A3) holds for every $l$, i.e., $\sqrt{n}\epsilon\sqrt{1-\epsilon}\sqrt{\eta}\geq 288\max(K_{1},\sqrt{2K_{2}})$ (recall $\epsilon=1/2$ and $\eta=2^{l}t$), we just need to ensure the inequality holds for $l=0$, i.e.

That is, $t\geq c_{1}n^{-1}$ where $c_{1}=8(288\max(K_{1},\sqrt{2K_{2}}))^{2}$.

Next we verify the condition (A4). It is straightforward to see that with $M=\sqrt{\frac{2^{l}t}{\mu_{n}}}$, we have

where $\eta(\Pi,\mathcal{F}_{M})=\{s\mapsto\sum_{a}\pi(a|s)q(s,a):H\in\mathcal{F}_{M},\pi\in\Pi\}$ is a class of state-only function depending on the policy class and the function class $\mathcal{F}_{M}$. Let $c_{2}=6+12F_{\max}+24G_{\max}$. As a result of the entropy condition in Assumption (4-4), we have

Or equivalently $x^{\frac{1+\alpha}{2}}\geq 4\cdot 96\sqrt{2}\max(K_{1},2K_{2})\sqrt{c_{3}}\left(\frac{2^{l}t}{\mu_{n}}\right)^{\alpha/2}n^{-1/2}$. Clearly we only need to ensure the inequality holds when $x$ is at the minimum. That is, below is sufficient for the condition (A4) to hold:

To ensure the above holds for all $l\geq 0$, we require $t$ to satisfy

Or, simply requiring $t\geq c_{4}(\mu_{n}^{\alpha}n)^{-1}$ where $c_{4}=c_{3}18(32)^{4}\max(K_{1}^{2},4K_{2}^{2})$.

To summarize, the conditions (A1-A4) would be satisfied for every $l$ as long as $t\geq c_{1}n^{-1},t\geq\mu_{n}$ and $t\geq c_{4}(\mu_{n}^{\alpha}n)^{-1}$. Applying Theorem 19.3 in (Györfi et al. 2006) for each term implies that

where $c_{5}=8\cdot 128\cdot 2304\max(K_{1}^{2},K_{2})$. For any $\delta>0$, when $t\geq\log(120/\delta)c_{5}n^{-1}$, we have both $\exp(-\frac{nt}{c_{5}})\leq 1/2$ and $120\exp(-nt/c_{5})\leq\delta$ and as a result

Collecting all the conditions on $t$ and combing with the bound of $I_{1}(q,\pi)$, we have shown that w.p. at least $1-\delta$, the following holds for all $q,\pi$:

where the leading constant can be chosen by $K=5+8(288\max(K_{1},\sqrt{2K_{2}}))^{2}+6\cdot 8\cdot 128\cdot 2304\max(K_{1}^{2},K_{2})+2(18(32)^{4}\max(K_{1}^{2},4K_{2}^{2}))(6+12F_{\max}+24G_{\max})^{2\alpha}\cdot(4C_{1}+(\text{diam}(\Theta)L_{\Theta}F_{\max})^{2\alpha}/(2\alpha))$.

Suppose Assumptions (4-1), (4-2), (6-1), (6-2), (6-4), (4-4) and (3-2) hold. Then, the following hold with probability at least $1-{2\delta}$: for all policy $\pi\in\Pi$:

and the remainder term, $\operatorname{Rem}(\pi)$ is given by

For $g_{1},g_{2}\in\mathcal{G}$, define the functionals $\boldsymbol{f}_{1},\boldsymbol{f}_{2}$,

Using the optimizing property of $\hat{H}^{\pi}_{n}$ in (4.8), the term in the first parentheses can be bounded

where in the last equality we use the fact that $g_{1}^{2}+(g_{2}-2g_{3})g_{2}=g_{1}^{2}+g_{2}^{2}-2g_{2}g_{3}=(g_{1}-g_{2})^{2}+2g_{1}g_{2}-2g_{2}g_{3}=(g_{1}-g_{2})^{2}+2g_{2}(g_{1}-g_{3})$. In summary, we have

Below we provide the upper bound for each of the three terms. Recall that by Lemma B.1, the event, $E_{n}$ holds with probability at least $1-\delta$. Let the leading constant specified in Lemma B.1 be ${K_{0}}$.

Step I: bounding $I_{1}(\pi)$ Under the event $E_{n}$, we have

In addition, under the event $E_{n}$, we have

For simplicity, let $\beta(n,\mu_{n},\delta,p)=p^{1/2}n^{-1/2}\mu_{n}^{-(1+\alpha)/2}+(1+\sqrt{\log(1/\delta)})(n\mu_{n})^{-1/2}$. Under $E_{n}$, we have

Now we have $\Pr(\exists\pi\in\Pi,I_{2}(\pi)>t)\leq\Pr(\{\exists\pi\in\Pi,I_{2}(\pi)>t\}\cap E_{n})+\delta$ and we bound the first term using peeling device on $\lambda_{n}J_{1}^{2}(\hat{H}^{\pi}_{n})$ in $I_{2}(\pi)$:

where $\mathcal{F}_{l}=\{\boldsymbol{f}_{1}(g):J_{2}(g)\leq c_{1}(1+\sqrt{(2^{l+1}t)/\lambda_{n}}+\beta(n,\mu_{n},\delta,p)),g\in\mathcal{G}\}$. In what follows we verify the conditions (A1-A4) in Theorem 19.3 in (Györfi et al. 2006) with $\mathcal{F}=\mathcal{F}_{l}$, $\epsilon=1/2$ and $\eta=2^{l}t$ to get an exponential inequality for each term in the summation.