Doubly robust off-policy evaluation with shrinkage

Yi Su, Maria Dimakopoulou, Akshay Krishnamurthy, Miroslav Dudík

Introduction

Many real-world applications, ranging from online news recommendation (Li et al., 2011), advertising (Bottou et al., 2013), and search engines (Li et al., 2015) to personalized healthcare (Zhou et al., 2017), are naturally modeled by the contextual bandit protocol (Langford & Zhang, 2008), where a learner repeatedly observes a context, takes an action, and accrues reward. In news recommendation, the context is any information about the user, such as history of past visits, the action is the recommended article, and the reward could indicate the user’s click on the article. The goal is to maximize the reward, but the learner can only observe the reward for chosen actions, and not for the others.

We study a fundamental problem in contextual bandits known as off-policy evaluation, where the goal is to use the data gathered by a past algorithm, known as the logging policy, to estimate the average reward of a new algorithm, known as the target policy. High-quality off-policy estimates help avoid costly A/B testing and can also be used as subroutines for optimizing a policy (Dudík et al., 2011).

The most accurate approaches to off-policy evaluation are variants of doubly robust (DR) estimators (Robins & Rotnitzky, 1995; Bang & Robins, 2005; Dudík et al., 2011). DR estimation begins by fitting a regression model to predict rewards as a function of context and action. The fitted model can be used to impute unobserved rewards of the target policy on the training data, but such a direct estimate is typically biased. Instead, DR adds a correction term obtained by importance weighting the difference between observed rewards and predicted rewards. The resulting approach is unbiased, and it is asymptotically optimal under weaker assumptions than other methods (Rothe, 2016). However, its finite-sample variance can still be quite high when importance weights (also known as inverse propensity scores) are large. Therefore, several works have developed variants of DR that clip or remove large importance weights. Although weight clipping incurs some bias, it substantially decreases the variance and can yield a lower mean squared error (Bembom & van der Laan, 2008; Bottou et al., 2013; Wang et al., 2017; Su et al., 2018). These works motivate weight shrinkage as a heuristic for trading off bias and variance, but they do not provide insight into when and how these different methods should be used.

In this paper, we ask: What are the systematic strategies for shrinking importance weights? We seek to answer this question without making strong assumptions about the quality of the reward predictor, but we would like to adapt to its quality. We make the following contributions:

We derive a general framework for shrinking the importance weights by optimizing a sharp bound on the mean squared error (MSE). We use two bounding techniques. The first is agnostic to the quality of the reward estimator and yields pessimistic shrinkage estimators. The second incorporates the quality of the reward predictor and yields optimistic shrinkage estimators.

We provide theoretical justification for the standard practice of weight clipping by showing that it corresponds to pessimistic shrinkage.

Using optimistic shrinkage, we derive new estimators, which are also applicable to combinatorial actions, arising, for example, when a news portal is recommending not just a single article, but a list of articles (Cesa-Bianchi & Lugosi, 2012; Swaminathan et al., 2017).

Apart from the conceptual and theoretical contributions above, we also carry out an extensive empirical evaluation. For atomic (i.e., non-combinatorial) actions, we consider 108 experimental conditions derived from 9 real-world data sets and covering a range of data set sizes, feature dimensions, policy overlap (i.e., the magnitude of importance weights), and quality of reward estimators. For combinatorial actions, we consider a standard learning-to-rank data set and vary the quality of reward estimators. In all instances, we demonstrate the efficacy of our shrinkage approach. Via extensive ablation studies, we also identify a robust configuration of our shrinkage approach that we recommend as a practical choice.

Comparison with related work. Off-policy estimation is studied in observational settings under the name average treatment effect (ATE) estimation, with many results on asymptotically optimal estimators (Hahn, 1998; Hirano et al., 2003; Imbens et al., 2007; Rothe, 2016), but only few that optimize MSE in finite samples. Most notably, Kallus (2017, 2018) develops the kernel optimal matching (KOM) approach that adjusts importance weights by optimizing MSE under smoothness (or parametric) assumptions on the reward function. This method is reminiscent of direct modeling, whose bias can be bounded under smoothness assumptions, but whose performance deteriorates if these assumptions are violated. In contrast, we optimize importance weights with essentially no modeling assumptions. Another difference is that KOM runs in time that is super-linear in the data set size, which prevents its use with large data sets, whereas our approach requires a single pass through the data and readily applies to large-scale scenarios.

Several recent works study how to improve DR estimators under similar assumptions as we make here (Wang et al., 2017; Farajtabar et al., 2018; Su et al., 2018), focusing either on weight shrinkage or on training of the reward predictor. However, to our knowledge, we are the first to provide a detailed theoretical and empirical investigation of the interplay between these two design components. For example, in Table 3, we show that the more robust doubly robust (MRDR) approach for training of the reward predictor (Farajtabar et al., 2018) performs poorly in combination with weight shrinkage. More generally, different estimators may require different reward predictors. This specific finding has practical implications that are missing in prior work.

Setup

In the off-policy evaluation problem, we are given a dataset {(xi,ai,ri)}i=1n∼μ\{(x_{i},a_{i},r_{i})\}_{i=1}^{n}\sim\mu consisting of context-action-reward triples collected by some logging policy μ\mu, and we would like to estimate the value of a target policy π\pi. The quality of an estimator V^(π)\hat{V}(\pi) is measured by the mean squared error

where the expectation is with respect to the data generation process. In analyzing the error of an estimator, we rely on the decomposition of MSE into the bias and variance terms:

We consider three standard approaches for off-policy evaluation. The first two are direct modeling (DM) and inverse propensity scoring (IPS). In DM, we train a reward predictor η^:X×A→\hat{\eta}:\mathcal{X}\times\mathcal{A}\to and use it to impute rewards. In IPS, we simply reweight the data. The two estimators are:

