EconBase
← Back to paper

Off-Policy Evaluation of Bandit Algorithm from Dependent Samples under Batch Update Policy

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

50,133 characters · 19 sections · 56 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Off-Policy Evaluation of Bandit Algorithm from Dependent Samples under Batch Update Policy

abstractThe goal of off-policy evaluation (OPE) is to evaluate a new policy using historical data obtained via a behavior policy. However, because the contextual bandit algorithm updates the policy based on past observations, the samples are not independent and identically distributed (i.i.d.). This paper tackles this problem by constructing an estimator from a martingale difference sequence (MDS) for the dependent samples. In the data-generating process, we do not assume the convergence of the policy, but the policy uses the same conditional probability of choosing an action during a certain period. Then, we derive an asymptotically normal estimator of the value of an evaluation policy. As another advantage of our method, the batch-based approach simultaneously solves deficient support problem. Using benchmark and real-world datasets, we experimentally confirm the effectiveness of the proposed method.

Introduction

As an instance of sequential decision-making problems, the multi-armed bandit (MAB) algorithms have attracted significant attention in various applications, such as ad optimization, personalized medicine, search engines, and recommendation systems. Recently, various methods for evaluating a new policy using historical data obtained via the MAB algorithms kdd2009_ads,li2010contextual have emerged. The goal of off-policy evaluation (OPE) is to evaluate a new policy by estimating the expected reward obtained from the new policy dudik2011doubly,wang2017optimal,narita2019counterfactual,pmlr-v97-bibaut19a,Kallus2019IntrinsicallyES,Oberst2019. Although an OPE algorithm estimates the expected reward from a new policy, most existing studies presume that the samples are independent and identically distributed (i.i.d.). However, the MAB algorithm policy updates the probability of choosing an action based on past observations, and samples are not i.i.d. owing to this update. In this case, such existing studies do not guarantee that their estimators have asymptotic normality and $\sqrt{T}$-consistency for a sample size $T$. Therefore, there is a strong motivation to establish a novel method for OPE from dependent samples.

Several pioneering studies address OPE from dependent samples Laan2008TheCA,Hahn2011,Laan2016onlinetml,Luedtke2016,hadad2019,Kato2020. We can group the methods for deriving asymptotic normality into the following three groups: (a) Laan2008TheCA, Laan2016onlinetml, hadad2019, and Kato2020 derive the asymptotic normality with the central limit theorem (CLT) of a martingale difference sequence (MDS) by assuming that the probability of choosing an action converges to a time-invariant probability; (b) Luedtke2016 derives the asymptotic normality by standardizing a MDS, which is also used in the first group; (c) Hahn2011 apply asymptotic theory for the batched probability update process.

This paper focuses on an approach of the third group; that is, there are sufficiently large sample sizes in each batch. Compared with Hahn2011, our proposed method is more general and applicable in practical applications. Our method has the following three advantages compared with existing studies: (i) it does not assume convergence of the probability of choosing an action; (ii) it allows the probability of choosing an action to be $0$ for some actions in some batches; (iii) we can also use non-Donsker nuisance estimators as well as Laan2016onlinetml.

This paper has three main contributions. First, we provide a solution for OPE from dependent samples obtained via the MAB algorithms. Second, under the batch update policy, the proposed estimator achieves the asymptotic normality with fewer assumptions. Third, the estimator also experimentally shows a lower mean squared error (MSE) in some cases.

Problem Setting

Here, we formulate OPE under a batch update.

Date-Generating Process

Let $A_t$ be an action taking variable in $\mathcal{A}=\{1,2,\dots,K\}$, $X_t$ the covariate observed by the decision maker when choosing an action , and $\mathcal{X}$ the space of covariate. Let us denote a random variable of a reward at period $t$ as $Y_t=\sum^K_{a=1}\mathbbm{1}[A_t = a]Y_t(a)$, where $Y_t:\mathcal{A}\to\mathbb{R}$ is a potential outcome\footnote{We can express the DGP without using the potential outcome variable kato_uehara_2020.}. In this paper, we have access to a dataset $\{(X_t, A_t, Y_t)\}^{T}_{t=1}$ with the following data-generating process (DGP):

align[align omitted — 116 chars of source]

where $\Omega_{t-1}\in\mathcal{M}_{t-1}$ denotes the history until $t-1$ period defined as $\Omega_{t-1}=\{X_{t-1}, A_{t-1}, Y_{t-1}, \dots, X_{1}, A_1, Y_{1}\}$ with the space $\mathcal{M}_{t-1}$, $p(x)$ denotes the density of the covariate $X_t$, $\pi_t(a\mid x, \Omega_{t-1})$ denotes the probability of choosing an action $A_t$ conditioned on $X_t$ and $\Omega_{t-1}$, and $p(y\mid a, x)$ denotes the density of an outcome $Y_t$ conditioned on $A_t$ and $X_t$. We assume that $p(x)$ and $p(y\mid a, x)$ are invariant across periods, but $\pi_t(a\mid x, \Omega_{t-1})$ can take different values across periods. Let us call a policy inducing $\pi_t(a\mid x, \Omega_{t-1})$ a behavior policy.

Off-Policy Evaluation

