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.
48,090 characters · 13 sections · 39 citation commands
Debiased Off-Policy Evaluation for Recommendation Systems
\email{[email removed]}
\email{[email removed]}
\email{[email removed]}
\ccsdesc[500]{Information systems Online advertising} \ccsdesc[500]{Computing methodologies Learning from implicit feedback}
Interactive bandit and reinforcement learning (RL) systems (e.g., ad/news/recommendation/search platforms, personalized education and medicine) produce log data valuable for evaluating and redesigning the systems. For example, the logs of a news recommendation system record which news article was presented and whether the user read it, giving the system designer a chance to make its recommendation more relevant li2010contextual.
Exploiting log data is, however, more difficult than conventional supervised machine learning: the result of each log is only observed for the action chosen by the system (e.g., the presented news) but not for all the other actions the system could have taken. Moreover, the log entries are biased in that the logs over-represent actions favored by the system.
A potential solution to this problem is an A/B test that compares the performance of counterfactual systems. However, A/B testing counterfactual systems is often technically or managerially infeasible, since deploying a new policy is time- and money-consuming, and entails a risk of failure.
This leads us to the problem of counterfactual (off-policy, offline) evaluation, where one aims to use batch data collected by a logging policy to estimate the value of a counterfactual policy or algorithm without deploying it. Such evaluation allows us to compare the performance of counterfactual policies to decide which policy should be deployed in the field. This alternative approach thus solves the above problem with the naive A/B test approach. Key prior studies include Strehl2010,li2010contextual,li2011unbiased,Dudik2014,wang2016optimal,swaminathan2017off, narita2018efficient,gilotte2018offline for bandit algorithms, and precup2000eligibility, Jiang16, Thomas16, liu2018representation, Farajtabar2018MoreRD, Irpan2019OffPolicyEV, kallus2019DRL, Uehara2020MWL for RL algorithms.
Method. For off-policy evaluation with log data of RL feedback, this paper develops and empirically implements a novel technique with desirable theoretical properties. To do so, we consider a class of RL algorithms, including contextual bandit algorithms as important special cases. This class includes most of the widely-used algorithms such as (deep) Q-learning, Actor Critic, contextual $\epsilon$-greedy, and Thompson Sampling, as well as their non-contextual analogs and random A/B testing. We allow the logging policy to be an unknown function of numerous potentially important state variables. This feature is salient in real-world applications. We also allow the evaluation target policy to be degenerate, again a key feature of real-life situations.
We consider an offline estimator for the expected reward from a counterfactual policy. Our estimator integrates a well-known Doubly Robust estimator (rotnitzky1995semiparametric and modern studies cited above) with “Double/Debiased Machine Learning” developed in econometrics and statistics Chernozhukov2018, Chernozhukov2018LR. Building upon these prior studies, we show that our estimator is “$\sqrt{N}$-consistent” in the sense that its prediction converges in probability to the true performance of a counterfactual policy at a rate of $1/\sqrt{N}$ as the sample size $N$ increases. Our estimator is also shown to be “asymptotically normal,” meaning that it has an approximately normal distribution as $N$ gets large. We also provide a consistent estimator of its asymptotic variance, thus allowing for the measurement of statistical uncertainty in our prediction. Moreover, for special cases in which data are generated by contextual bandit algorithms, our estimator has the lowest variance in a wide class of estimators, achieving variance reduction relative to standard estimators. Importantly, these properties hold even when the analyst does not know which among a large number of potentially important state variables are actually important.
Simulation Experiment. We evaluate the performance of our estimator by conducting an experiment in a slightly different version of the OpenAI Gym CartPole-v0 environment brockman2016openai. In this version, there are many more state variables than the original one and some of them are irrelevant to the reward. In this challenging environment, our estimator produces smaller mean squared errors than widely-used benchmark methods (Doubly Robust estimator in the spirit of Jiang16 and Thomas16 and Inverse Probability Weighting estimator).
Real-World Experiment. We empirically apply our estimator to evaluate and optimize the design of online advertisement formats. Our application is based on proprietary data provided by CyberAgent Inc., the second largest Japanese advertisement company with about 5 billion USD market capitalization (as of February 2020). This company uses randomly chosen bandit algorithms to determine the visual design of advertisements assigned to users. This A/B test of randomly choosing an algorithm produces logged data and the ground truth for the performance of alternative algorithms. We use this data to examine the performance of our proposed method. We use the log data from an algorithm to predict the click through rates (CTR) of another algorithm, and assess the accuracy of our prediction by comparing it with the ground truth. This exercise shows that our estimator produces smaller mean squared errors than widely-used benchmark methods (Doubly Robust estimator in the spirit of Dudik2014 and Inverse Probability Weighting estimator). This improvement is statistically significant at the 5% level. Importantly, this result holds regardless of whether we know the data-generating logging policy or not, which shows that our estimator can substantially reduce bias and uncertainty we face in real-world decision-making.
This improved performance motivates us to use our estimator to optimize the advertisement design for maximizing the CTR. We estimate how much the CTR would be improved by a counterfactual policy of choosing the best action (advertisement) for each context (user characteristics). This exercise produces the following bottom line: Our estimator predicts the hypothetical policy to statistically significantly improve the CTR by 30% (compared to the logging policy) in one of the three campaigns we analyze. Our approach thus generates valuable managerial conclusions.
We consider historical data from a Markov Decision Process (MDP) as a mathematical description of RL and bandit algorithms. An MDP is given by ${\mathcal M}=\langle {\mathcal S},{\mathcal A},P_{S_0},P_S,P_R\rangle$, where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space, $P_{S_0}: {\mathcal S}\rightarrow [0,1]$ is the initial state distribution, $P_S: {\mathcal S}\times {\mathcal A}\rightarrow \Delta({\mathcal S})$ is the transition function with $P_S(s'|s,a)$ being the probability of seeing state $s'$ after taking action $a$ given state $s$, and $P_R:{\mathcal S}\times {\mathcal A}\times \mathbb{R}\rightarrow [0,1]$ be the conditional distribution of the immediate reward with $P_R(\cdot|s,a)$ being the immediate reward distribution conditional on the state and action being $(s,a)$. We assume that the state and action spaces ${\mathcal S}$ and ${\mathcal A}$ are finite. Given $P_R$, we define the mean reward function $\mu:{\mathcal S}\times {\mathcal A}\rightarrow \mathbb{R}$ as $\mu(s,a)=\int r dP_R(r|s,a)$.
We call a function $\pi: {\mathcal S}\rightarrow \Delta({\mathcal A})$ a {\it policy}, which assigns each state $s\in {\mathcal S}$ a distribution over actions, where $\pi(a|s)$ is the probability of taking action $a$ when the state is $s$. Let $H=(S_0,A_0,R_0,...,S_T,A_T,R_T)$ be a {\it trajectory}, where $S_t$, $A_t$, and $R_t$ are the state, the action, and the reward in step $t$, respectively, and $T$ denotes the last step and is fixed. We say that a trajectory $H$ is generated by a policy $\pi$, or $H\sim \pi$ in short if $H$ is generated by repeating the following process for $t=0,...,T$: (1) When $t=0$, the state $S_t$ is drawn from the initial distribution $P_{S_0}$. When $t\ge 1$, $S_t$ is determined based on the transition function $P_S(\cdot|S_{t-1},A_{t-1})$. (2) Given $S_t$, the action $A_t$ is randomly chosen based on $\pi(\cdot|S_t)$. (3) The reward $R_t$ is drawn from the conditional reward distribution $P_R(\cdot|S_t,A_t)$.
Suppose that we observe {\it historical data} $\{H_i\}_{i=1}^N$ where trajectories are independently generated by a fixed {\it behavior policy} $\pi_b$, i.e., $H_i\sim \pi_b$ independently across $i$. The historical data is a collection of i.i.d. trajectories. Importantly, we allow the components of the data generating process ${\mathcal M}$ and $\pi_b$ to vary with the sample size $N$. Specifically, let ${\mathcal M}_N$ and $\pi_{bN}$ be the MDP and the behavior policy, respectively, when the sample size is $N$, and let $\mathbb{P}_N$ denote the resulting probability distribution of $H_i$. $\mathbb{P}_N$ is allowed to vary with $N$ in a way that the functions $P_{S_0N}$, $P_{RN}$ and $\pi_{bN}$ are high dimensional relative to sample size $N$ even when $N$ is large. In some RL problems, for example, there are a large number of possible states. To capture the feature that the number of states $|{\mathcal S}_N|$ is potentially large relative to sample size $N$, we may consider a sequence of $\mathbb{P}_N$ such that $|{\mathcal S}_N|$ is increasing with $N$. For the sake of notational simplicity, we make implicit the dependence of ${\mathcal M}$ and $\pi_{b}$ on $N$.
We assume that we know the state space ${\mathcal S}$ and the action space ${\mathcal A}$ but know none of the functions $P_{S_0}$, $P_S$ and $P_R$. In some environments, we know the function $\pi_b$ or observe the probability vector $(\pi_b(a|S_{it}))_{a\in {\mathcal A}, t=0,...,T}$ for every trajectory $i$ in the historical data. Our approach is usable regardless of the availability of such knowledge on the behavior policy.
With the historical data $\{H_i\}_{i=1}^N$, we are interested in estimating the discounted value of the {\it evaluation policy} $\pi_e$, which might be different from $\pi_b$: with $\gamma\in [0,1]$ as the discount factor, $$ V^{\pi_e}\coloneqq\mathbb{E}_{H\sim \pi_e}\left[\sum_{t=0}^T \gamma^t R_t\right]. $$
The estimation of $V^{\pi_e}$ involves estimation of the behavior policy $\pi_b$ (if unknown), the transition function $P_S$, and the mean reward function $\mu$. These functions may be high dimensional in the sense of having a large number of arguments, possibly much larger than the number of trajectories $N$. To handle the issue, we use Double/Debiased Machine Learning (DML) by Chernozhukov2018. DML is a general method for estimation and inference of semiparametric models in the presence of a high-dimensional vector of control variables, and is characterized by two key features, cross-fitting and Neyman orthogonality, which we will discuss in detail later. These two features play a role in reducing the bias that may arise due to the estimation of high-dimensional parameters, which makes DML suitable for the off-policy evaluation problem with potentially high dimensions of $\pi_b$, $P_S$ and $\mu$.
Before presenting our estimator, we introduce some notation. $H_t^{s,a}=(S_0,A_0,...,S_t,A_t)$ is a trajectory of the state and action up to step $t$. $\rho_t^{\pi_e}: ({\mathcal S}\times {\mathcal A})^{t+1}\rightarrow \mathbb{R}_+$ is the importance weight function:
This equals the probability of $H$ up to step $t$ under the evaluation policy $\pi_e$ divided by its probability under the behavior policy $\pi_b$. Viewing $\rho_t^{\pi_e}$ as a function of $\pi_b$, define $\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)$ as the value of $\rho_t^{\pi_e}(H_t^{s,a})$ with the true behavior policy $\pi_b$ replaced with a candidate function $\tilde \pi_b$: $$\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)\coloneqq \prod_{t'=0}^t\frac{\pi_e(A_{t'}|S_{t'})}{\tilde \pi_b(A_{t'}|S_{t'})}.$$ We can think of $\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)$ as the estimated importance weight function when we use $\tilde \pi_b$ as the estimate of $\pi_b$. By definition, $\rho_t^{\pi_e}(H_t^{s,a};\pi_b)=\rho_t^{\pi_e}(H_t^{s,a})$, where the left-hand side is $\rho_t^{\pi_e}(H_t^{s,a};\tilde \pi_b)$ evaluated at the true $\pi_b$ and the right-hand side is the true importance weight function. Finally, let $q_t^{\pi_e}:{\mathcal S}\times {\mathcal A}\rightarrow \mathbb{R}$ be the state-action value function under policy $\pi_e$ at step $t$, where $q_t^{\pi_e}(s,a)\coloneqq\mathbb{E}_{H\sim \pi_e}[\sum_{t'=t}^T\gamma^{t'-t}R_{t'}|S_t=s,A_t=a]$. Using $P_S$, $\mu$ and $\pi_e$, $q_t^{\pi_e}$ can be obtained recursively: let
and for $t=0,...,T-1$,
Our estimator is based on the following expression of $V^{\pi_e}$ Thomas16:
where for any candidate tuple $\tilde \eta=(\tilde \pi_b,\{\tilde q_t^{\pi_e}\}_{t=0}^T)$, we define
where $\rho_{-1}^{\pi_e}=1$. To give an intuition behind the expression, we arrange the terms as follows:
The first term is the well-known Inverse Probability Weighting (IPW) estimator. The second term serves as a control variate that has zero mean as long as we plug in the true $\pi_b$; the mean of the term remains zero regardless of the function we plug in for $\{\tilde{q}_t^{\pi_e}\}_{t=0}^T$. This is the key to make our estimator insensitive to $\{\tilde{q}_t^{\pi_e}\}_{t=0}^T$. We construct our estimator as follows.
Possible estimation methods for $\pi_b$, $\mu$ and $P_S$ in Step 2 are (i) classical nonparametric methods such as kernel and series estimation, (ii) off-the-shelf machine learning methods such as random forests, lasso, neural nets, and boosted regression trees, and (iii) existing methods developed in the off-policy policy evaluation literature such as representation balancing MDPs liu2018representation. These methods, especially (ii) and (iii), are usable even when the analyst does not know which among a large number of potentially important state variables are actually important. This DML estimator differs from the Doubly Robust (DR) estimator developed by Jiang16 and Thomas16 in that we use the cross-fitting procedure, as explained next.
{\bf A. Cross-Fitting.} The above method uses a sample-splitting procedure called {\it cross-fitting}, where we split the data into $K$ folds, take the sample analogue of Eq. ((ref)) using one of the folds ($I_k$) with $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ estimated from the remaining folds ($I_k^c$) plugged in, and average the estimates over the $K$ folds to produce a single estimator. Cross-fitting has two advantages. First, if we use instead the whole sample both for estimating $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ and for computing the final estimate of $V^{\pi_e}$ (the “full-data” variant of the DR estimator of Thomas16), substantial bias might arise due to overfitting Chernozhukov2018,newey2018crossfitting. Cross-fitting removes the potential bias by making the estimate of $V^{\pi_e}$ independent of the estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$.
Second, a standard sample splitting procedure uses a half of the data to construct estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ and the other half to compute the estimate of $V^{\pi_e}$ (the DR estimator of Jiang16 and the “half-data” variant of the DR estimator of Thomas16). In contrast, cross-fitting swaps the roles of the main fold ($I_k$) and the rest ($I_k^c$) so that all trajectories are used for the final estimate, which enables us to make efficient use of data.
{\bf B. Neyman Orthogonality.} There is another key ingredient important for DML to have desirable properties. The DML estimator is constructed by plugging in the estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$, which may be biased due to regularization if they are estimated with machine learning methods. However, the DML estimator is robust to the bias, since $\psi$ satisfies the {\it Neyman orthogonality} condition Chernozhukov2018. The condition requires that for any candidate tuple $\tilde \eta=(\tilde \pi_b,\{\tilde q_t^{\pi_e}\}_{t=0}^T)$, $$ \left. \frac{\partial \mathbb{E}_{H\sim \pi_b}[\psi(H;\eta+r(\tilde \eta-\eta))]}{\partial r}\right|_{r=0}=0, $$ where $\eta=(\pi_b,\{q_t^{\pi_e}\}_{t=0}^T)$ is the tuple of the true functions (see Appendix (ref) for the proof that DML satisfies this). Intuitively, the Neyman orthogonality condition means that the right-hand side of Eq. ((ref)) is locally insensitive to the value of $\pi_b$ and $\{q_t^{\pi_e}\}_{t=0}^T$. More formally, the first-order approximation of the bias caused by using $\tilde \eta$ instead of the true $\eta$ is given by
which is exactly zero by the Neyman orthogonality condition. As a result, plugging in noisy estimates of $\pi_b$ and $\{q_t^{\pi_e}\}_{t=0}^T$ does not strongly affect the estimate of $V^{\pi_e}$.
In contrast, IPW is based on the following expression of $V^{\pi_e}$: $V^{\pi_e}=\mathbb{E}_{H\sim \pi_b}[\psi_{{\rm IPW}}(H;\pi_b)]$, where for any candidate $\tilde \pi_b$, we define
The Neyman orthogonality condition does not hold for IPW: for some $\tilde \pi_b\neq \pi_b$, $ \left. \frac{\partial \mathbb{E}_{H\sim \pi_b}[\psi_{{\rm IPW}}(H;\pi_b+r(\tilde \pi_b-\pi_b))]}{\partial r}\right|_{r=0}\neq0. $ Therefore, IPW is not robust to the bias in the estimate of $\pi_b$.
Let $ \sigma^2 =\mathbb{E}_{H\sim \pi_b}[(\psi(H;\eta)-V^{\pi_e})^2]$ be the variance of $\psi(H;\eta)$. To derive the properties of the DML estimator, we make the following assumption.
Assumption (ref) (ref) assumes that the variance of $\psi(H;\eta)$ is nonzero and finite. Assumption (ref) (ref) states that the estimator $(\hat \pi_{b,k},\{\hat q_{t,k}^{\pi_e}\}_{t=0}^T)$ belongs to the set ${\mathcal T}_N$, a shrinking neighborhood of the true $(\pi_b,\{q_{t}^{\pi_e}\}_{t=0}^T)$, with probability approaching one. It requires that $(\hat \pi_{b,k},\{\hat q_{t,k}^{\pi_e}\}_{t=0}^T)$ converge to $(\pi_b,\{q_{t}^{\pi_e}\}_{t=0}^T)$ at a sufficiently fast rate so that the rate conditions in Assumption (ref) (ref) are satisfied.
The following proposition establishes the $\sqrt{N}$-consistency and asymptotic normality of $\hat V^{\pi_e}_{{\rm DML}}$ and provides a consistent estimator for the asymptotic variance. Below, $\rightsquigarrow$ denotes convergence in distribution, and $\stackrel{p}{\rightarrow}$ denotes convergence in probability.
The proof is an application of Theorems 3.1 and 3.2 of Chernozhukov2018, found in Appendix (ref). The above convergence result holds under any sequence of probability distributions $\{\mathbb{P}_N\}_{N\ge 1}$ as long as Assumption (ref) holds. Therefore, our approach is usable, for example, in the case where there are a growing number of possible states, that is, $|{\mathcal S}|$ is increasing with $N$.
When $T=0$, a trajectory takes the form of $H=(S_0,A_0,R_0)$. Regarding $S_0$ as a context, it is possible to consider $\{H_i\}_{i=1}^N$ as batch data generated by a contextual bandit algorithm. In this case, the DML estimator becomes
where $(\hat \pi_{b,k},\hat \mu_k)$ is the estimator of $(\pi_b,\mu)$ using the subset of data $\{H_i\}_{i\in I_k^c}$. This estimator is the same as the DR estimator of Dudik2014 except that we use the cross fitting procedure. Proposition 1 has the following implication for the contextual bandit case. Let $\sigma^2_R(s,a)=\int (r-\mu(s,a))^2dP_R(r|s,a)$.
The variance expression coincides with the “semiparametric efficiency bound” obtained by narita2018efficient, where the semiparametric efficiency bound is the smallest possible asymptotic variance among all consistent and asymptotically normal estimators. Hence $\hat V^\pi_{{\rm DML}}$ is the lowest variance estimator.
We demonstrate the effectiveness of our proposed estimator in the OpenAI Gym CartPole-v0 environment brockman2016openai. In this environment, the agent decides between two actions, moving the cart left or right, so that the pole attached to it stays upright. We make the following two changes to this environment to make it a more high-dimensional setting. First, the inputs to the agent in the original environment are 4 values representing the state (cart position, cart velocity, pole angle, and pole velocity at the tip), but we instead use a patch of the screen centered on the cart as an input, using Reinforcement Learning (DQN) Tutorial NEURIPS2019_9015.\footnote{\url{https://pytorch.org/tutorials/intermediate/reinforcement_q_learning.html}} We then represent the state ($S_t$) as the difference between the current screen patch and the previous one. This leads to an environment with a larger number of state variables. Second, at the beginning of every episode, we randomly select $10\%$ of pixels in the screen patch and black them out. We then move the pixels painted black horizontally over time. In this environment, some pixels that have nothing to do with the movement of the cart and pole become black or white during the episode, and the state variables representing these pixels are completely irrelevant to the reward. This mimics a situation where we do not know which among a large number of state variables are actually important, which makes the environment suitable for testing the performance of our method.
We use a greedy policy from a learned Q function as the evaluation policy ($\pi_e$) and an $\epsilon$-greedy policy as the behavior policy ($\pi_b$). Following Reinforcement Learning (DQN) Tutorial NEURIPS2019_9015, we use a convolutional neural network to learn the Q function. Our convolutional network has the following structure: a convolution layer with 3 $16 \times 16$ filters and stride 2, followed by a $16 \times 16$ batch normalization layer, followed by a convolution layer with 16 $32 \times 32$ filters and stride 2, followed by a $32 \times 32$ batch normalization layer, followed by a convolution layer with 32 $32 \times 32$ filters and stride 2, followed by a $32 \times 32$ batch normalization layer, followed by a dense layer, and finally the binary output layer. All hidden layers employ ReLU activation functions.
Our experiment proceeds as follows. Firstly, we run the behavior policy to collect 100 trajectories and use them to estimate the expected value of the evaluation policy. Secondly, we run the evaluation policy to observe the ground truth value. We iterate this procedure 100 times and then report the mean squared error (MSE) as a performance metric.
We compare DML (our proposal) with four baseline methods: Direct Method (DM), also known as Regression Estimator, IPW, full-data DR, and half-data DR. DM, IPW and full-data DR are given by
where $\psi$ and $\psi_{{\rm IPW}}$ are defined in Eq. ((ref)) and Eq. ((ref)), respectively, and $(\hat \pi_{b},\{\hat q_t^{\pi_e}\}_{t=0}^T)$ is the estimator of $(\pi_b,\{q_t^{\pi_e}\}_{t=0}^T)$ using the {\it whole} data. Half-data DR is given by $$ \hat V_{{\rm DR-half}}^{\pi_e}=\frac{1}{(N/2)}\sum_{i\in I_1} \psi(H_i;\hat\pi_{b,1}, \{\hat q_{b,1}^{\pi_e}\}_{t=0}^T), $$ where $I_1$ is a random subset of trajectory indices $\{1,...,N\}$ of size $N/2$, and $(\hat \pi_{b,1},\{\hat q_{t,1}^{\pi_e}\}_{t=0}^T)$ is the estimator of $(\pi_b,\{q_t^{\pi_e}\}_{t=0}^T)$ using the other half of the data. $\hat V_{{\rm DM}}^{\pi_e}$ and $\hat V^{\pi_e}_{{\rm IPW}}$ do not satisfy the Neyman orthogonality condition nor use cross-fitting. $\hat V^{\pi_e}_{{\rm DR-full}}$ and $\hat V^{\pi_e}_{{\rm DR-half}}$ satisfy the Neyman orthogonality condition but do not use cross-fitting. DML is therefore expected to outperform these four baseline methods.
We use a convolutional neural network to estimate $\pi_b$ and $\{q_t^{\pi_e}\}_{t=0}^T$. We use the same model structure as the one used to construct the evaluation policy and behavior policy. We set the learning rate to $0.0001$ for estimating the Q function and to $0.0002$ for estimating the behavior policy. The batch size is $128$ for both models.
Note that for DML, we hold out a fold and use the rest of the data to obtain the estimators of $\pi_b$ and $\{q_t^{\pi_e}\}_{t=0}^T$. We use $K = 2$ folds. We verify that increasing $K$ to 3 or 4 changes the results little.
The MSEs of the five estimators are reported in Table (ref). Regardless of the value of $\epsilon$ (the probability that the $\epsilon$-greedy algorithm chooses an action randomly), DML outperforms the other estimators. Since the CartPole-v0 environment tends to have long trajectories, the importance weight (Eq. ((ref))) tends to be very volatile due to the curse of the horizon. As a result, the MSE of IPW is much larger than the other four. When $\epsilon$ is smaller, the performance of full-data DR also seems to be severely affected by the importance weight and is worse than DM. As $\epsilon$ increases, this effect gets weaker, and full-data DR becomes better than DM. Half-data DR performs much better than DM and IPW for any $\epsilon$ thanks to Neyman Orthogonality and sample splitting, but performs worse than DML due to the lack of cross-fitting. This simulation result thus highlights the validity of the DML estimator even in the challenging situation where there are a large number of potentially important state variables.
We apply our estimator to empirically evaluate the design of online advertisements. This application uses proprietary data provided by CyberAgent Inc., which we described in the introduction. This company uses bandit algorithms to determine the visual design of advertisements assigned to user impressions in a mobile game.
Our data are logged data from a 7-day A/B test on ad “campaigns,” where each campaign randomly uses either a multi-armed bandit (MAB) algorithm or a contextual bandit (CB) algorithm for each user impression. The parameters of both algorithms are updated every 6 hours. As the A/B test lasts for 7 days, there are 28 batches in total. The parameters of both algorithms stay constant across impressions within each batch.
In the notation of our theoretical framework, $T$ is zero, and a trajectory takes the form of $H=(S_0,A_0,R_0)$. Reward $R_0$ is a click while action $A_0$ is an advertisement design. The number of possible designs (the size of the action space) is around 55 and slightly differs depending on the campaign. State (or context) $S_0$ is user and ad characteristics used by the algorithm. $S_0$ is high dimensional and has tens of thousands of possible values.
We utilize this data to examine the performance of our method. For each campaign and batch, we regard the MAB algorithm as the behavior policy $\pi_b$ and the CB algorithm as the evaluation policy $\pi_e$. We use the log data from users assigned to the MAB algorithm to estimate the expected value of the CB algorithm. We also compute the actual click rate using the log data from users assigned to the CB algorithm, which serves as the ground truth. We then compare the predicted value with the actual one.
We consider the relative root mean squared error (relative-RMSE) as a performance metric. Let $C$ and $B$ denote the number of campaigns ($C=4$) and the number of batches ($B=28$), respectively. Let $N_{c,b}$ denote the number of user impressions assigned to the CB algorithm for campaign $c$ in batch $b$. Let $\hat V^{\pi_e}_{c,b}$ and $\bar V^{\pi_e}_{c,b}$ be the estimated value and the actual value (the click rate) for campaign $c$ in batch $b$. We define relative-RMSE as follows:
As the actual click rate $\bar V^{\pi_e}_{c,b}$ varies across campaigns and batches, we normalize the prediction error $\hat V^{\pi_e}_{c,b}-\bar V^{\pi_e}_{c,b}$ by dividing it by the actual value $\bar V^{\pi_e}_{c,b}$ to equally weight every campaign-batch combination.
We examine the performance of four estimators, IPW Strehl2010, full-data DR, half-data DR, and DML (our proposal). We use LightGBM ke2017lightgbm to estimate the mean reward function $\mu$. First, we split the data used for estimating $\mu$ into training and validation sets, train the reward model with the training set, and tune hyperparameters with the validation set. We then use the tuned model and the original data to obtain the reward estimator. We estimate the behavior policy $\pi_b$ by taking the empirical shares of ad designs in the log data. We also try using the true behavior policy, computed by the Monte Carlo simulation from the beta distribution used in Thompson Sampling the MAB algorithm uses. For DML, we use $K = 2$ folds to perform cross-fitting. We verify that increasing $K$ to 3 or 4 changes the results little.
We present our key empirical result in Table (ref). Regardless of whether the true behavior policy is available or not, DML outperforms IPW by more than 20 % and outperforms both of the two variants of DR by more than 10 % in terms of relative-RMSE. These differences are statistically significant at 5% level. This result supports the practical value of our proposed estimator.
Finally, we use our method to measure the performance of a new counterfactual policy of choosing the ad design predicted to maximize the CTR. We obtain the counterfactual policy by training a click prediction model with LightGBM on the data from the first $B-1$ batches. In the training, we set the number of leaves to $20$ and the learning rate to $0.01$, and decide the number of boost rounds by cross-validation. We then use our proposed estimator to predict the performance of the counterfactual policy on the data from the last batch ($b=B$). The last batch contains three campaigns. The resulting off-policy evaluation results are reported in Figure (ref). The results show that the counterfactual policy performs better than the existing algorithms in two of the three campaigns, with a statistically significant improvement in one of the campaigns.
In this section, we discuss the relationships between several most closely related papers and ours. Our methodological building block is the Chernozhukov2018's DML framework for estimation and inference of semiparametric models. Our paper is among the first to show the applicability of their framework to the problem of off-policy evaluation for RL and bandits. In addition, we validate the practical importance of the two features of DML, cross-fitting and Neyman orthogonality, by both simulation and real world experiments.
Within the off-policy evaluation literature, our paper mostly builds on the work by Jiang16 and thomas2016data, who propose DR estimators. Jiang16 are the first to apply the DR technique to off-policy evaluation for RL. thomas2016data extend Jiang16's method to the settings where the time horizon is infinite and/or unknown, and propose a new estimator which combines a model-based estimator with a DR estimator. As explained in Section 3, the DR estimators of these two papers and our DML are based on the same expression of $V^{\pi_e}$ (Eq. ((ref))). The difference between our estimator and theirs is that we use the cross-fitting procedure while they either use the full sample both for estimating $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ and for computing the final estimate of $V^{\pi_e}$, or use a half of the data to construct estimates of $\pi_b$ and $\{ q_t^{\pi_e}\}_{t=0}^T$ and the other half to compute the estimate of $V^{\pi_e}$. In their experiments, Jiang16 also implement cross-fitting as one variant of the DR estimator, which is the same as our DML estimator if we use the true behavior policy $\pi_b$. However, we (i) allow the behavior policy to be unknown, (ii) allow for having a large number of potentially important state variables, and (iii) present statistical properties of the estimator while they do not.
The most closely related paper to ours is kallus2019DRL, who also apply DML to off-policy evaluation for RL. Their work and our work are independent and simultaneous. More substantively, there are several key differences between their paper and ours. Empirically, their paper does not have an application to a real dataset. In contrast, we show the estimator's practical performance in a real product setting. In addition, we provide a consistent estimator of the asymptotic variance, and allow for having a large number of potentially important state variables while they do not. Methodologically, their estimator differs from ours in that they estimate the marginal distribution ratio $\frac{\Pr_{H\sim \pi_e}(S_t,A_t)}{\Pr_{H\sim \pi_b}(S_t,A_t)}$ and plug it in place of $\rho_t^{\pi_e}(H_t^{s,a};\pi_b)$ of our estimator. There are tradeoffs between the two approaches. Their estimator is better in that it has a smaller asymptotic variance. In the important case with a known behavior policy, however, our estimator has the following advantages. (1) Our approach only needs to estimate $\{q_t^{\pi_e}\}_{t=0}^T$, while their approach needs to estimate the marginal distribution ratio as well as $\{q_t^{\pi_e}\}_{t=0}^T$. (2) Our estimator is $\sqrt{N}$-consistent and can be asymptotic normal even when the nonparametric estimator of $\{q_t^{\pi_e}\}_{t=0}^T$ does not converge to the true one, or when the parametric model of $\{q_t^{\pi_e}\}_{t=0}^T$ is misspecified. In contrast, their estimator may not be even consistent in such a case unless the nonparametric estimator of the marginal distribution ratio is consistent, or the parametric model of the marginal distribution ratio is correctly specified.
This paper proposes a new off-policy evaluation method, by marrying the DR estimator with Double/Debiased Machine Learning. Our estimator has two features. First, unlike the IPW estimator, it is robust to the bias in the estimates of the behavior policy and of the state-action value function (Neyman orthogonality). Second, we use a sample-splitting procedure called cross-fitting. This removes the overfitting bias that would arise without sample splitting but still makes full use of data, which makes our estimator better than DR estimators. Theoretically, we show that our estimator is $\sqrt{N}$-consistent and asymptotically normal with a consistent variance estimator, thus allowing for correct statistical inference. Our experiments show that our estimator outperforms the standard DR, IPW and DM estimators in terms of the mean squared error. This result not only demonstrates the capability of our estimator to reduce prediction errors, but also suggests the more general possibility that the two features of DML (Neyman orthogonality and cross-fitting) may improve many variants of the DR estimator such as MAGIC Thomas16, SWITCH wang2016optimal and MRDR Farajtabar2018MoreRD.