Let w(x,a)≔π(a∣x)/μ(a∣x)w(x,a)\coloneqq\pi(a\mathbin{|}x)/\mu(a\mathbin{|}x) denote the importance weight. We make a standard assumption that π\pi is absolutely continuous with respect to μ\mu, meaning that μ(a∣x)>0\mu(a\mathbin{|}x)>0 whenever π(a∣x)>0\pi(a\mathbin{|}x)>0. This ensures that the importance weights are well defined and V^IPS(π)\smash{\hat{V}_{\textup{IPS}}(\pi)} is an unbiased estimator of V(π)V(\pi). If there is a substantial mismatch between π\pi and μ\mu, then the importance weights will be large and V^IPS(π)\smash{\hat{V}_{\textup{IPS}}(\pi)} will have large variance. On the other hand, given any fixed reward predictor η^\hat{\eta} (fit on a separate dataset), V^DM(π)\smash{\hat{V}_{\textup{DM}}(\pi)} has low variance, but it can be biased due to approximation errors in fitting η^\hat{\eta}.

The third approach, called the doubly robust (DR) estimator, combines DM and IPS:

The DR estimator applies IPS to a shifted reward, using η^\hat{\eta} as a control variate to decrease the variance of IPS, while preserving its unbiasedness. DR is asymptotically optimal, as long as it is possible to derive sufficiently good reward predictors η^\hat{\eta} given enough data (Rothe, 2016).

However, even when the reward predictor η^\hat{\eta} is perfect, stochasticity in the rewards may cause the terms ri−η^(xi,ai)r_{i}-\hat{\eta}(x_{i},a_{i}), appearing in the DR estimator, to be far from zero. Multiplied by large importance weights w(xi,ai)w(x_{i},a_{i}), these terms yield large variance for DR in comparison with DM. As mentioned in Section 1, several approaches seek a more favorable bias–variance trade-off by shrinking the importance weights. Our work also seeks to systematically replace the weights w(xi,ai)w(x_{i},a_{i}) with new weights w^(xi,ai)\hat{w}(x_{i},a_{i}) to bring the variance of DR closer to that of DM.

where F\mathcal{F} is some function class of reward predictors. Natural choices of the weighting function zz, explored in our experiments, include z(x,a)=1z(x,a)=1, z(x,a)=w(x,a)z(x,a)=w(x,a) and z(x,a)=w2(x,a)z(x,a)=w^{2}(x,a). We stress that the assumption on how we fit η^\hat{\eta} only serves to guide our derivations, but we make no specific assumptions about its quality. In particular, we do not assume that F\mathcal{F} contains a good approximation of η\eta.

Our Approach: DR with Shrinkage

We assume that 0≤w^≤w0\leq\hat{w}\leq w, justifying the terminology “shrinkage”. For a fixed choice of π\pi and η^\hat{\eta}, we will seek the mapping w^\hat{w} that minimizes the MSE of V^DRs(π;η^,w^)\hat{V}_{\textup{DRs}}(\pi;\hat{\eta},\hat{w}), which we simply denote as MSE(w^)\textup{MSE}(\hat{w}). We similarly write Bias(w^)\textup{Bias}(\hat{w}) and Var⁡(w^)\operatorname*{{\rm Var}}(\hat{w}) for the bias and variance of this estimator.

We treat w^\hat{w} as the optimization variable and consider two upper bounds on MSE: an optimistic one and a pessimistic one. In both cases, we separately bound Bias(w^)\textup{Bias}(\hat{w}) and Var⁡(w^)\operatorname*{{\rm Var}}(\hat{w}). To bound the bias, we use the following expression, derived from the fact that V^DRs\smash{\hat{V}_{\textup{DRs}}} is unbiased when w^=w\hat{w}=w:

To bound the variance, we rely on the following proposition, which states that it suffices to focus on the second moment of the terms \hat{w}(x_{i},a_{i})\bigl{(}r_{i}-\hat{\eta}(x_{i},a_{i})\bigr{)}:

See appendix for the proof of Proposition 1 (as well as other mathematical statements from this paper).

We derive estimators for two different regimes depending on the quality of the reward predictor η^\hat{\eta}. Since we do not know the quality of η^\hat{\eta} a priori, in Section 5 we derive a model selection procedure to select between these two estimators.

Our first family of estimators is based on an optimistic MSE bound, which adapts to the quality of η^\hat{\eta}, and which we expect to be tighter when η^\hat{\eta} is more accurate. Recall that η^\hat{\eta} is trained to minimize weighted square loss with respect to some weighting function zz, which we denote as

The loss L(η^)L(\hat{\eta}) quantifies the quality of η^\hat{\eta}. We use it to bound the bias by applying the Cauchy–Schwarz inequality to (4):

where the first inequality follows by the Cauchy-Schwarz inequality, and the second from the fact that w^2(x,a)≤w2(x,a)\hat{w}^{2}(x,a)\leq w^{2}(x,a) and ∣r−η^(x,a)∣≤1\left\lvert r-\hat{\eta}(x,a)\right\rvert\leq 1.

Combining the bounds (5) and (6) with Proposition 1 yields the following bound on MSE(w^)\textup{MSE}(\hat{w}):

A direct minimization of this bound appears to be a high dimensional optimization problem. Instead of minimizing the bound directly, we note that it is a strictly increasing function of the two expectations that appear in it. Thus, its minimizer must be on the Pareto front with respect to the two expectations, meaning that for some choice of λ∈[0,∞]\lambda\in[0,\infty], it can be obtained by minimizing

with respect to w^\hat{w}. This objective decomposes across contexts and actions. Taking the derivative with respect to w^(x,a)\hat{w}(x,a) and setting it to zero yields the solution