This paper considers estimating the value of an evaluation policy using samples obtained under the behavior policy. Let an evaluation policy $\pi^\mathrm{e}:\mathcal{A}\times\mathcal{X}\to[0,1]$ be a probability of choosing an action $A_t$ conditioned on a covariate $X_t$. We are interested in estimating the expected reward from any pre-specified evaluation policy $\pi^{\mathrm{e}}(a \mid x)$. Then, we define the expected reward under an evaluation policy as $R(\pi^\mathrm{e}) := \mathbb{E}\left[\sum^K_{a=1}\pi^\mathrm{e}(a \mid x)Y_t(a)\right]$. For brevity, we also denote $R(\pi^\mathrm{e})$ as $\theta_0$. The goal of this paper is to estimate $\theta_0=R(\pi^{\mathrm{e}})$ using dependent samples under a batch update policy. To identify $\theta_0$, we assume overlaps in policy and the boundedness of the outcome.

assumptionThere exists a constant $C_1$ such that $0\leq \frac{\pi^\mathrm{e}(a\mid x)}{\pi_t(a\mid x, \Omega_{t-1})}\leq C_1$.
assumptionThere exists a constant $C_2$ such that $|Y| \leq C_2$.
remark[Existing Methods for OPE] We review three types of standard estimators of $R(\pi^{\mathrm{e}})$ under the case where $\pi_1(a\mid x, \Omega_{0})=\pi_2(a\mid x, \Omega_{1})=\cdots=\pi_T(a\mid x, \Omega_{T-1})=p(a\mid x)$ in the DGP defined in ((ref)). The first estimator is an inverse probability weighting (IPW) estimator given by $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\frac{\pi^\mathrm{e}(a\mid X_t)\mathbbm{1}[A_t=a]Y_t}{p(a\mid X_t)}$ rubin87,hirano2003efficient,swaminathan15a. Although this estimator is unbiased when the behavior policy is known, it suffers from high variance. The second estimator is a direct method (DM) estimator $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\hat{f}_{T}(a, X_t)$, where $\hat{f}_{T}(a, X_t)$ is an estimator of $f^*(a, X_t)$ HahnJinyong1998OtRo. This estimator is known to be weak against model misspecification for $f^*(a, X_t)$. The third estimator is an augmented IPW (AIPW) defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\Bigg(\frac{\pi^\mathrm{e}(a\mid X_t)\mathbbm{1}[A_t=a]\big(Y_t - \hat{f}_{T}(a, X_t)\big)}{p(a\mid X_t)} + \pi^\mathrm{e}(a\mid X_t)\hat{f}_{T}(a, X_t)\Bigg)$ robins94,ChernozhukovVictor2018Dmlf. Under certain conditions, it is known that this estimator achieves the efficiency bound (a.k.a semiparametric lower bound), which is the lower bound of the asymptotic MSE of OPE among regular $\sqrt{T}$-consistent estimators VaartA.W.vander1998As.
remark[Semiparametric Lower Bound] The lower bound of the variance is defined for an estimator of $\theta_0$ under some posited models of the DGP. If this posited model is a parametric model, it is equal to the Cram\'{e}r-Rao lower bound. When this posited model is a semiparametric model, we can define a corresponding Cram\'{e}r-Rao lower bound bickel98. narita2019counterfactual gives the semiparametric lower bound of the DGP ((ref)) under $\pi_1(a\mid x, \Omega_{0})=\cdots=\pi_T(a\mid x, \Omega_{T-1})=p(a\mid x)$ as $\mathbb{E}\Big[\sum^{K}_{a=1}\frac{\big(\pi^\mathrm{e}(a\mid X)\big)^2v^*(a, X_t)}{p(a\mid X)}+\left(\sum^{K}_{a=1}\pi^\mathrm{e}(a\mid X)f^*(a, X_t) - \theta_0\right)^2\Big]$.

\paragraph{Notations:} Let us denote $\mathbb{E}[Y_t(a)\mid x]$ and $\mathrm{Var}(Y_t(a)\mid x)$ as $f^*(a, x)$ and $v^*(a, x)$, respectively. Let $\mathcal{F}$ be the class of $f^*(a, x)$. Let $\hat{f}_{t}(a, x\mid \Omega_{t-1})$ be an estimator of $f^*(a, x)$ constructed from $\Omega_{t-1}$, respectively. Let $\mathcal{N}(\mu, \mathrm{var})$ be the normal distribution with the mean $\mu$ and the variance $\mathrm{var}$. For a random variable $Z$ and function $\mu$, let $\|\mu(Z)\|_2=\int |\mu(z)|^2 p(z) dz $ be the $L^{2}$-norm.

Patterns of Probability Update