where “o” above is a mnemonic for optimistic shrinkage. We refer to the DRs estimator with w^=w^o,λ\hat{w}=\hat{w}_{\textup{o},\lambda} as the doubly robust estimator with optimistic shrinkage (DRos) and denote it by V^DRos(π;η^,λ)\smash{\hat{V}_{\textup{DRos}}(\pi;\hat{\eta},\lambda)}. Note that this estimator does not depend on zz, although it was included in the optimization objective. When λ=0\lambda=0, we have w^(x,a)=0\hat{w}(x,a)=0 corresponding to DM. As λ→∞\lambda\to\infty, the weights increase and in the limit become equal to w(x,a)w(x,a), corresponding to standard DR.

2 DR with Pessimistic Shrinkage

Our second estimator family makes no assumptions on the quality of η^\hat{\eta} beyond the range bound η^(x,a)∈\hat{\eta}(x,a)\in, which implies ∣η^(x,a)−r∣≤1\left\lvert\hat{\eta}(x,a)-r\right\rvert\leq 1 and yields the bounds

As before, we do not optimize the resulting MSE bound directly and instead solve for the Pareto front points parameterized by λ∈[0,∞]\lambda\in[0,\infty] (we scale λ\lambda by a factor of two to obtain the solution that more cleanly matches the clipping estimator):

The objective again decomposes across context-action pairs, yielding the solution

which recovers (and justifies) existing weight-clipping approaches (Kang et al., 2007; Strehl et al., 2010; Su et al., 2018) (see Appendix A for detailed calculations). We refer to the resulting estimator as V^DRps(π;η^,λ)\smash{\hat{V}_{\textup{DRps}}(\pi;\hat{\eta},\lambda)}, for doubly robust with pessimistic shrinkage. Similarly to optimistic shrinkage, we recover DM for λ=0\lambda=0, and DR as λ→∞\lambda\to\infty.

Shrinkage for Combinatorial Actions

We showcase the generality of our optimization-based approach by deriving a shrinkage estimator for combinatorial actions (also called slates), which arise, for example, when recommending a ranked list of items.

The detailed derivation is in Appendix B. To our knowledge this is the first weight-shrinkage estimator for contextual combinatorial bandits.

Assume that μ(⋅∣x)\mu(\cdot\mathbin{|}x) is supported on a linearly independent set of actions for every xx. If w^(x)⊤a=c(x,a)w(x)⊤a\hat{\mathbf{w}}(x)^{\top}\mathbf{a}=c(x,\mathbf{a})\mathbf{w}(x)^{\top}\mathbf{a} for some c(x,a)∈c(x,\mathbf{a})\in, then

Note that the quantity ∥vπ,x∥1\lVert\mathbf{v}_{\pi,x}\rVert_{1} on the right-hand side only depends on the set Bx\mathcal{B}_{x}, but not on the probabilities with which μ\mu chooses a∈Bx\mathbf{a}\in\mathcal{B}_{x}. Non-combinatorial setting of Section 3 is a special case of the linearly independent setting, where d=∣A∣d=\lvert\mathcal{A}\rvert and actions are represented by standard basis vectors. In this case, ∥vπ,x∥1=1\lVert\mathbf{v}_{\pi,x}\rVert_{1}=1 and we recover Proposition 1. We can always select Bx\mathcal{B}_{x} to be an (approximate) barycentric spanner and achieve ∥vπ,x∥1=O(d)\lVert\mathbf{v}_{\pi,x}\rVert_{1}=O(d) (Awerbuch & Kleinberg, 2008; Dani et al., 2008).

Model Selection

All of our shrinkage estimators have hyperparameters which we condense into a tuple θ\theta. For example θ=(η^,o,λ)\theta=(\hat{\eta},\textup{o},\lambda) denotes that we are using a reward predictor η^\hat{\eta} and optimistic shrinkage with the parameter λ\lambda. To select among these hyperparameters, we propose and analyze a simple model selection procedure.

Let V^θ\hat{V}_{\theta} denote the estimator parameterized by θ\theta. We consider the procedure that estimates the variance of V^θ\hat{V}_{\theta} by sample variance Var⁡^(θ)\smash{\widehat{\operatorname*{{\rm Var}}}(\theta)}, and bounds the bias of V^θ\hat{V}_{\theta} by a data-dependent upper bound BiasUB(θ)\textup{BiasUB}(\theta). The only requirement is that for all θ\theta, Bias(θ)≤BiasUB(θ)\textup{Bias}(\theta)\leq\textup{BiasUB}(\theta) (with high probability), and that BiasUB(θ)=0\textup{BiasUB}(\theta)=0 whenever Bias(θ)=0\textup{Bias}(\theta)=0; this holds for both bias bounds from Section 3, as they become zero when w^=w\hat{w}=w. Now, to choose θ\theta from a set of hyperparameters Θ\Theta, we optimize the estimate of the MSE:

The next theorem shows that this procedure always compares favorably with all the unbiased estimators included in Θ\Theta, up to an asymptotically negligible term O(n−3/2)O(n^{-3/2}). In particular, the procedure is asymptotically optimal whenever Θ\Theta includes a standard (non-shrunk) DR.

Let Θ\Theta be a finite set of hyperparameter values and let Θ0≔{θ∈Θ: Bias(θ)=0}\Theta_{0}\coloneqq\{\theta\in\Theta:\>\textup{Bias}(\theta)=0\} denote the subset of unbiased estimators. Assume that with probability 1−δ/21-\delta/2 we have Bias(θ)≤BiasUB(θ)\textup{Bias}(\theta)\leq\textup{BiasUB}(\theta) for all θ∈Θ\theta\in\Theta. Then there exists a universal constant CC such that with probability at least 1−δ1-\delta we have