In this paper, based on the update of $\pi_t(a\mid x, \Omega_{t-1})$, we classify the policies into two patterns, sequential update policy and batch update policy. For the sequential update policy, the policy updates $\pi_t(a\mid x, \Omega_{t-1})$ at each period Laan2008TheCA,Laan2016onlinetml,Kato2020. Under the batch update policy, after the policy continues using a fixed probability $\pi_t(a\mid x, \Omega_{t-1})$ for some periods without updates, the policy updates $\pi_t(a\mid x, \Omega_{t-1})$ Hahn2011,narita2019counterfactual. Although the sequential update is standard in the MAB problem, we often apply batch updates in industrial applications such as ad-optimization narita2019counterfactual. For OPE under the sequential update, Laan2008TheCA, Laan2016onlinetml, Luedtke2016, hadad2019, and Kato2020 proposed estimators with the asymptotic normality. For instance, an adaptive AIPW (A2IPW) estimator Laan2016onlinetml,hadad2019,Kato2020 has asymptotic normality if the behavior policy converges. On the other hand, we consider OPE under a batch update. Let $M$ denote the number of updates and $\tau\in I = \{1,2,\dots,M\}$ denotes the batch index. For $\tau\in I$, the probability is updated at a period $t_\tau$, where $t_\tau - t_{\tau-1} = Tr_\tau$, using samples $\{(X_t, Y_t, A_t)\}^{t_{\tau}}_{t=t_{\tau-1}}$, where $r_1 + r_2 + \cdots + r_M = 1$ and $t_0 = 0$. Thus, in addition to the DGP ((ref)), we assume

align*[align* omitted — 147 chars of source]

where $\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})$ denotes the assignment probability updated based on samples until the period $t_{\tau-1}$.

comment\begin{remark}[Policy and Probability] To transform a function $\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})$ to a probability, we consider the following process. First, we draw a random variable $\xi_t$ following the uniform distribution on $[0,1]$ at the $t$-th round. Then, at each round $t$, we select an action $A_t$ such that $\sum^{A_t-1}_{a=0}\pi_\tau(a\mid X_t, \Omega_{t_{\tau-1}}) \leq \xi_t \leq \sum^{A_t}_{a=0}\pi_\tau(a\mid X_t, \Omega_{t_{\tau-1}})$, where $\pi_t(0 \mid X_t, \Omega_{t-1})=0$. Under this process, we regard the probability $p_t(a\mid x)$ as $\pi_\tau(a\mid x, \Omega_{t_{\tau-1}})$. \end{remark}

Related Work

For the sequential update, Laan2008TheCA, Laan2016onlinetml, hadad2019, and Kato2020 assume that the probability of choosing an action converges to a time-invariant function almost certainly; that is, $\pi_t(a\mid x, \Omega_{t-1})\xrightarrow{\mathrm{p}} \alpha(a\mid x)$, where $\alpha:\mathcal{X}\to (0,1)$. This assumption enables us to apply the CLT for MDS. Laan2016onlinetml proposed constructing step-wise nuisance estimators, which enables us to derive asymptotic normality without Donsker's conditions of nuisance estimators. This technique is a generalization of sample-splitting, which is also called cross-fitting in a context of double/debiased machine learning klaassen1987,ZhengWenjing2011CTME,ChernozhukovVictor2018Dmlf. For the A2IPW estimator, hadad2019 proposed using an adaptive weight for stabilizing the behavior, and Kato2020 derived concentration inequality based on the law of iterated logarithms. On the other hand, we also construct an MDS and apply the CLT, but do not assume $\pi_t(a\mid x, \Omega_{t-1})\xrightarrow{\mathrm{p}} \alpha(a\mid x)$ by using batch update policy. Instead, we assume a sufficient sample size for each batch.

For such a non-stationary setting, Luedtke2016 also proposed an estimator with asymptotically normality for sequential update policy without using batch update policy. For deriving the asymptotic normality, Luedtke2016 used standardization for an MDS. Although the method enables us to construct an asymptotically normal estimator for various estimators, the proposed estimator only has $\sqrt{T - \ell}$-consistency for another sample size $\ell > 0$, not $\sqrt{T}$, to estimate the variances of the MDS.

As other related work, in the MAB problem, perchet2016 considered the setting of batch policy updates. In OPE, narita2019counterfactual also discuss a similar problem setting, but they assume that samples are i.i.d. Independently, kelly2020 provided a method for deriving a confidence interval of an ordinary least squares estimator, which is a different parameter of what we want to estimate.

OPE under Batch Update Policy

This section introduces a concept for conducting OPE under the batch update policy and a method based on the concept with its theoretical properties.

Strategy for OPE

For OPE under batch update policy, we consider asymptotic properties based on the assumption of $t_\tau - t_{\tau-1} \to \infty$ as $T\to \infty$ for fixed $\tau$. Because $\big\{(X_t, A_t, Y_t)\big\}^{t_\tau}_{t=t_{\tau-1}}$ is i.i.d., we can use the standard limit theorems for the partial sum of the samples to obtain an asymptotically normal estimator of $\theta_0 = R(\pi^\mathrm{e})$. However, we also have the motivation to use all samples together to increase the efficiency of the estimator. Therefore, based on the idea of generalized method of moments (GMM), we propose an estimator of $\theta_0$ considering the sample averages of each block as an empirical moment conditions. The main difference from the standard GMM is the assumption that the samples are not i.i.d. However, for the case under the batch update, we can apply the central limit theorems (CLT) for the martingale difference sequences (MDS) by appropriately constructing an estimator. We describe the proposed method as follows.

Estimator for OPE