There are many strategies to construct data-dependent bias bounds with the required properties. The three bounds in our experiments take form of sample averages that approximate expectations in: (i) the expression for the bias given in (4), (ii) the optimistic bias bound in (5), and (iii) the pessimistic bias bound in (7). In our theory, these estimates need to be adjusted to obtain high-probability confidence bounds. In our experiments, we evaluate both the basic estimates and adjusted variants where we add twice the standard error.

Our model selection procedure is related to MAGIC (Thomas & Brunskill, 2016) as well as the procedure for the SWITCH estimator (Wang et al., 2017). Unlike MAGIC, we pick a single hyperparameter value θ\theta rather than aggregating several, and we use different bias and variance estimates. SWITCH uses our pessimistic bias bound (7), but with no theoretical justification. We use two additional bounding strategies, which are empirically shown to help, and provide theoretical justification in the form of an oracle inequality.

Experiments

We evaluate our new estimators on the tasks of off-policy evaluation and off-policy learning and compare their performance with previous estimators. Our secondary goal is to identify the configuration of the shrinkage estimator that is most robust for use in practice.

Datasets. Following prior work (Dudík et al., 2014; Wang et al., 2017; Farajtabar et al., 2018; Su et al., 2018), we simulate bandit feedback on 9 UCI multi-class classification datasets. This lets us evaluate estimators in a broad range of conditions and gives us ground-truth policy values (see Table 4 in the appendix for the dataset statistics). Each multi-class dataset with kk classes corresponds to a contextual bandit problem with kk possible actions coinciding with classes. We consider either deterministic rewards, where on multiclass example (x,y∗)(x,y^{*}), the action yy yields the reward r=1{y=y∗}r={\bf 1}\{y=y^{*}\}, or stochastic rewards where r=1{y=y∗}r={\bf 1}\{y=y^{*}\} with probability 0.750.75 and r=1−1{y=y∗}r=1-{\bf 1}\{y=y^{*}\} otherwise. For every dataset, we hold out 25%25\% of the examples to measure ground truth. On the remaining 75%75\% of the dataset, we use logging policy μ\mu to simulate nn bandit examples by sampling a context xx from the dataset, sampling an action y∼μ(⋅∣x)y\sim\mu(\cdot\mathbin{|}x) and then observing a deterministic or stochastic reward rr. The value of nn varies across experimental conditions.

Policies. We use the 25%25\% held-out data to obtain logging and target policies as follows. We first obtain two deterministic policies π1,det\pi_{1,\textrm{det}} and π2,det\pi_{2,\textrm{det}} by training two logistic models on the same data, but using either the first or second half of the features. We obtain stochastic policies parameterized by (α,β)(\alpha,\beta), following the softening technique of Farajtabar et al. (2018). Specifically, π1,(α,β)(a∣x)=(α+βu)\pi_{1,(\alpha,\beta)}(a\mathbin{|}x)=(\alpha+\beta u) if a=π1,det(x)a=\pi_{1,\textrm{det}}(x) and π1,(α,β)(a∣x)=1−α−βuk−1\pi_{1,(\alpha,\beta)}(a\mathbin{|}x)=\frac{1-\alpha-\beta u}{k-1} otherwise, where u∼Unif([−0.5,0.5])u\sim\textrm{Unif}([-0.5,0.5]). In off-policy evaluation experiments, we consider a fixed target and several choices of logging policy (see Table 3). In off-policy learning we use π1,(0.9,0)\pi_{1,(0.9,0)} as the logging policy.

Baselines. We include a number of estimators in our evaluation: the direct modeling approach (DM), doubly-robust approach (DR) and its self-normalized variant (snDR), our approach (DRs), and the doubly-robust version of the switch estimator of Wang et al. (2017), which also performs a form of weight clipping.For simplicity we call this estimator switch, although Wang et al. call it switch-DR. Note that DR with η^≡0\hat{\eta}\equiv 0 is identical to inverse propensity scoring (IPS); we refer to its self-normalized variant as snIPS. Our estimator and switch have hyperparameters, which are selected by their respective model selection procedures (see Appendix D for details about the hyperparameter grid).

We begin by evaluating different configurations of DRs via an ablation analysis. Then we compare DRs with baseline estimators. We have a total of 108108 experimental conditions: for each of the 99 datasets we use 66 logging policies and consider stochastic or deterministic rewards. Except for the learning curves below, we always take nn to be all available bandit data (75%75\% of the overall dataset).

Ablation analysis. We conduct two ablation studies: one evaluating different reward predictors and the other evaluating the optimistic and pessimistic shrinkage types.

In Table 3, for each fixed estimator type (e.g., DR) we evaluate each reward predictor by reporting the number of conditions where it is statistically indistinguishable from the best and the number of conditions where it statistically dominates all other predictors. For DRs we use oracle tuning for the shrinkage type and coefficient λ\lambda. The table shows that weight shrinkage strongly influences the choice of regressor. For example, z≡1z\equiv 1 and z=wz=w are top choices for DR, but with the inclusion of shrinkage in DRs, z=w2z=w^{2} emerges as the best choice. In our comparison experiments below, we run each method with its best reward predictor: DM with z≡1z\equiv 1, snDR with z=wz=w, and DRs and switch with z=w2z=w^{2}. For DRs and switch, we additionally also consider η^≡0\hat{\eta}\equiv 0, because it allows including IPS as their special case. Somewhat surprisingly, in our experiments, MRDR is dominated by other reward predictors (except for η^≡0\hat{\eta}\equiv 0), and this remains true even with a deterministic target policy (see Table 5 in the appendix).

In Table 3, we compare optimistic and pessimistic shrinkage when paired with a fixed reward predictor (using oracle tuning for λ\lambda). We report how many times each estimator statistically dominates the other. The results suggest that both shrinkage types are important for robust performance across conditions, so we consider both choices going forward.

Comparisons. In Figure 1 (left two plots), we compare our new estimator with the baselines. We visualize the results by plotting the cumulative distribution function (CDF) of the normalized MSE of each method (normalized by the MSE of snIPS) across the experimental conditions. Better performance corresponds to CDF curves towards the top-left corner, meaning the method achieves a lower MSE more frequently. The first plot summarizes 54 conditions where the reward is deterministic, while the second plot considers the 54 stochastic reward conditions. For DRs we consider two model selection procedures outlined in Section 5 that differ in their choice of BiasUB. DRs-direct estimates the expectations in the expressions in Eqs. (4), (5), and (7) (corresponding to the bias and bias bounds) by empirical averages and takes their pointwise minimum. DRs-upper adds to these estimates twice their standard error, before taking minimum, more closely matching our theory. For DRs, we use the zero reward predictor and the one trained with z=w2z=w^{2}, and we always select between both shrinkage types. Since switch also comes with a model selection procedure, we use it to select between the same two reward predictors as DRs.

In the deterministic case (the first plot), we see that DRs-upper has the best aggregate performance, by a large margin. DRs-direct also has better aggregate performance than the baselines on most of the conditions. In the stochastic case (the second plot), DRs-direct has similarly strong performance, but DRs-upper degrades considerably, suggesting this model selection scheme is less robust to stochastic rewards. We illustrate this phenomenon in the right two plots of Figure 1, plotting the MSE as a function of the number of samples for one choice of a logging policy and dataset, first with deterministic rewards and then with stochastic rewards. Because of a more robust performance, we therefore advocate for DRs-direct as our final method.

1.2 Off-policy Learning

In Figure 3, we show the performance of four methods (DM, DR, IPS, and DRs-direct) on four of the UCI datasets. For each method, we compute the average value of the learned policy on the test set (averaged over 10 replicates) and report this value normalized by the average value for IPS. For DM and DR, we select the hyperparameter γ\gamma and reward predictor optimally in hindsight, while for DRs we use our model selection. Note that we do not compare with switch here as it is not amenable to gradient-based optimization (Su et al., 2018). We find that off-policy learning using DRs-direct always outperforms the baselines, with the exception of the optdigits dataset, where all the methods perform similarly.

2 Combinatorial Setting

Reward predictors. We consider two reward predictors η^\hat{\eta} trained on logged data. Both are trained via ridge regression, but differ in feature sets they consider: ridge(all) is trained on all features, ridge(5) is trained on the five features that are most correlated with the reward.

Baselines. We compare our method (DRs-PI) with DM and DR-PI.DR-PI dominates self-normalized version of DR-PI as well as the standard pseudo-inverse estimator (i.e., with η^≡0\hat{\eta}\equiv 0). In DRs-PI we select the hyperparameter λ\lambda from a geometrically spaced grid using our model selection procedure with the empirical version of Eq. (4) in place of bias bound and also consider the oracle tuning of λ\lambda from the same grid (details in Appendix D.2).

Results and discussion. In Figure 2 we show the MSE of all the methods as a function of sample size, averaged over 20 replicates. Across all conditions, DRs-PI outperforms DR by a factor of 1.5 or more (note that MSE is reported on log scale). A more striking result is the superior quality of the oracle-tuned DRs-PI. It shows that the shrinkage strategy is highly effective in achieving a good bias–variance trade-off, but to unlock its potential in combinatorial settings requires improvements in model selection.

Conclusion

In this paper, we have derived shrinkage-based doubly-robust estimators for off-policy evaluation using a principled optimization-based framework. Our approach recovers the weight-clipping estimator from prior work and also yields novel optimistic shrinkage estimators for both atomic and combinatorial settings. Extensive experiments demonstrate the efficacy of these estimators and highlight the role of model selection in achieving good performance. Thus, the next step is to develop model selection procedures for off-policy evaluation that can close the gap with oracle tuning. We look forward to pursuing this direction in future work.

Acknowledgements

This work was partially completed during Yi’s and Maria’s internships at Microsoft Research. Yi is also supported by the Bloomberg Data Science Fellowship. All content represents the opinion of the authors, which is not necessarily shared or endorsed by their respective employers or sponsors.

References

Appendix A Derivation of Shrinkage Estimators for Non-combinatorial Setting

In this section we provide detailed derivations for the two estimators in the non-combinatorial setting.

We first derive the pessimistic version. Recall that the optimization problem decouples across (x,a)(x,a), so we focus on a single (x,a)(x,a) pair such that μ(a∣x)>0\mu(a\mathbin{|}x)>0 since only such pairs can appear in the data. For conciseness, we omit the dependence on (x,a)(x,a) and simply write w=w(x,a)w=w(x,a), w^=w^(x,a)\hat{w}=\hat{w}(x,a) and μ=μ(a∣x)\mu=\mu(a\mathbin{|}x). Fixing λ≥0\lambda\geq 0, we must solve

Since μ>0\mu>0, the first equation can be rewritten as

Now a simple case analysis shows that if w>λw>\lambda then the choice w^=λ\hat{w}=\lambda, v=−1v=-1 satisfies Eq. (11), and if 0≤w≤λ0\leq w\leq\lambda then the choice w^=w\hat{w}=w, v=−w/λv=-w/\lambda satisfies Eq. (11), yielding

For the optimistic version, the optimization problem is

where z=z(x,a)z=z(x,a). The optimality conditions are

Notice that this estimator does not depend on the weighting function zz, so it does not depend on how we train the regression model.

Appendix B Derivation of the Shrinkage Estimator for Combinatorial Setting