We propose an estimator of $\theta_0$ based on a idea of GMM. For an index of batch $\tau\in I$, a function $f\in\mathcal{F}$ such that $f:\mathcal{A}\times \mathcal{X} \to \mathbb{R}$ and an evaluation policy $\pi^\mathrm{e} \in \Pi$, we define $h^{\mathrm{OPE}}_t:\mathcal{X}\times\mathcal{A}\times\mathbb{R}\times I \times \Theta\times\mathcal{F}\times\Pi \to \mathbb{R}$ as $h^{\mathrm{OPE}}_t(x, k, y; \tau, \theta, f, \pi^\mathrm{e}) = \frac{1}{r_\tau}\eta_t(x, k, y; \tau, \theta, f, \pi^\mathrm{e})\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]$, where $\eta_t(x, k, y; \tau, \theta, f, \pi^\mathrm{e}) := \phi_t(x, k, y; \tau, f, \pi^\mathrm{e}) - \theta$ and $\phi_t(x, k, y; \tau, f, \pi^\mathrm{e}):=$

align*[align* omitted — 165 chars of source]

Let us note that, for $\tau\in I$, $\theta_0\in \Theta$, $f_{t-1}\in \mathcal{F}$, and $\pi^\mathrm{e}\in\Pi$, the sequence $\big\{h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; \tau, \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e})\big\}^{T}_{t=1}$ is an MDS: for $h^{\mathrm{OPE}}_t(X_t, A_t, Y_t; \tau, \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e})$, by $\mathbb{E}[\mathbbm{1}[A_t=a]\mid X_t, \Omega_{t-1}]=\pi_{\tau}(a\mid X_t, \Omega_{t_{\tau - 1}})$, we have

align*[align* omitted — 570 chars of source]

Let us also define $\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta, \hat{f}_{t-1}, \pi^\mathrm{e}\right):=$

align*[align* omitted — 289 chars of source]

Then, the sequence $\left\{\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\right\}^{T}_{t=1}$ is an MDS with respect to $\big\{\Omega_t\big\}^{T-1}_{t=0}$; that is, $\mathbb{E}\left[\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\mid \Omega_{t-1}\right] = \bm{0}$. Using the sequence $\left\{\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\right\}^{T}_{t=1}$, we define an estimator of OPE as $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) := $

align[align omitted — 191 chars of source]

where $\hat{\bm{q}}^{\mathrm{OPE}}_T(\theta) = \frac{1}{T}\sum^{T}_{t=1}\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta, \hat{f}_{t-1}, \pi^\mathrm{e}\right)$ and $\hat{W}_T$ is a data-dependent $(M\times M)$-dimensional positive semi-definite matrix. Let us note that the estimator defined in Eq. ((ref)) is an application of GMM with the moment condition $\bm{q}^{\mathrm{OPE}}(\theta_0) = \mathbb{E}\left[\frac{1}{T}\sum^T_{t=1}\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\right]=0$. For the minimization problem defined in Eq. ((ref)), we can analytically calculate the minimizer as $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) = w^\top_T D_T(\pi^\mathrm{e})$, where $w_T = (w_{T, 1}\ \cdots\ w_{T, M})^\top$ is an $M$-dimensional vector such that $\sum^M_{\tau=1}w_{T, \tau}=1$, and $D_T(\pi^\mathrm{e})$ is

align*[align* omitted — 350 chars of source]

We call the estimator a Batch-based Adaptive AIPW (BA2IPW) estimator. In Appendix (ref), we discuss the GMM perspective in more detail.

Asymptotic Properties

Here, we show the consistency and asymptotic normality of the proposed BA2IPW estimator $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})$.

theorem[Consistency of the BA2IPW Estimator] Suppose that there exists a constant $C_f > 0$ such that $\big|f_{t_{\tau-1}}(a, x)\big| < C_f$ for $\tau \in I$. Then, under Assumptions (ref) and (ref), $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})\xrightarrow{\mathrm{p}}\theta_0$.
proofWe use the law of large numbers for an MDS from the boundedness of $\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)$, and we have $\frac{1}{T}\sum^T_{t=1}\bm{h}^{\mathrm{OPE}}_t\left(X_t, A_t, Y_t; \theta_0, \hat{f}_{t-1}, \pi^\mathrm{e}\right)\xrightarrow{\mathrm{p}}0$ (Proposition (ref) in Appendix (ref)). This result means that $D_T(\pi^\mathrm{e})\xrightarrow{\mathrm{p}}I\theta_0$, where $I = (1\ 1\ \cdots\ 1)^\top$ is an $M$-dimensional vector. Therefore, $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) = w_T D_T(\pi^\mathrm{e}) \xrightarrow{\mathrm{p}} w_T I \theta_0 = \theta_0$.
theorem[Asymptotic Distribution of the BA2IPW Estimator] Suppose that (i) $w_T = (w_{T,1}\ \cdots\ w_{T,M})^\top \xrightarrow{\mathrm{p}} w = (w_{1}\ \cdots\ w_{M})^\top$; (ii) $w_{T,\tau} > 0$ and $\sum^M_{\tau=1}w_{T,\tau} = 1$; (iii) $\hat{f}_{t-1}(a, x) \xrightarrow{\mathrm{p}} f^*(a, x)$ for all $a \in \mathcal{A}$ and $x \in \mathcal{X}$; (iv) There exists a constant $C_f > 0$ such that $\big|\hat{f}_{t-1}(a, x)\big| < C_f$. Then, under Assumptions (ref) and (ref), $\sqrt{T}\big(\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}) - \theta_0\big) \xrightarrow{\mathrm{d}}\mathcal{N}\big(0, \sigma^2\big)$, where $\sigma^2 = \sum^M_{\tau=1}w_{\tau}\sigma^2_{\tau}$ and $\sigma^2_{\tau}=\frac{1}{r_\tau}\mathbb{E}\Big[\sum^{K}_{a=1}\frac{\big(\pi^\mathrm{e}(a\mid X)\big)^2\nu^*(a, X)}{\pi_{\tau}(a\mid X, \Omega_{t_{\tau - 1}})} + \left(\sum^{K}_{a=1}\pi^\mathrm{e}(a\mid X)f^*(a, X) - \theta_0\right)^2\Big]$.