We provide a complete derivation of the shrinkage estimator for combinatorial actions. We use the notation w(x,a)=w(x)⊤aw(x,\mathbf{a})=\mathbf{w}(x)^{\top}\mathbf{a}. We assume that the regression model takes form η^(x,a)=η^(x)⊤a\hat{\eta}(x,\mathbf{a})=\hat{\bm{\eta}}(x)^{\top}\mathbf{a}, satisfies η^(x,a)∈\hat{\eta}(x,\mathbf{a})\in, and is trained to minimize

for some z(x,a)>0z(x,\mathbf{a})>0. We assume that the linearity assumption holds, so we can write η(x,a)=η(x)⊤a\eta(x,\mathbf{a})=\bm{\eta}(x)^{\top}\mathbf{a}. And we also assume that \operatorname{span}\bigl{(}\operatorname{supp}\pi(\cdot\mathbin{|}x)\bigr{)}\subseteq\operatorname{span}\bigl{(}\operatorname{supp}\mu(\cdot\mathbin{|}x)\bigr{)}, so, as shown by Swaminathan et al. (2017), the pseudo-inverse estimator is unbiased:

Therefore, if we replace ww by an arbitrary function w^\hat{w} (not necessarily linear in a\mathbf{a}), we obtain the expression for the bias

Using Cauchy–Schwarz inequality, we can bound the bias in terms of L(η^)L(\hat{\bm{\eta}}):

For the variance bound, we begin with a proxy based on Proposition 2, and then bound it using Cauchy–Schwarz inequality, the fact that η^(x)⊤a\hat{\bm{\eta}}(x)^{\top}\mathbf{a} and rr are bounded in $,andanadditionalassumptionthat, and an additional assumption that\left\lvert\hat{w}(x,\mathbf{a})\right\rvert\leq\left\lvert w(x,\mathbf{a})\right\rvert$ (which we will show is true for the specific optimistic estimator that we derive below):

Similar to non-combinatorial setting, the solutions of the resulting MSE bound must lie on the Pareto front parameterized by a single scalar λ∈[0,∞]\lambda\in[0,\infty]:

This decomposes across (x,a)(x,\mathbf{a}) and by first-order optimality, we obtain the same solution as in non-combinatorial setting:

Appendix C Proofs

For T1T_{1}, since ∑a′∈Aπ(a′∣x)η^(x,a′)\sum_{a^{\prime}\in\mathcal{A}}\pi(a^{\prime}\mathbin{|}x)\hat{\eta}(x,a^{\prime}) does not depend on a,ra,r, it does not contribute to the conditional variance, and we get

The first term is our variance proxy. To bound the second term, write w^(x,a)=c(x,a)w(x,a)\hat{w}(x,a)=c(x,a)w(x,a) for some c(x,a)∈c(x,a)\in, which is possible since 0≤w≤w^0\leq w\leq\hat{w} by assumption. The second term can then be rewritten and bounded as

where the second equality follows by the unbiasedness of inverse-propensity scoring, and the final bound follows because c(x,a),r,η^(x,a)∈c(x,a),r,\hat{\eta}(x,a)\in.

For T2T_{2}, of course we have T2≥0T_{2}\geq 0, and further

where we again write w^(x,a)=c(x,a)w(x,a)\hat{w}(x,a)=c(x,a)w(x,a) for some c(x,a)∈c(x,a)\in, then appeal to the unbiasedness of the inverse-propensity scoring, and finally use the bounds c(x,a),η(x,a),η^(x,a)∈c(x,a),\eta(x,a),\hat{\eta}(x,a)\in.

and we have just shown that the right hand side is in \bigl{[}-\frac{1}{n},\frac{1}{n}\bigr{]}. This proves the proposition.

C.2 Proof of Proposition 2

To obtain its psedo-inverse, we use tho following fact:

where K−1K^{-1} is well defined thanks to the linear independence of columns of BB.

Let G≔BDB⊤G\coloneqq BDB^{\top} and G′≔BK−1D−1K−1B⊤G^{\prime}\coloneqq BK^{-1}D^{-1}K^{-1}B^{\top}. To show that G†=G′G^{\dagger}=G^{\prime}, it suffices to argue that GG′G=GGG^{\prime}G=G and G′GG′=G′G^{\prime}GG^{\prime}=G^{\prime}:

We are now ready to start the proof of Proposition 2. Similarly to the proof of Proposition 1, we first apply the law of total variance

To analyze T2T_{2}, we first rewrite and bound the inner expectation for a fixed xx. We drop the dependence on xx from the notation and write η^=η^(x)\hat{\bm{\eta}}=\hat{\bm{\eta}}(x), qπ=qπ,x\mathbf{q}_{\pi}=\mathbf{q}_{\pi,x}, w^=w^(x)\hat{\mathbf{w}}=\hat{\mathbf{w}}(x), η=η(x)\bm{\eta}=\bm{\eta}(x), and Γμ=Γμ,x\Gamma_{\mu}=\Gamma_{\mu,x}:

where in the last step we introduced shorthands B=BxB=B_{x}, Dμ=Dμ,xD_{\mu}=D_{\mu,x} and vπ=vπ,x\mathbf{v}_{\pi}=\mathbf{v}_{\pi,x}.

To continue with the derivation, observe that by assumption, we have w^⊤a=c(x,a)w⊤a\hat{\mathbf{w}}^{\top}\mathbf{a}=c(x,\mathbf{a})\mathbf{w}^{\top}\mathbf{a} for all a∈Bx\mathbf{a}\in\mathcal{B}_{x}, and so we can write w^⊤B=w⊤BC\hat{\mathbf{w}}^{\top}B=\mathbf{w}^{\top}BC where CC is a diagonal matrix with entries c(x,a)c(x,\mathbf{a}) across a∈Bx\mathbf{a}\in\mathcal{B}_{x}. Next, using the fact that w=Γμ†qπ\mathbf{w}=\Gamma_{\mu}^{\dagger}\mathbf{q}_{\pi} and then plugging in Eq. (14), we obtain

where we introduced the shorthand K=KxK=K_{x}. Now combining Eqs. (15) and (16), we obtain

To bound T1T_{1}, we first note that η^(x)⊤qπ,x\hat{\bm{\eta}}(x)^{\top}\mathbf{q}_{\pi,x} is independent of a\mathbf{a} and rr, and so it does not contribute to the variance, and so

where we applied similar reasoning as in Eq. (17). Combining this bound with the bound on T2T_{2} completes the proof:

C.3 Proof of Theorem 3

The main technical part of the proof is a deviation inequality for the sample variance. For this, let us fix θ\theta, which we drop from notation, and focus on estimating the variance

Let Z1,…,ZnZ_{1},\ldots,Z_{n} be iid random variables, and assume that ∣Zi∣≤R|Z_{i}|\leq R almost surely. Then there exists a constant C>0C>0 such that for any δ∈(0,1)\delta\in(0,1), with probability at least 1−δ1-\delta

We work with the second term first. Let Z1′,…,Zn′Z_{1}^{\prime},\ldots,Z_{n}^{\prime} be an iid sample, independent of Z1,…,ZnZ_{1},\ldots,Z_{n}. Now, by Theorem 3.4.1 of De la Pena & Giné (2012), we have

for a universal constant C>0C>0. Thus, we have decoupled the U-statistic. Now let us condition on Z1,…,ZnZ_{1},\ldots,Z_{n} and write Xj=1n−1∑i≠j(Zi−μ)X_{j}=\frac{1}{n-1}\sum_{i\neq j}(Z_{i}-\mu), which conditional on Z1,…,ZnZ_{1},\ldots,Z_{n} is non-random. We will apply Bernstein’s inequality on 1n∑j=1nXj(Zj′−μ)\frac{1}{n}\sum_{j=1}^{n}X_{j}(Z^{\prime}_{j}-\mu), which is a centered random variable, conditional on Z1:nZ_{1:n}. This gives that with probability at least 1−δ1-\delta

This bound holds with high probability for any {Xj}j=1n\{X_{j}\}_{j=1}^{n}. In particular, since ∣Xj∣≤R|X_{j}|\leq R almost surely, we get that with probability 1−δ1-\delta

The factors of CC arise from working through the decoupling inequality.

Next, by a standard application of Bernstein’s inequality, with probability at least 1−δ1-\delta, we have

Therefore, with probability 1−2δ1-2\delta we have

Let us now address the first term, a simple application of Bernstein’s inequality gives that with probability at least 1−δ1-\delta

Combining the two inequalities, we obtain the result. ∎

Since we are estimating the variance of the sample average estimator, we divide by another factor of nn. Meanwhile the range and the variance terms themselves are certainly O(1)O(1), so the error terms in Lemma 5 are O(n−3/2)O(n^{-3/2}) and O(n−2)O(n^{-2}) respectively. Formally, there exists a universal constants C1,C2>0C_{1},C_{2}>0 such that for any δ∈(0,1)\delta\in(0,1) with probability at least 1−δ1-\delta we have

By adjusting the constant, we can simplify the expression by removing the n−2n^{-2} term. In other words, there exists a different universal constant C>0C>0 such that

holds with probability at least 1−δ1-\delta.

For the model selection result, first apply Lemma 5 for all θ∈Θ\theta\in\Theta, taking a union bound. Further take a union bound over the event that Bias(θ)≤BiasUB(θ)\textup{Bias}(\theta)\leq\textup{BiasUB}(\theta) for all θ∈Θ\theta\in\Theta, if it is needed. Then, observe that for any θ0∈Θ0\theta_{0}\in\Theta_{0} we have

The first inequality uses Lemma 5 and the fact that Bias≤BiasUB\textup{Bias}\leq\textup{BiasUB}. The second uses that θ^\hat{\theta} optimizes this quantity, and the third uses the property that BiasUB(θ0)=0\textup{BiasUB}(\theta_{0})=0 by assumption. Note that the universal constant here is slightly different from the one in the variance bound, since we have also taken a union bound for the bias term.

C.4 Construction of Upper Bounds on Bias

In this section we give detailed construction of bias upper bounds that we use in the model selection procedure. Recall that this is for the analysis only. Empirically we found that using the estimators alone — not the upper bounds — leads to better performance.

Throughout, we fix a set of hyperparameters θ\theta, which we suppress from the notation.

The most straightforward bias estimator is to simply approximate the expectation with a sample average.

Inflating the estimate by the right hand side gives BiasUB, which is a high probability upper bound on Bias.

The bias bound used in the pessimistic estimator and its natural sample estimator are

Note that since we have already eliminated the dependence on the reward, we can analytically evaluate the expectation over actions, which will lead to lower variance in the estimate.

Again we perform a fairly naive analysis. Since 0≤w^(x,a)≤w(x,a)0\leq\hat{w}(x,a)\leq w(x,a), the random variables, equal to the inner sum over aa, take values in $.Therefore,Hoeffding’sinequalitygivesthatwithprobability. Therefore, Hoeffding’s inequality gives that with probability1-\delta$

and we use the right hand side for our high probability upper bound.

For the optimistic bound, we must estimate two terms, one involving the regressor and one involving the importance weights. We use

Note here that the former uses sampled actions from μ\mu, but does not involve the importance weight, while the latter involves the importance weight but analytically evaluates the expectation over μ\mu. Thus we can expect that both are fairly low variance.

For T2T_{2}, we similarly to the pessimistic case convert the inner expectation w.r.t. μ(a∣xi)\mu(a\mathbin{|}x_{i}) to an expectation w.r.t. π(a∣xi)\pi(a\mathbin{|}x_{i}), obtaining a random variable bounded between 0 and max⁡x,aw(x,a)/z(x,a)\max_{x,a}w(x,a)/z(x,a). Using Hoeffding’s inequality, we obtain that with probability 1−2δ1-2\delta