The proof is shown in Appendix (ref). Readers might consider that the use of MDS for deriving the asymptotic normality is unnecessary. We discuss the necessity of MDS in Appendix (ref). We can also define a corresponding Batch-based Adaptive IPW (BAdaIPW). For the BAdaIPW estimator, the variance of a batch $\tau$ is $\sigma^2_{\mathrm{IPW}, \tau}=\frac{1}{r_\tau}\mathbb{E}\left[\sum^{K}_{a=1}\frac{\big(\pi^\mathrm{e}(a\mid X)\big)^2\mathbb{E}[Y^2_t\mid X_t]}{\pi_{\tau}(a\mid X, \Omega_{t_{\tau - 1}})} - \theta^2_0\right]$.

remark[Construction of $f_{t_m}$ and Donsker Condition] As well as the cross-fitting of double/debiased machine learning proposed by klaassen1987,ZhengWenjing2011CTME,Laan2016onlinetml, and ChernozhukovVictor2018Dmlf, the proposed estimator does not require Donsker's condition for asymptotic normality. This property comes from the MDS as pointed by Laan2016onlinetml. On the other hand, because the samples are not independent, we cannot use the standard regression to obtain a consistent estimator of $f^*$. For example, yang2002 propose a nonparametric method for the bandit process under some mild conditions.

Weight of the Proposed Estimator

Next, we discuss the choice of weight $w_T$.

\paragraph{Equal Weight:} A naive choice is weighting the moment conditions equality; that is, $w_{T, \tau} = \frac{1}{M}$. In this case, the proposed estimator boils down to $I^\top D_T(\pi^\mathrm{e})$, which is almost the same as the A2IPW estimator. Although the estimator itself is similar to the A2IPW estimator, the theoretical guarantee for the asymptotic normality is different. While the A2IPW estimator uses the assumption that the policy converges to a time-invariant policy, the proposed BA2IPW estimator uses the assumption of the batch update. We call the BA2IPW estimator with the equal weight a Plain BA2IPW (PBA2IPW) estimator.

\paragraph{Efficient Weight:} First, we consider an efficient weight $w_T$ that minimizes the asymptotic variance of $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})$. As well as the standard GMM, the $\tau$-th element of the efficient weight is given as $w^*_{\tau}=\frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$ GVK126800421. Here, we use the orthogonality among moment conditions; that is, zero covariance. In this case, the asymptotic variance becomes $1/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$. Therefore, for gaining efficiency, we use a weight $\hat{w}_{T,\tau} = \frac{1}{ \hat{\sigma}^2_{T, \tau}}/\sum^M_{\tau'=1}\frac{1}{\hat{\sigma}^2_{T, \tau'}}$, where $\hat{\sigma}^2_{T, \tau}$ is an estimator of $\sigma^2_\tau$. We call the BA2IPW estimator with the efficient weight an Efficient BA2IPW (EBA2IPW) estimator.

remark[Estimation of $w^*$] If $\hat{\sigma}^2_{T, \tau}\xrightarrow{\mathrm{p}}\sigma^2_\tau$, we also have $\hat{w}_{T,\tau} \xrightarrow{\mathrm{p}}w^*_\tau$ from the continuous mapping theorem. In this paper, we propose two estimators defined as $\hat{\sigma}^2_{T, \tau}=\frac{1}{r_\tau T}\sum^{T}_{t=1}\left\{h^{\mathrm{OPE}}_t(x, a, y; \tau, \widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}), \hat{f}_{t-1}, \pi^\mathrm{e})\right\}^2\times$ $\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]$ and $\tilde{\sigma}^2_{T, \tau}=\frac{1}{r_\tau T}\sum^{T}_{t=1}\left\{h^{\mathrm{OPE}}_t(x, a, y; \tau, \widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e}), \hat{f}_{T}, \pi^\mathrm{e})\right\}^2\times$ $\mathbbm{1}\big[t_{\tau - 1} < t \leq t_{\tau} \big]$.

\paragraph{Weight for Numerical Stability:} As explained above, we can obtain an efficient weight that minimizes the asymptotic variance. However, as hadad2019 pointed out, when using time-variance estimators $f_t$, early variable $\phi_t(X_t, A_t, Y_t; \tau, f_t, \pi^\mathrm{e})$ might be unstable due to the existence of an inaccurate estimators $f_t$ at early stages. hadad2019 proposed an adaptive weight that puts more weight on later variables $\phi_t(X_t, A_t, Y_t; \tau, f_t, \pi^\mathrm{e})$ in the estimator. In our problem setting, we can also introduce such weights. However, unlike the adaptive weight proposed by hadad2019, which must be martingale, we do not require martingales on the weights. This property is a benefit of the batch update. For instance, for stabilization, we define a weight $\frac{1}{\ddot{\sigma}^2_{T, \tau}}/\sum^M_{\tau'=1}\frac{1}{\ddot{\sigma}^2_{T, \tau'}}$, where $\ddot{\sigma}^2_{T, \tau} = \hat{\sigma}^2_{T, \tau}+ \alpha\frac{1}{t_\tau}\sum^{t_{\tau}}_{t=t_{\tau-1}+1} \big(f_t(A_t, X_t) - f_T(A_t, X_t)\big)^2$ and $\alpha > 0$ is a constant. The first term is an efficient weight described above. The second term reflects the deviation between $f_t$ and $f_T$, which would be more accurate because it uses more samples.

Main Algorithm

As discussed in Section (ref), we can minimize the asymptotic variance of the proposed estimator $\widehat{R}^{\mathrm{BA2IPW}}_T(\pi^\mathrm{e})$ by choosing $w_T$ appropriately. However, to obtain $w_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$, which is the optimal weight that minimizes the asymptotic variance, we need a consistent estimator of $\theta_0$, which is what we want to estimate. On the other hand, we have a consistent estimator of $\theta_0$ without using an optimal weight matrix $w_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$. Based on these properties, we propose two-step estimation. First, using an arbitrary positive definite weight $w^{(0)}_T$, such as the identity matrix, we obtain an initial estimate $\widehat{R}^{\mathrm{BA2IPW}, (1)}_T(\pi^\mathrm{e})$. Then, using $\widehat{R}^{\mathrm{BA2IPW}, (1)}_T(\pi^\mathrm{e})$, we construct $w^{(1)}_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$. We can obtain an efficient estimator $\widehat{R}^{\mathrm{BA2IPW}, (2)}_T(\pi^\mathrm{e})$ of $\theta_0$, as discussed. More generally, we consider an algorithm with iteration such that after obtaining $\widehat{R}^{\mathrm{BA2IPW}, (i-1)}_T(\pi^\mathrm{e})$, we estimate $w^{(i-1)}_{T,\tau}\xrightarrow{\mathrm{p}} \frac{1}{\sigma^2_\tau}/\sum^M_{\tau'=1}\frac{1}{\sigma^2_{\tau'}}$ and obtain a next estimator $\widehat{R}^{\mathrm{BA2IPW}, (i)}_T(\pi^\mathrm{e})$ by using $w^{(i-1)}_{T,\tau}$. We refer to this algorithm with $N$-iterations as $N$-step BA2IPW estimation. We can use sufficiently large $N$ because, at each iteration, we only calculate the weighted average of the moment conditions using $\widehat{R}^{\mathrm{BA2IPW}, (i)}_T(\pi^\mathrm{e})$, which is not time-consuming. Although the asymptotic properties of the iterated estimator are the same as those of the two-step estimator, we report that the iteration improves the empirical performance in some cases. For $N\geq 2$, we summarize the $N$-Step BA2IPW Estimation in Algorithm (ref). We can use any method to construct $f_{t_\tau}$ and $\pi_\tau$ as long as they are consistent for bandit data and satisfy some regularity conditions needed for Theorem (ref).

algorithm[algorithm omitted — 762 chars of source]

Deficient Support Problem

As an application of BA2IPW, we consider an OPE without Assumption (ref), which assumes that there exists $C_1$ such that $0\leq \frac{\pi^\mathrm{e}(a\mid x)}{p_t(a\mid x)}\leq C_1$. Instead of Assumption (ref), we consider a situation in which we are allowed to change the support of actions in each batch. For example, in the first batch, we choose an action from a set $\{1,2,3\}$ with a probability larger than $0$, but we choose an action from a set $\{1,2,4\}$ with a probability larger than $0$ in the second batch. In this case, the probability of choosing the action $4$ is $0$ in the first batch, while the probability of choosing the action $3$ is $0$ in the second batch. This situation is a common in practice and called deficient support problem Sachdeva2020. For this problem, instead of Assumption (ref), we use the following assumption.

assumptionFor $a\in\{1,2,\dots,K\}$, there exist $\tau\in\{1,2,\dots,M\}$ and $C_1$ such that $0\leq \frac{\pi^\mathrm{e}(a\mid x)}{\pi_{\tau}(a\mid x, \Omega_{t_{\tau - 1}})}\leq C_1$.

Under this assumption, if $\pi_{\tau}(a\mid x, \Omega_{t_{\tau - 1}})>0$ for at least one batch, we are allowed to use $\pi_{\tau'}(a\mid x, \Omega_{t_{\tau' - 1}})=0$ for $\tau'\neq \tau$. With this assumption, we derive the asymptotic normality in Appendix (ref). Thus, our approach provides a new solution to this problem.

table*[table* omitted — 3,163 chars of source]

Estimation of the Behavior Policy

In the proposed BA2IPW method, we assume that the true behavior policy is known. However, in many real-world applications, the assumption does not hold. To solve this problem, by using an estimator $\hat{g}_{t-1}$ of $\pi_\tau$, which is constructed from $\Omega_{t-1}$ as well as $\hat{f}_{t-1}$, we also propose a Batch-based Adaptive Doubly Robust (BADR) as $\hat{R}^{\mathrm{BADR}}_T(\pi^\mathrm{e}) = w^\top_T \widetilde{D}_T(\pi^\mathrm{e})$, where $w_T = (w_{T, 1}\ \cdots\ w_{T, M})^\top$ is an $M$-dimensional vector such that $\sum^M_{\tau=1}w_{T, \tau}=1$, $\tilde{D}_T(\pi^\mathrm{e})=$