The high probability upper bound follows by multiplying the two right hand sides together and taking square root.

Appendix D Experimental Details and Additional Results

We use datasets from the UCI Machine Learning Repository (Dua & Graff, 2017). Dataset statistics are displayed in Table 4.

For our shrinkage estimators and switch, we choose the shrinkage coefficients from a grid of 30 geometrically spaced values. For the pessimistic estimator and switch, the largest and smallest values in the grid are the 0.050.05 quantile and 0.950.95 quantile of the importance weights. For the optimistic estimator, the largest and smallest values are 0.01×(w0.05)20.01\times(w_{0.05})^{2} and 100×(w0.95)2100\times(w_{0.95})^{2} where w0.05w_{0.05} and w0.95w_{0.95} are the 0.05 and 0.95 quantile of the importance weights.

For the off-policy learning experiments, we only consider the shrinkage coefficients in {0.0,0.1,1,10,100,1000,∞}\{0.0,0.1,1,10,100,1000,\infty\} during training, while for model selection, we use the same grid as in the evaluation experiments.

Farajtabar et al. (2018) propose training the regression model with a specific choice of weighting zz, which we also use in our experiments. When the evaluation policy π\pi is deterministic, they set z(x,a)=1{π(x)=a}⋅1−μ(a∣x)μ(a∣x)2z(x,a)={\bf 1}\{\pi(x)=a\}\cdot\frac{1-\mu(a\mathbin{|}x)}{\mu(a\mathbin{|}x)^{2}}. For stochastic policies, following the implementation of Farajtabar et al., we sample ai∼π(⋅∣xi)a_{i}\sim\pi(\cdot\mathbin{|}x_{i}) for each example in the dataset used to train the reward predictor. Then we proceed as if the evaluation policy deterministically chooses aia_{i} on example xix_{i}.

Since MRDR is more suited to deterministic policies, we also report the results of our regressor and shrinkage ablations for a deterministic target policy π1,det\pi_{1,\textrm{det}} in Table 5. As with the stochastic policies, the estimator influences the choice of reward predictor, but note that z=1z=1 and z=wz=w are more favorable here. This is likely due to high variance suffered from training with z=w2z=w^{2}, because the importance weights are larger with a deterministic policy. Our shrinkage ablation reveals that both estimator types are important also when the target policy is deterministic.

In Table 6, we show the comparison of different model selection methods under different reward predictors and different shrinkage types. In most cases, Dir-all (DRs-direct where the bias bound is estimated as the pointwise minimum of (1) the bias, (2) the optimistic bound and (3) the pessimistic bound) and Up-all (all bias estimates are adjusted by adding twice standard error before taking pointwise minimum) are most frequently statistically indistinguishable from the best, which suggests that our proposed bias estimate (by taking pointwise minimum of the three) is robust and adaptive.

In Figure 4 and Figure 5, we compare our new estimators, DRs-direct and DRs-upper, with baselines across various conditions (apart from deterministic versus stochastic rewards from the main paper). We first investigate the performance under friendly logging (logging and evaluation policies are derived from the same deterministic policy, π1,det\pi_{1,det}), adversarial logging (logging and evaluation policies are derived from different policies π1,det\pi_{1,det}, π2,det\pi_{2,det}), and uniform logging (logging policy is uniform over all actions). Then we plot the performance in the small sample regime, where we aggregate the 108 conditions (6 logging policies, 9 datasets, deterministic/stochastic reward) at just 200 bandit samples.

In Table 7–Table 11, we compare the performance of DRs-direct and DRs-upper against baselines across various choices of reward predictors. We begin with using the best reward predictor type for each method (matching the setting of the main paper), and then consider each reward predictor in turn, across all estimators. We report the number of conditions where each estimator is statistically indistinguishable from the best, and the number of conditions where each estimator statistically dominates all others. DRs-upper is most often in the top group and most often the unique winner. DRs-direct is also better than snIPS, snDR, and switch. These results suggest that our shrinkage estimators are robust to different choices of reward predictors, and not just limited to the recommended set {η^≡0,z=w2}\{\hat{\eta}\equiv 0,z=w^{2}\}.

In Figure 6, we test the robustness of our proposed methods as we incorporate more reward predictors. Our practical suggestions is to use {η^≡0,z=w2}\{\hat{\eta}\equiv 0,z=w^{2}\} (shown as DRs-direct and DRs-upper in the figure). Here we also evaluate these methods when selecting from all reward predictors in the set {η^≡0,z≡1,z=w,z=w2,MRDR}\{\hat{\eta}\equiv 0,z\equiv 1,z=w,z=w^{2},\text{MRDR}\} (shown as DRs-direct (all) and DRs-upper (all) in the figure). For DRs-direct, the curves almost match, suggesting that it is quite robust. However, DRs-upper is less robust to including additional reward predictors.

At the end of appendix, we provide learning curves across all conditions. Dataset glass is excluded since we only ran it for a single sample size n=214n=214.

D.2 Experimental Details for Combinatorial Actions

We select the hyperparameter λ\lambda from the grid of 15 geometrically spaced values, with the smallest value 0.01×(w0.05)20.01\times(w_{0.05})^{2} and the largest value 100×(w0.95)2100\times(w_{0.95})^{2}, where w0.05w_{0.05} and w0.95w_{0.95} are the 0.05 and 0.95 quantiles of the weights w(xi,ai)w(x_{i},\mathbf{a}_{i}) on the logged data. We also add two boundary values λ=10−50\lambda=10^{-50} and λ=1030\lambda=10^{30} to include DM and DR-PI as special cases.