align*[align* omitted — 419 chars of source]

and $\tilde{\phi}_t(x, k, y; \tau, f, g, \pi^\mathrm{e}):=$

align*[align* omitted — 133 chars of source]

For the BADR estimator, we show the asymptotic normality as follows. The proof is shown in Appendix (ref).

theorem[Asymptotic Distribution of the BADR Estimator] Suppose that (i) $w_T = (w_{T,1}\ \cdots\ w_{T,M})^\top \xrightarrow{\mathrm{p}} w = (w_{1}\ \cdots\ w_{M})^\top$; (ii) $w_{T,\tau} > 0$ and $\sum^M_{\tau=1}w_{T,\tau} = 1$; (iii) for $\alpha\beta=\mathrm{o}_{p}((t-t_{\tau-1})^{-1/2}),\alpha=\mathrm{o}_{p}(1), \beta=\mathrm{o}_{p}(1)$, and each $\tau\in I$, the nuisance estimators satisfy $\|\hat{g}_{t-1}(a\mid X_t) - \pi_{\tau}(a\mid X_t, \Omega_{t_{\tau-1}})\|_{2}=\alpha$, and $\|\hat{f}_{t-1}(a,X_t)-f^*(a,X_t)\|_2=\beta$, where the expectation of the norm is over $X_t$; (iv) there exit constants $C_f$ and $C_g$ such that $|\hat{f}_{t-1}(a, x)| \leq C_f$ and $0 < \left|\frac{\pi^\mathrm{e}(a\mid x)}{\hat{g}_{t-1}(a\mid x)}\right| \leq C_g$ for all $a\in\mathcal{A}$ and $x\in\mathcal{X}$. Then, under Assumptions (ref) and (ref), $\sqrt{T}\big(\hat{R}^{\mathrm{BADR}}_T(\pi^\mathrm{e}) - \theta_0\big) \xrightarrow{\mathrm{d}}\mathcal{N}\big(0, \sigma^2\big)$.

Off-Policy Learning

An important application of OPE is Off-Policy Learning (OPL), which attempts to determine the optimal policy maximizing the expected reward. Let us define the optimal policy $\pi^*$ as $\pi^* = \operatorname*{arg\,max}_{\pi\in \Pi} R(\pi)$, where $\Pi$ is a policy class. By applying each OPE estimator, we estimate the optimal policy as $\hat{\pi} = \operatorname*{arg\,max}_{\pi\in \Pi}\widehat{R}^{\mathrm{BA2IPW}}_T(\pi)$.

table*[table* omitted — 1,257 chars of source]
table*[table* omitted — 1,474 chars of source]

Experiments

Using benchmark datasets and real-world logged data, we demonstrate the effectiveness of the BA2IPW estimator with an equal weight (PBA2IPW) and efficient weight using variance estimators $\hat{\sigma}^2_{T,\tau}$ (EBA2IPW) and $\tilde{\sigma}^2_{T,\tau}$ (EBA2IPW'), and BAdaIPW estimator with an equal weight (BAdaIPW). Note that although the forms of the several estimators are the same as the existing studies, the theoretical guarantees are different.

Experiments with Benchmark Dataset

Following dudik2011doubly and Chow2018, we evaluate the proposed estimators using classification datasets by transforming them into contextual bandit data. From the LIBSVM repository, we use the satimage, pendigits, mnist, letter, sensorless, and connect-4 datasets \footnote{\url{https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/}}. For a batched update behavior policy, we use the random walk (RW) and LinearUCB (UCB) sutton1998reinforcement,li2010contextual,Wei2011 policies. When using the RW policy, we first decide the probability of choosing an action from uniform distribution. Then, we add a noise $0.01\mathcal{N}(0,1)$ at each batch, i.e., the policy is a random walk. At each batch, we standardize the values of random walk to be probability, i.e., all values are positive and the sum is $1$. When using the UCB policy, we choose an estimated best arm $\hat{A}_t$ firstly. Then, we create an adaptive policy. Then, we construct a behavior policy as a policy that chooses $\hat{A}_t$ with probability $0.8$ and the other arms with equal probability. While the probability of choosing an action of the UCB policy converges, that of the RW policy does not converge.

For each dataset, we compare the performances of the following estimators of policy value: PBA2IPW, EBA2IPW, EBA2IPW', BAdaIPW, Adaptive DM estimator (AdaDM) defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\pi^\mathrm{e}(a\mid X_t)\hat{f}_{t}(a, X_t)$, an AIPW defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\frac{\pi^\mathrm{e}(a\mid X_t)\mathbbm{1}[A_t=a]\big(Y_t - \hat{f}_{T}(a, X_t)\big)}{\pi_t(a\mid X_t, \Omega_{t-1})} + \frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\pi^\mathrm{e}(a\mid X_t)\hat{f}_{T}(a, X_t)$, and DM estimator (DM) defined as $\frac{1}{T}\sum^T_{t=1}\sum^K_{a=1}\pi^\mathrm{e}(a\mid X_t)\hat{f}_{T}(a, X_t)$. When estimating $f^*$, we use the Nadaraya-Watson regression (NW) estimator yang2002.

\paragraph{MSEs:} To construct an evaluation policy, we create a deterministic policy $\pi^d$ by training a logistic regression classifier on historical data and set the output as $\pi^d$. Through experiments, the behavior policy $\pi_t$ is assumed to be known. More details are in Appendix (ref). Let us construct the evaluation policy $\pi^\mathrm{e}$ as a mixture of $\pi^d$ and the uniform random policy $\pi^u$, defined as $\pi^\mathrm{e} = 0.9\pi^d+0.1\pi^u$. We construct the evaluation policy $\pi^\mathrm{e}$ as a mixture of $\pi^d$ and the uniform random policy $\pi^u$ defined as $\pi^\mathrm{e} = 0.9\pi^d+0.1\pi^u$. We compare the MSEs of six estimators, the PBA2IPW, EBA2IPW, EBA2IPW', BAdaIPW, AdaDM, AIPW, and DM estimators. For estimating the weights of EBA2IPW and EBA2IPW' estimators, we iteratively estimate the wights and the value $\theta_0$ $10$ times. In each experiment, we have historical data with a sample size $T=1500$. When estimating $f^*$, we use the Nadaraya-Watson regression (NW) estimator yang2002.

The resulting MSEs and their standard deviations (SDs) over $100$ replications of each experiment are shown in Tables (ref) and (ref). In many cases, the proposed methods show the preferable the existing methods. When using the RW policy, the policy does not converge to a time-invariant policy. Therefore, the proposed method is theoretically preferable for the situation. More importantly, we can construct confidence intervals from the proposed methods, but cannot construct it from the A2IPW estimator. On the other hand, when using the UCB policy, the policy converges to a time-invariant policy, but the proposed methods still show higher performance for various datasets in some cases. We consider that this result is based on the fact that, even though the policy approaches to a time-invariant policy, the update is only allowed in batch, and does not converge sufficiently. In Appendix (ref), we show the additional results.

\paragraph{Coverage Ratio of Confidence Interval:} In Table (ref), we show the coverage ratio of the confidence intervals derived in the previous experiments together with the MSEs. The coverage ratio of the confidence interval is a percentage at which it covers the true value $\theta_0$ in the confidence interval. For the $100$ trials of the previous experiment, we calculate the coverage ratio (CR) of $95\%$ confidence interval, which is constructed as $\left[\hat{\theta}_T - 1.96\sqrt{\frac{\hat{\sigma}^2}{1500}}, \hat{\theta}_T + 1.96\sqrt{\frac{\hat{\sigma}^2}{1500}}\right]$, where $\hat{\theta}_T$ is an estimator of $\theta_0$ and $\hat{\sigma}^2$ is its estimated asymptotic variance. In the results, the PBA2IPW and BAdaIPW estimator shows CR close to $0.95$ in many cases. The BAdaIPW estimator does not require an estimator of $f^*$. Therefore, compared with the other estimators, it shows more preferable performances. EBA2IPW estimator requires several variance estimators, and their estimation error worsen the results compared with the PBA2IPW and BAdaIPW estimators.

\paragraph{OPL:} In the experiments of OPL, we compare the performances of estimated policy maximizing expected reward obtained from the PBA2IPW, EBA2IPW, BAdaIPW, AdaDM, AIPW, and DM estimators. We conducted $5$ trials for each experiment with $T=1500$. The resulting expected rewards over the evaluation data (RWDs) and the SDs are shown in Table (ref), where we highlight in bold the best two estimators.

Experiments with Real-World Data

We apply our estimators to evaluate a policy using the real-world dataset in CyberAgent Inc., which is the second-largest Japanese advertisement company with about $7$ billion USD market capitalization (as of August 2020). This company simultaneously runs Thompson sampling and uniformly random sampling to determine the design of advertisements. The Thompson sampling updates the parameter every $30$ minute; therefore, there are batches with samples obtained during the $30$ minute. We use the logged data produced by the algorithms to confirm the empirical performance of the proposed estimators. To check the performance, we calculate the estimation error between the estimates of the value of the uniformly random sampling policy estimated from the dataset obtained by the Thompson sampling and the observed average reward of the uniformly random policy. More details are shown in Appendix (ref). We apply the PBA2IPW, EBA2IPW, EBA2IPW', BAdaIPW, and AdaDM estimators. The results are shown in Table (ref) and Figure (ref). We show the Bias, MSE, and averaged confidence intervals. While the PBA2IPW and BAdaIPW estimators suffer the high variance, the EBA2IPW, EBA2IPW', and AdaDM estimators show the preferable performances. Although the AdaDM estimator also shows the effectiveness for this dataset, the BA2IPW estimator is theoretically more robust because it is consistent even if $\hat{f}_t$ does not converge to $f^*$.

table[table omitted — 608 chars of source]
figure[figure omitted — 175 chars of source]

Conclusion

This study presented solutions for causal inference from dependent samples obtained via a batch-based bandit algorithm. By using the asymptotic property in batch, we applied the CLT for an MDS to obtain the asymptotic normality without requiring the convergence of the assignment probability. Additionally, we showed that the proposed method is applicable when the support of arms is incomplete, which is a notorious problem in OPE. In experiments, the proposed batch-based estimators showed theoretically expected performances for benchmark and real-world datasets.

\onecolumn