This paper introduces the confounded pure exploration transductive linear bandit () problem.
As a motivating example,  often online services cannot directly assign users to specific control or treatment experiences either for business or practical reasons. In these settings, naively comparing treatment and control groups that may result from self-selection can lead to biased estimates of underlying treatment effects.
Instead, online services can employ a properly randomized encouragement that incentivizes users toward a specific treatment.
Our methodology provides online services with an adaptive experimental design approach for learning the best-performing treatment for such encouragement designs.
We consider a more general underlying model captured by a linear structural equation and formulate pure exploration linear bandits in this setting. Though pure exploration has been extensively studied in standard adaptive experimental design settings, we believe this is the first work considering a setting where noise is confounded.  Elimination-style algorithms using experimental design methods in combination with a novel finite-time confidence interval on an instrumental variable style estimator are presented with sample complexity upper bounds nearly matching a minimax lower bound. Finally, experiments are conducted that demonstrate the efficacy of our approach.

§ INTRODUCTION

In this study, we present a methodology for adaptive experimentation in scenarios characterized by potential confounding.  Online services routinely conduct thousands of A/B tests annually [Kohavi et al., 2020]. In most online A/B/N experimentation, meticulous user-level randomization is essential to ensure unbiased estimates of treatment effects at the population level, commonly known as average treatment effects (ATE). In this setting, firms are often interested in understanding the treatment with the highest average outcome if presented to all members of the population.
However, in many settings firms may not be able to randomize, for example if a feature must be rolled out to all users for various business reasons [Mummalaneni et al., 2022].

In such instances,  users may choose to engage with a feature or not based on potentially unobservable preferences. Thus the resulting measured outcome may be correlated with the decision to engage in a specific feature. I.e. the underlying characteristics of the user confound the relationship between the decision to use the feature being evaluated and the effect of the feature. Thus, naively comparing the average outcome for users who engage with a feature with those who do not suffers from a (selection) bias. This setting is captured in Figure <ref>.

[state, align=center, rectangle] (x) at (0,0) $X$
User Choice;

[state, align=center, rectangle] (y) [right =of x] $Y$

[state, align=center, rectangle] (z) [left =of x] $Z$

[state, align=center] (u) [below =of x] $U$

(x) edge (y);
(z) edge (x);
[dashed] (u) edge (x);
[dashed] (u) edge (y);

Causal graph of the model.

A potential solution is for services to employ encouragement designs where users are presented with incentives that encourage users to engage with a specific feature [Bakshy et al., 2014, Bradlow, 1998, Elbers, 2023]. As a concrete example, many online

services have introduced membership levels with different offerings and prices available to all users.
Given a set of membership level options, the service is interested in knowing the counterfactual of which level has the optimal outcome (e.g., total revenue) if every user chooses to join that membership level.
In this setting, encouragements could be coupons or trials for corresponding membership levels.

In these settings, the firm can use intent-to-treat (ITT) estimates for the treatment effect which naively compare the average outcomes between the groups given different encouragements. In practice, given an encouragement a user may not engage with the corresponding feature choosing either the control or a different feature. Hence, the resulting ITT estimate may be a diluted estimate of the ATE [Angrist et al., 1996].
However, all is not lost: if the encouragement presented to a user is properly randomized, and the service guarantees that the encouragement only affects the outcome through the choice of user treatment, then the encouragement acts as an instrumental variable. Standard analysis from the econometrics and compliance literature show that two-stage least squares () estimators can then be used to provide consistent estimates of treatment effects.

At the same time, firms are also increasingly utilizing adaptive experimentation techniques, often known as pure exploration multi-armed bandits () algorithms [Lattimore and Szepesvári, 2020, Fiez et al., 2019], to accelerate traditional A/B/N testing. Pure exploration  techniques promise to deliver accurate inferences in a fraction of the time and cost as traditional methods. Similar to A/B testing, bandit methods assume users are properly randomized  and can fail to learn the optimal treatment if naively used and may be sample inefficient if they fail to take the confounded structure into account.

Contributions. In this work, we provide a methodology for experimenters seeking to use adaptive experimentation in settings with confounding where encouragements are available. We formulate this work in the more general and novel setting of confounded pure exploration transductive linear bandits () (Section <ref>). We present algorithms using experimental design for the  problem and analyze the resulting sample complexity. As we demonstrate, even in the simple multi-armed bandit setting described above, computing an effective sampling pattern requires using the machinery of linear bandits. Without knowledge of the underlying structural model, existing linear bandit approaches could lead to sub-optimal sampling. The main technical challenge we face is simultaneously improving our estimate of the structural model while designing with inaccurate estimates (Section <ref>). This approach crucially relies on our development of novel finite-time confidence bounds for two-stage least squares () style estimators that may be of independent interest (Section <ref>).
Moreover, we provide worst-case sample complexity lower bounds that are nearly matched
by our sample complexity upper bounds (<ref>).  Though the goal of this work is primarily theoretical, we empirically show the efficacy of our method over existing solutions (<ref>).

Related Works. Due to space limitations, we defer the related works in the linear/causal bandit, and econometrics literature to <ref>. We remark that there are only a handful of works in the multi-armed bandit literature that even consider confounding, and none that have studied the framework we pose.

§.§ General Problem Formulation

A confounded pure exploration transductive linear bandits () instance consists of finite collections of measurement vectors $\mathcal{Z}\subset \mathbb{R}^d$ and evaluation vectors $\mathcal{W}\subset \mathbb{R}^d$. At each time $t\in \mathbb{N}$, the learner selects $z_{t}\in \mathcal{Z}$ and observes a pair of noisy responses $x_t\in  \mathcal{X}\subseteq\mathbb{R}^d$ and $y_t\in\mathbb{R}$ generated via the structural equation model

\begin{equation}
\begin{aligned}
x_t=\Gamma^\top z_{t}+\eta_t,\quad y_t =x_t^\top\theta+\eps_t,
\end{aligned}
\label{eq:sem}
\end{equation}

where $\Gamma\in \mathbb{R}^{d\times d}$ and $\theta\in \mathbb{R}^d$ are model parameters.[We assume throughout that $\Gamma \in \mathbb{R}^{d\times d}$ is an invertible matrix.]
We define the history $\mc{H}_{t-1} = \{(z_s,x_s,y_s)\}_{s<t}$ and $\mathbb{E}_{t-1}[\cdot] = \mathbb{E}[\cdot|\mc{H}_{t-1}]$ denoting the conditional expectation under the filtration generated by $\mc{H}_{t-1}$. The noise  $\{\eta_t\}_{t=1}^{\infty}$ and $\{\eps_t\}_{t=1}^{\infty}$ satisfy the following set of assumptions unless otherwise noted.

We assume $\eps_t \mid \cH_{t-1}$ is 1-sub-Gaussian (and thus $\EE[\eps_t \mid \cH_{t-1}] = 0$).
Furthermore,  $\eta_t \mid \cH_{t-1}$ is $\sig_\eta^2$-sub-Gaussian vectors (and thus $\EE[\eta_t \mid \cH_{t-1}] = 0$), i.e.,

\begin{align*}
\forall \beta\in \RR, \max_{a:\|a\|_2\le 1} \EE[\exp(\beta\langle a, \eta_t \rangle)] \le \exp\del[2]{\fr{\beta^2\sig_\eta^2}{2}} ~.
\end{align*}

Goal. The objective is to identify $w^{\ast}:=\arg\max_{w\in\mc{W}}w^\top\theta$ with probability at least $1-\delta$ for $\delta \in (0, 1)$ while taking a minimum number of measurements.

In the setting where $\Gamma = I, \eta = 0$ and $\mathbb{E}_{t-1}[\epsilon_t|x_t] =0$,  our setting reduces to the standard pure exploration transductive linear bandit problem [Soare et al., 2014, Fiez et al., 2019].
In general, the joint noise process $\begin{bmatrix}\eta_t,\eps_t\end{bmatrix}$ may be dependent across the entries. In particular, we are allowing for the data generating process to be endogenous, meaning that $\mathbb{E}_{t-1}[\epsilon_t|x_t]\neq 0$. That is, $\epsilon_t$ can affect not just $y_t$, but also $x_t$ given a choice of $z_t$. The presence of endogeneity  is a key challenge in the  problem.

[Exclusion Restriction]
We assume that $\mathbb{E}_{t-1}[z_t\epsilon_t]=0$, or alternatively that $z_t$ is uncorrelated with $\epsilon_t$.

The variable $z_t$ is commonly referred to as an instrumental variable [Angrist et al., 1996]. We consider algorithms for the  problem that stop at a $\mc{H}_t$-measurable time $\tau \in \mathbb{N}$, and produce a recommendation $\widehat{w}\in \mathcal{W}$. The goal is $\delta$-PAC algorithms with efficient sample complexity guarantees.

We say an algorithm is $\delta$-PAC for a  problem with $\mathcal{W} ,\mathcal{Z}\subset \mathbb{R}^d$ if for all $\theta \in \mathbb{R}^d$ and $\Gamma\in \mathbb{R}^{d\times d}$, it holds that $\mathbb{P}_{\theta, \Gamma}(\widehat{w}\neq w^{\ast})\leq \delta$ for $\delta \in (0, 1)$.

§.§ Encouragement Designs

The  feedback model generalizes the classical compliance setting.

Compliance as a Special Case. In compliance problems, a decision-maker has access to a set of treatment that can be offered to users, while the users themselves have the option to accept the treatment they are presented or instead opt-in to a different treatment. The goal is to identify the treatment with the optimal average outcome if all users were to accept it.
Specifically, given a finite set $\mathcal{A}=\{1, 2, \dots, d\}$, a decision-maker presents user $t\in \mathbb{N}$ with an encouragement for a treatment $i\in \mathcal{A}$, the user then selects treatment $j\in \mathcal{A}$ where potentially $j\neq i$, and an outcome $y_t$ results. To capture compliance with the  framework, we set $\mc{Z}=\mc{X}=\mc{W}=\cbr{e_1,\cdots,e_d}$ and the parameter $\Gamma$ captures the probability of accepting a treatment given an encouragement for a potentially different treatment. Specifically, $\Gamma(i,j)=\pp{x_t=e_i \mid z_t=e_j}$, and a straightforward computation shows that
$x_t = \Gamma^{\top}z_t + \eta_t$ where $\mathbb{E}[\eta_t|z_t] = 0$ with

\begin{equation}
\eta_t=x_t-\begin{bmatrix}\pp{e_1\mid z_t},\cdots, \pp{e_d\mid z_t}\end{bmatrix}^\top.
\label{eq:compliance_noise}
\end{equation}

Moreover, $y_t = x_t^{\top}\theta+\epsilon_t$ gives the resulting reward, which is clearly correlated with the user choice so that $\cov(\eta_t, \epsilon_t)\neq 0$. Finally, $e_i^{\top}\theta = \theta_i$ gives the expected value of treatment $i$ over the population and our goal is to identify $w^{\ast}=\argmax_{e_i\in \mathcal{W}}e_i^{\top}\theta$.

[Our setting differs slightly from the traditional compliance setting based on a potential outcomes framework [Angrist et al., 1996]. Our setting is equivalent to one where we assume that there is a constant treatment effect. See Chapter 4 of [Angrist et al., 1996] for further discussion of the differences. Thus in our setting, learning $\theta_i$ and the local average treatment effect (LATE) on compilers are the same. ]
Note that when $\Gamma=I$, we automatically have that $\eta_t = 0$ and there is no confounding. This reduces to the standard MAB setting.

[$w^{\top}\theta$ and $w^{\top}\Gamma\theta$ for $w\in \{e_1,\dots, e_d\}$ ]
[$\mathbb{P}(\widehat{w}_t=w^{\ast})$ for the instance]
(a) A bar chart showing $\mathbb{E}[y|x=w]=w^{\top}\theta$ and $\mathbb{E}[y|z=w]=z^{\top}\Gamma\theta$ for all $w\in \mathcal{W}$. This chart shows that the optimal evaluation vector is $w^{\ast}=e_1=\argmax_{w\in \mathcal{W}}\mathbb{E}[y|x=w]$, while $e_6=\argmax_{w\in \mathcal{W}}\mathbb{E}[y|z=w]$ and consequently estimation based on this quantity is problematic. (b)  The probability of identifying $w^{\ast}=e_1$ for a collection of algorithms on the  instance from Section <ref>.

Motivating Compliance Example.
As a motivating compliance example representing the membership level discussion from the introduction, consider a location model that assumes each user $t\in \mathbb{N}$ arriving online has an underlying unobserved one-dimensional preference $u_t \sim \mathcal{N}(0,\sigma^2_u)$. If an algorithm presents the user with encouragement $z_{I_t} = e_{I_t}$ for $I_t\in \mathcal{A}$, then the user selects into the membership level given by $J_t = \min_{j\in \mathcal{A}}|I_t + u_t - j|$ so that $x_t=e_{J_t}$.
This process captures a user being more likely to opt-in to membership levels that are closer to the encouragement that they were presented. The outcome is then given by $y_t = x_t^{\top}\theta + u_t$.

We conduct an experiment with this problem instance (see Fig. <ref> and <ref> for more details).
Specifically, $d=6$, $\theta=\begin{bmatrix}1& -0.95&0&0.45 &0.95& 0.99\end{bmatrix}$, and $\sigma_u^2=0.35$. Observe that $w^{\ast}=e_1=\argmax_{w\in \mathcal{W}}w^{\top}\theta$.
An upper confidence bound (UCB) selection strategy is simulated that maintains estimates of the mean reward of each encouragement $i\in \mathcal{A}$, namely $\widehat{\mu}_{i,t} = \sum_{s=1}^t \mathbf{1}\{z_t = e_i\} y_t$, and then pulls the one with the highest UCB. The UCB selection strategy is combined with a pair of recommendation strategies.
The UCB-OLS algorithm estimates the mean reward of each treatment using an OLS estimator, namely $\widehat{\theta}_{\texttt{LS}}^{i, t}  = \sum_{s=1}^t \mathbf{1}\{x_t = e_i\} y_t / \sum_{s=1}^t \mathbf{1}\{x_t = e_i\}$, and recommends $\argmax_{a\in \mathcal{A}}\widehat{\theta}_{\texttt{LS}}^{i, t}$. Moreover, the UCB-IV algorithm uses an instrumental variable-estimator (see the next section) that incorporates knowledge of $\Gamma$ similar to 2SLS to deconfound estimates of the mean rewards of each treatment and recommends the treatment with the maximum estimate. The results over 100 simulations are shown in Fig. <ref>.
UCB-OLS completely fails to identify $\theta_1 = \arg\max_{i\in d} \theta_i$ due to a biased estimate, whereas UCB-IV does better. However, UCB-IV methods seem to have a constant probability of error. To see why, note that the expected reward from pulling $z = e_{i}$ is $e_i^{\top}\Gamma\theta$. These values are plotted in orange in Figure <ref>. In particular, with some constant probability, UCB zeroes in on arm 6 becauses of the mean estimates on the $z$'s, and as a result fails to give enough samples to learn that arm 1 is indeed the best. In contrast, our proposed method CPEG, Algorithm <ref> manages to find the best arm with significantly higher probability.

Notation. Let $\Delta(\mathcal{Z})=\{\lambda \in \mathbb{R}^{|\mathcal{Z}|}: \lambda \geq 0, \sum_{z\in \mathcal{Z}}\lambda_z =1\}$ denote the set of probability distributions over the set $\mathcal{Z}$. Given a distribution $\lambda \in \Delta(\mathcal{Z})$ and matrix $\Gamma \in \mathbb{R}^{d\times d}$, define the operator $A(\lambda, \Gamma):=\sum_{z\in \mathcal{Z}}\lambda_z \Gamma^{\top} zz^{\top}\Gamma$. Given  $Z \in \mathbb{R}^{T\times d}$ and  $\Gamma \in \mathbb{R}^{d\times d}$,  define the operator $\bar{A}(Z, \Gamma):=\sum_{t=1}^T \Gamma^{\top} z_tz_t^{\top}\Gamma = \Gamma^{\top}Z^{\top}Z\Gamma$ where $z_t\in \mathbb{R}^d$ denotes row $t$ of $Z$. Given a vector $x\in \mathbb{R}^d$ and a symmetric positive-definite matrix $A\in \mathbb{R}^{d\times d}$ we let $\|x\|^2_{A} = x^{\top} A x$. We adopt the standard notation that $(a\vee b )\equiv \max\{a, b\}$ and $(a\wedge b)\equiv \min\{a, b\}$ for $a, b\in \mathbb{R}$. $\sigma_{\min}(A), \sigma_{\max}(A)$ denote the minimum and maximum singular value of a matrix $A$.
We denote by $\polylog(x_1,\ldots,x_n)$ any polylogarithmic factors of $x_1,\ldots,x_n$.

§ ESTIMATORS AND INFERENCE

We now present estimators for the unknown parameter $\theta$ and prove the associated statistical properties. The estimators discussed in this section are critical to our algorithmic solution outlined in Section <ref>.

§.§ Estimators

Before describing our solution concept, we quickly review potential options for estimating $\theta$ based on a dataset $Z_T = [z_1, \cdots, z_T]^{\top}\in\mathbb{R}^{T\times d}, X_T = [x_1, \cdots, x_T]^{\top}\in\mathbb{R}^{T\times d}$, $Y_T = [y_1, \cdots, y_T]^{\top}\in \mathbb{R}^T$, assumed to be generated according to the model in Eq. (<ref>).
Recall that the ordinary-least-squares () estimator for $\theta$ is given by

\begin{align}
\textstyle \widehat{\theta}_{\texttt{LS}}:=\argmin_{\widehat{\theta}\in \mathbb{R}^d}\sum_{t=1}^{T}(y_t-x_t^{\top}\widehat{\theta})^2= (X_T^\top X_T)^{-1}X_T^\top Y_T=\theta + (X_T^\top X_T)^{-1}X_T^\top\varepsilon_T.
\label{eq:ls}
\end{align}

Observe that $\widehat{\theta}_{\texttt{LS}}$ is potentially a biased and inconsistent estimator for $\theta$ in the presence of endogenous noise since $\mathbb{E}[\epsilon_t|x_t]\neq 0$.
To remediate this problem, we define a general class of estimators that includes several standard estimators.
Given an invertible matrix $\Psi\in \mathbb{R}^{d\times d}$, let $\bar{X}_T := Z_T\Psi$, and consider corresponding estimators termed $\Psi$- estimators of the form

\begin{equation}
\widehat{\theta}_{\Psi}:=(\bar{X}_T^\top \bar{X}_T)^{-1}\bar{X}_T^\top Y_T
=( \Psi^{\top} Z_T^{\top} Z_T\Psi)^{-1} \Psi^{\top}Z_T^\top Y_T = ( Z_T^{\top} Z_T\Psi)^{-1} Z_T^\top Y_T.
\end{equation}

When $\Psi = I$ we recover the  estimator. In the rest of the paper, we will focus on two different potential options for $\Psi$.

Case 1: . $\Psi = \Gamma$.
To begin, observe that the structural equation model from Eq. (<ref>) can be combined by substituting the second equation into the first to obtain the reduced form

\begin{equation}
y_t= z_{t}^{\top}\Gamma\theta + \eta_t^{\top}\theta+\epsilon_t.
\label{eq:combined_struct}
\end{equation}

Since $z_t$ is independent of the i.i.d. process $\eta_t^{\top}\theta+\epsilon_t$, the least squares estimator which regresses $y_t$ onto $z_t^{\top}\Gamma$ is unbiased for estimation of $\theta$ and given by

\begin{equation}
\widehat{\theta}_{\texttt{oracle}}
=( \bar{X}_T^{\top} \bar{X}_T)^{-1} \bar{X}_T^\top Y_T=( Z_T^{\top} Z_T\Gamma)^{-1} Z_T^\top Y_T.
\label{eq:oracle_estimator}
\end{equation}

This estimator will be used to design our general solution concept presented in Section <ref>. Of course in practice we cannot expect to know $\Gamma$, but we may be able to estimate it.

Case 2: . $\Psi = \widehat{\Gamma}$. We consider a setting where $\widehat{\Gamma}$ is an (unbiased) estimator of $\Gamma$, learned using least-squares from an independent dataset $Z_{T_1} = [z'_{1}, \cdots, z'_{T_1}], X_{T_1}= [z'_{1}, \cdots, z'_{T_1}]$ collected non-adaptively.[Formally, we say that a set of data $(Z_T, X_T, Y_T)$ generated via the model in Eq. (<ref>) is collected non-adaptively from an experimental design if $z_t$ is $\mathcal{H}_0$ measurable for all $1\leq t\leq T$.] That is, $\widehat{\Gamma} = (Z_{T_1}^{\top}Z_{T_1})^{-1}Z_{T_1}^{\top}X_{T_1}$ and:

\begin{align*}
\widehat{\theta}_{\psls} = (\widehat{\Gamma}^{\top}Z_{T}^{\top}Z_{T}\widehat{\Gamma})^{-1}\widehat{\Gamma}^{\top}Z_{T}^{\top}Y_{T}.
\end{align*}

We refer to the resulting estimator as a pseudo two stage least squares () estimator. The main advantage of the  estimator over standard  (given in <ref>) is easier inference since now $\{\epsilon_{t}\}_{t\leq T}$ of our dataset is independent of the measurements of the first dataset $Z_{T_1}, X_{T_1}$. In the econometrics literature, such an estimator is referred to as a two-sample  estimator [Inoue and Solon, 2010].

§.§ Confidence Intervals

In the section that follows, we develop a general algorithmic approach that relies on experimental design aimed at reducing the uncertainty in our estimates of the optimal treatment. To this end, we first develop finite-time confidence intervals for estimators presented in the previous section given data generated according to the model in Eq. (<ref>) and collected from non-adaptive designs.

We begin by characterizing the properties of the noise structure in the combined model of Eq. (<ref>) with the following set of results.

Under Assumption <ref>, the noise process $\nu:=\eta^{\top}\theta+\epsilon$ is $\sigma_\nu^2$-sub-Gaussian where $\sigma_\nu^2=2(\sig_\eta^2\|\th\|_2^2 + 1)$, specifically when the instance is compliance, $\sigma_\nu^2=2(4\|\theta\|_2^2+1)$.

Oracle Confidence Interval.
As in the last section, we assume that we have access to a dataset $(Z_T, X_T, Y_T)$ generated according to Eq. <ref> and collected non-adaptively.
Given Lemma <ref>, it can be shown that $w^{\top}\widehat{\theta}_{\texttt{oracle}}$ is a sub-Gaussian random variable satisfying the following.

With probability at least $1-\delta$ for $\delta\in (0, 1)$ and $w \in \mathbb{R}^d$,

\begin{align*}
&|w^{\top}(\widehat{\theta}_{\texttt{\emph{oracle}}} - \theta)| \leq \sqrt{2\sigma_\nu^2\|w\|^2_{\bar{A}(Z_T, \Gamma)^{-1}} \log\left(2/\delta\right)},
\end{align*}

where $\sigma_\nu^2$ is the sub-Gaussian parameter of the noise $\nu:=\eta^{\top}\theta+\epsilon$ characterized in Lemma <ref>.

The proof of this result is in Appendix <ref>.
P-2SLS Confidence Interval.
We now present a novel finite-time confidence interval for the  estimator. As discussed in the previous section with respect to this estimator, we assume access a set of data $(Z_{T_1}, X_{T_1})$ generated according to Eq. (<ref>) and collected non-adaptively for the purpose of estimating $\Gamma$. Moreover, assume access to a separate set of data $(Z_{T_2}, X_{T_2}, Y_{T_2})$ generated according to Eq. (<ref>) and collected non-adaptively for the purpose of estimating $\theta$.

Suppose that $\widehat{\Gamma} = (Z_{T_1}^{\top}Z_{T_1})^{-1}Z_{T_1}^{\top}X_{T_1}$ and $\widehat{\theta}_{\psls} = (\widehat{\Gamma}^{\top}Z_{T_2}Z_{T_2}\widehat{\Gamma})^{-1}\widehat{\Gamma}Z_{T_2}^{\top}Y_{T_2}$. Then, for any $w \in \mathcal{W}$, with probability at least $1-\delta$ for $\delta \in (0,1)$,

\begin{align*}
|w^{\top}(\widehat{\theta}_{\emph{\texttt{P-2SLS}}} - \theta)| \leq\vcn{w}_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4}{\delta}}}+\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}_2\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta/4}},
\end{align*}

\begin{align*}
\overline{\log}\rbr{Z_{T},\delta}:=8d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(Z_{T_1}^\top Z_{T_1}))} } + 16\ln\del{\fr{2\cd6^d}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(Z_{T_1}^\top Z_{T_1})}}}.
\end{align*}

The proof is presented in Appendix <ref>.
Observe that the first term in the  estimator confidence interval given by $\sqrt{2\sigma_\nu^2\|w\|^2_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}} \log(4/\delta)}$ matches the  estimator confidence interval in Lemma <ref> when $\widehat{\Gamma} = \Gamma$. The second term scaling like $\mathcal{O}(\|w\|_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\|\theta\|_2\sigma_{\eta}\sqrt{d + \log(1/\delta)})$, is an upper bound on the approximation error $w^{\top}(\widehat{\Gamma}^{-1}\Gamma - I)\theta$ for any $w\in \mathbb{R}^d$, assuming that $\widehat{\Gamma}$ is learned from an  estimator (see <ref> for details).

We will see that the form of this confidence interval is particularly convenient for our algorithmic approach given in Section <ref>. In particular, the form of the variance $\|w\|^2_{\bar{A}(Z_{T_2}, \Gamma)^{-1}}$ on the first term only depends on a design over instruments. Thus, we can choose an experimental design over $Z$'s which reduces this variance optimally.

In practice we expect the first stage of samples, $(Z_{T_1}, X_{T_1})$ to be collected from either a burn-in period or from existing historical data. We remark that assuming two stages of samples is common in the orthogonal and double machine learning for estimating nuisance parameters in the data generating process (e.g. $\Gamma$) [Chernozhukov et al., 2018]. Our result matches the existing literature on the asymptotic variance of two sample  estimators (e.g., Theorem 1 of [Inoue and Solon, 2010]).

The asymptotic variance of standard  is known to involve a factor $\sigma_{\epsilon}^2$, instead of $\sigma_{\nu}^2$ as we have [Greene, 2003]. Recent work by [Della Vecchia and Basu, 2023] shows a variance involving $d\sigma_{\epsilon}^2$. However, it's unclear how to use the form of their confidence interval directly for experimental design. In addition, their work is not sufficiently general to handle the general forms of noise that we consider in <ref>.

§ ADAPTIVE EXPERIMENTAL DESIGN ALGORITHMS

We now present adaptive experimental design algorithms for the  problem.
Our main insight utilizes Eq. <ref> by plugging the model for $x$ into the top equation resulting in the relationship

\[y = z^{\top}\Gamma\theta + \theta^{\top}\eta + \epsilon.\]

When $\Gamma$ is known, by Eq. <ref>, we see that  reduces to a standard pure exploration transductive linear bandit problem where the measurement set is given by $\{\Gamma^{\top}z\}_{z\in \mathcal{Z}}\subset \mathbb{R}^d$, the evaluation set is $\mathcal{W}\subset \mathbb{R}^d$, and the feedback model is given by $y=v^{\top}\theta+\nu$ where the noise $\nu = \theta^{\top}\eta + \epsilon$ is sub-Gaussian and as before the goal is to identify $\arg\max_{w\in \mc{W}} w^{\top}\theta$.
An existing approach to this problem is given by the  algorithm [Fiez et al., 2019], which we use as the basis of our approach. Addressing the case of unknown $\Gamma$ is our major algorithmic contribution, where we develop solutions to improve our estimate of $\Gamma$ and learn $w^{\ast}$ simultaneously. As a warm-up to this approach, we first consider the setting when $\Gamma$ is known.

§.§ Warm-Up: Known Structural Model

Algorithm <ref> assumes a parameter $L_\nu$, which acts as an upper bound on the sub-Gaussian constant of the noise $\nu = \theta^{\top}\eta + \epsilon$. In each round $k$, an active set of potentially optimal vectors $\widehat{\mathcal{W}}_k\subset \mathcal{W}$ is maintained. CPEG aims to sample in such a way that reduces the uncertainty of the estimates on the gaps $(w-w')^{\top}\theta$ for each pair $w, w'\in \widehat{\mathcal{W}}_k$ maximally each round. In any given round the algorithm takes $N_k$ samples $Z_{N_k}$, the confidence interval of <ref> shows that the error in estimating $(w-w')^{\top}\theta$ scales with $\|w-w'\|^2_{(\Gamma^{\top}Z_{N_k}^{\top}Z_{N_k}\Gamma)^{-1}}$. This motivates utilizing an experimental design approach where we choose a distribution $\lambda_k\in \Delta(\mathcal{Z})$ to minimize $\max_{w, w'\in \widehat{\mathcal{W}}_k}\|w-w'\|_{(\sum_{z\in \mc{Z}} \lambda_z \Gamma^{\top}zz^{\top}\Gamma)^{-1}}$. The number of resulting samples taken from this design $N_k$ is chosen to guarantee that the confidence interval of Lemma <ref> is less than $2^{-k}$. Then, the elimination step in Line 8 guarantees that all $w\in \mathcal{W}$ such that $(w^{\ast}-w)^{\top}\theta> 2\cdot 2^{-k}$ are then eliminated from the active set by round $k+1$ of the procedure. To actually choose our samples, as is common in this literature [Fiez et al., 2019], we use an efficient rounding procedure, ROUND that requires a minimum number of samples $r(\omega)$.

Sample Complexity Guarantee.
The sample complexity of Algorithm <ref> depends on the following problem-dependent quantity $\rho^{\ast}(\gamma)$ that captures the underlying hardness of a problem instance in terms of $(\mathcal{W}, \mathcal{Z}, \Gamma, \theta)$, when $\gamma=0$, we abbreviate $\rho^{\ast}(0)=\rho^{\ast}$,

\begin{equation}
\rho^{\ast}(\gamma)=\min_{\lambda\in\Delta(\mc{Z})}\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}}\frac{\|w^{\ast}-w\|^2_{(\sum_{z\in \mathcal{Z}}\lambda_z\Gamma^{\top} zz^{\top}\Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2\vee\gamma^2}.
\label{eq:rhostar}
\end{equation}

[H]CPEG:Confounded pure exploration with $\Gamma$

[1]
Input $\mc{Z}, \mc{W}, \Psi=\Gamma, \delta, L_\nu\ge\sigma_\nu^2, \omega$,
Initialize: $k=1, \mc{W}_1=\mc{W}, \zeta_1=1$
$f(w,w',\Gamma,\lambda):=\|w-w'\|^2_{(\sum_{z\in \mathcal{Z}}\lambda_z\Gamma^{\top} zz^{\top}\Gamma)^{-1}}$
$N_k:=\lceil{2(1+\omega)2^{2k}\rho(\mc{W}_k)L_\nu\log\rbr{4k^2\abr{\mc{W}}/\delta}}\rceil\vee r(\omega)$

Pull arms in $Z_{N_k}=\text{ROUND}(\lambda_k,N_k)$ and observe $Y_{N_k}$. Compute $\widehat{\theta}_{\Gamma}^k=\rbr{Z_{N_k}^\top Z_{N_k}\Gamma}^{-1}Z_{N_k}^\top Y_{N_k}$

$\mc{W}_{k+1}=\mc{W}_k\backslash\{w\in\mc{W}_k|\exists w'\in\mc{W}_k,\langle w'-w, {\widehat{\theta}_{\Gamma}^k\rangle }>2^{-k}\}$, $k\gets k+1$
Output: $w\in\mc{W}_k$

Algorithm <ref> is $\delta$-PAC and terminates in at most $c(1+\omega)L_\nu\rho^{\ast}\log(1/\delta)+cr(\omega)$ samples, where $c$ hides logarithmic factors of $\Delta:=\min_w\langle w^{\ast}-w, \theta\rangle$ and $\abr{\mc{W}}$, as well as constants.

The proof of this result is in Appendix <ref>.
In the unconfounded case when $\Gamma=I$ and $\eta=0, \mathbb{E}_{t-1}[\epsilon_t|x_t] = 0$ this matches the sample complexity of  [Fiez et al., 2019].  In particular, for the case where $\mc{Z} = \mc{X} = \mc{W}$, the problem further reduces to a standard multi-armed bandit, and if $\epsilon$ is 1-sub-Gaussian noise, [Soare et al., 2014] shows that $\rho^{\ast} = O(\sum_{i=2}^d (\theta_{1} - \theta_i)^{-2}))$, which is the optimal sample complexity of best-arm identification for multi-armed bandits. The following lemma shows that the conditioning of $\Gamma$ can have a strong impact on the resulting sample complexity.

For the compliance setting, we have
$\min_{\lambda\in \Delta^d}\max_{j,j'} \|e_{j} - e_{j'}\|^2_{(\sum_{i=1}^d \lambda_i \Gamma^{\top} e_i e_i^{\top}\Gamma)^{-1}} \leq d\max_{j,j'}\|\Gamma^{-1}(e_{j} - e_{j'})\|_2^2$. Furthermore,
$\rho^{\ast} \leq \frac{d \sigma_{\min}^2(\Gamma)^{-1}}{\Delta_{\min}^2}$.

To further illustrate the impact of $\Gamma$, imagine an extreme setting where $\Gamma = (1-\epsilon)/d \mathbf{11}^{\top} + \epsilon I$ and $\epsilon\approx 0$, i.e. $\Gamma$ is a perturbation of $1/d \mathbf{11}^{\top}$. It's straightforward to show that the upper bound in the first display of Lemma <ref> is of the order $O(d\epsilon^{-2})$ (this is also a lower bound - see Appendix <ref>). In particular, the upper bound on the sample complexity is of the form $d\epsilon^{-2}/\Delta_{\min}^2$. This is in sharp contrast to the linear bandit case, when $\Gamma=I$ and we are guaranteed a sample complexity of no more than $d/\Delta_{\min}^2$ samples.
To gain some intuition, regardless of the choice of $\lambda$, $\sum_{i \leq d} \lambda_i \Gamma^{\top}e_ie_i^{\top}\Gamma \approx \Gamma$. As a result, $\rho^{\ast} \rightarrow \infty$ as $\epsilon\rightarrow 0$. Intuitively in the limit, regardless of which instrument $i \leq d$ is being pulled, the resulting distribution on the treatments is uniform (the instruments are weak). Thus, it is impossible to deconfound the measurement noise, and recover an estimate of $\theta$. This is a phenomenon which does not arise in the standard multi-armed bandit case with unconfounding.

We also consider a setting where instead of given $\Gamma$ directly, we are given an estimate $\widehat{\Gamma}$ of $\Gamma$ based on offline data. We discuss such an adaptation of Algorithm <ref> to this setting in Appendix <ref> and provide a sample complexity which reflects the error in $\Gamma$ (scaling with $\rho^{\ast}(\gamma)$ for $\gamma>0$). We remark that this result is subsumed by the approach of Section <ref> and so we omit it in the main text.

Lower bound. Due to the noise model from confounding and the dependence of the noise $\theta^{\top}\eta + \epsilon$, the instance-dependent lower bounds of [Fiez et al., 2019] do not immediately apply. We develop a lower bound tailored for the confounding setting that nearly match the upper bounds of our algorithms. What's more, our lower bound illustrates the additional difficulty that arises from confounding by an additional factor of $d^2$ compared to the standard transductive linear bandit problem in the most general setting where entries of $\eta$ are sub-Gaussian, but not necessarily independent nor bounded. Due to space limit, we defer it to <ref>.

§.§ Fully Unknown Structural Model

We now consider the setting where $\Gamma$ is fully unknown. The difficulty of this setting is that the data collection process needs to support both estimation of $\Gamma$ and $\theta$ simultaneously.

Our algorithm, built upon Algorithm <ref>, is summarized in Algorithm <ref>.
At its core, each phase of the algorithm is divided into two sub-phases, for estimating $\Gamma$ and $\theta$ respectively.
Specifically, the second sub-phase is essentially same as Algorithm <ref> with $\widehat{\Gamma}_{k}$ in place of $\Gamma$ where $\widehat{\Gamma}_k$ is estimated from the first sub-phase.
The main novelty of our algorithmic design lies in the first sub-phase, which resolves the challenge of performing the optimal design for estimating $\Gamma$.
To explain this challenge, the confidence interval for P-2SLS estimators of <ref> indicates that one should pull arms so that we control both $D_2 := \max_{w,w'}\|w-w'\|^2_{\bar A(Z_{T_2},\kGam)}$ (error from $\hth_{\mathsf{P-2SLS}}$) and $D_1 := \max_{w,w'}\|w-w'\|^2_{\bar A(Z_{T_1},\kGam)}$ (error from $\kGam$) to be below the target error $O(\zeta_k^2)$ at each phase (ignoring unimportant factors for discussion).
Controlling $D_2$ is trivial, which is done in the second sub-phase as we described above.

However, for $D_1$, a similar strategy cannot be done because the estimate $\kGam$ is computed directly by sampling arms in $Z_{T_1}$. That is, the ideal design, based on which we will collect data points $z_1,\ldots, z_n$, requires access to the random matrix $\kGam$ that can only be computed after sampling $z_1,\ldots,z_n$,.
This creates a cycle that seems impossible to resolve.
Such an issue, to our knowledge, has not been seen in existing work on pure exploration, and thus resolving it is our key technical contribution.

Our solution is to compute the design based on $\kGam$ from the previous phase.
We then perform a doubling trick where we double the sample size (while following the computed design) until $D_1$ becomes smaller than the target error $O(\zeta_k^2)$.
The intuition is that in later phases the estimate $\kGam$ from the previous phase will be accurate enough to ensure that the design is efficient.
Note that this novel algorithm induces extra randomness in how many samples we end up collecting in the first sub-phase, which remains random even after conditioning on the history, unlike the second sub-phase.
This makes the analysis challenging, which we describe after the main result.

Input $\mc{W}, \delta, \zeta, \widehat{\Gamma}, \omega,L_\nu$
$N_{2}=\ur{2(1+\omega)\zeta^{-2}\rho(\mc{W})L_\nu\log\rbr{\fr{4\abr{\mc{W}}}{\delta}}}\vee r(\omega)$
get $N_{2}$ samples per design $\widehat{\lambda}$ denoted as $\cbr{Z_{2},X_{2},Y_{2}}$ via ROUND
update $\widehat{\theta}_{\texttt{P-2SLS}}=(\widehat{\Gamma}^{\top}Z_{2}^{\top}Z_{2}\widehat{\Gamma})^{-1}\widehat{\Gamma}^{\top}Z_{2}^{\top}Y_{2}$
Output: $\widehat{\theta}_{\texttt{P-2SLS}}$

Our algorithm additionally employs the so-called E-optimal design to ensure that the covariance matrix of the collected data used to estimate $\Gamma$ is well-conditioned.
This conditioning is required to ensure that $\kGam$ concentrates fast enough to $\Gamma$ as shown in the analysis.
The E-optimal design is a well-known design objective in experimental design that aims to maximize the smallest singular value: $\lam^*_E := \arg \min_{\lam\in \Delta(\cZ)} \sig_\tmax(V^{-1}(\lam))$, where $V = \sum_{z\in\cZ} \lam_z z z^\T$.
We denote $\kappa_0^{-1} := \sigmax(V^{-1}(\lam^*_E)) = \sigmin^{-1}(V(\lam^*_E))$ as the smallest singular value achieved by the E-optimal design.

[t]CPEUG: Confounded pure exploration with with unknown $\Gamma$

Input $\mc{Z}, \mc{W}, \delta, L_\nu\ge\sigma_\nu^2, L_\eta\ge\sigma_{\eta}^2, \omega, \gammin\le\lammin(\Gam), \lambda_E,\kappa_0$
Initialize: $k=1, \mc{W}_1=\mc{W},\widehat{\Gamma}_0=\perp, \zeta_1=1$
Define $f(w,w',\Gamma,\lambda):=\|w-w'\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}$,
$M:=\fr{32L_\eta}{\gammin^2\sigma_{\min}\rbr{A(\lambda_E, I)}}\vee1$, $\dt_{\ell}:=\fr{\delta}{4\ell^2}$
$\widehat{\Gamma}_k=\Gamma-\mathtt{estimator}\rbr{\mc{W}_k, \widehat{\Gamma}_{k-1}, \zeta_k, \delta/k^2, \omega,\lambda_E, M, L_\eta}$ Step 1: update $\widehat{\Gamma}$
$\widehat{\theta}_{\texttt{P-2SLS}}^k=\theta-\mathtt{estimator}\rbr{\mc{W}_k, \delta/k^2, \zeta_k, \widehat{\Gamma}_{k}, \omega,L_\nu}$ Step 2: update $\widehat{\theta}$
$\mc{W}_{k+1}=\mc{W}_k\backslash\cbr{w\in\mc{W}_k\mid \exists w'\in\mc{W}_k,\text{s.t.},\ibr{w'-w}{\widehat{\theta}_{\texttt{P-2SLS}}^k}>\zeta_k}$ Step 3: elimination
$k\leftarrow k+1$, $\zeta_k=2^{-k}$
Output: $\mc{W}_k$

Input $\mc{W}, \widehat{\Gamma}, \zeta, \delta, \omega,\lambda_E, M, L_\eta$
Define $\mathtt{Stop}\rbr{\mc{W},Z,\Gamma,\dt}:=\max_{w,w'\in\mc{W}}\vcn{w-w'}_{\bar{A}(Z, \Gamma)^{-1}}\vcn{\theta}_2\sqrt{L_\eta\overline{\log}\rbr{Z,\delta}}$
Initialize $\ell=1$, $N_{0,0}=0$ doubling trick initialization
$\ell=1$ or $\mathtt{Stop}\rbr{\mc{W},Z_{0,\ell},\widehat{\Gamma}',\dt_{\ell}}>1$
get $2^{\ell-1}\rbr{r(\omega)\vee\fr{2}{\kappa_0}}$ samples denoted as $\cbr{Z_{0,\ell},X_{0,\ell},Y_{0,\ell}}$ per design $\lambda_E$ via ROUND
Update $\widehat{\Gamma}'$ by  on $\cbr{Z_{0,\ell},X_{0,\ell}}$, $\ell\leftarrow \ell+1$
$N'=\drd{4gdM\ln\del{1+2M\rbr{d+L_z^2}+2M2gdM}+8M\ln\del{\fr{2\cd6^d}{\dt}}\vee r(\omega)}$
$\ell=1$ or $\mathtt{Stop}\rbr{\mc{W},Z_{0,\ell}\cup Z_{1,\ell},\widehat{\Gamma},\dt_{\ell}}>\zeta$
$N_{1,\ell}=2^{\ell}N'$ doubling trick update
get $N_{1,\ell}$ samples per $\tilde{\lambda}$ denoted as $\cbr{Z_{1,\ell},X_{1,\ell},Y_{1,\ell}}$ via ROUND
$N_{0,\ell}=\urd{2gdM\ln\del{M\rbr{d+N_{1,\ell}+L_z^2}}+4M\ln\del{\fr{2\cd6^d}{\dt_{\ell}}}\vee r(\omega)\vee\fr{2}{\kappa_0}}$
get $(N_{0,\ell}-N_{0,\ell-1})$ samples per $\lambda_E$ augmented to $\cbr{Z_{0,\ell-1},X_{0,\ell-1}}$ and get $\cbr{Z_{0,\ell},X_{0,\ell}}$
Update $\widehat{\Gamma}'$ by  on $\cbr{Z_{0,\ell}\cup Z_{1,\ell}, X_{0,\ell}\cup X_{1,\ell}}$, $\ell\leftarrow \ell+1$
Output: $\widehat{\Gamma}'$

We present our analysis result <ref> where we show that, even without knowledge of $\Gamma$, the sample complexity scales with the key problem difficulty $\rho^*$ almost matching the sample complexity of Algorithm <ref> which relies on knowledge of $\Gamma$.

Algorithm <ref> is $\delta$-PAC and terminates in at most

\begin{align*}
\end{align*}

pulls, ignoring both of the additive and multiplicative logarithms of $\Delta, \abr{\mc{W}}, \rho^*, \rho_0, M$, where

\begin{align*}
\rho_0=\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}}\|w^{\ast}-w\|^2_{(\sum_{z\in \mathcal{Z}}\lambda_{E,z}\Gamma^{\top} zz^{\top}\Gamma)^{-1}},\text{ and }
M=\fr{32L_\eta}{\gammin^2\sigma_{\min}\rbr{A(\lambda_E, I)}}\vee1.
\end{align*}

Note that $\rho_0$ does not get hurt by $\langle w^{\ast}-w, \theta\rangle$, ($\rho^*$ does). It comes from the fact that in the first phase, we initialize that algorithm with E-optimal design.

The challenge of the analysis can be summarized in two-fold.
First, since the concentration result in <ref> is w.r.t. $\kGam$, we need to analyze how the random matrix $\kGam$ concentrates around $\Gam$ and how this impacts the sample complexity.
For this, we develop a novel concentration inequality that relates the confidence width involving $\hGam$ from <ref> with the same quantity involving $\Gam$ in place of $\hGam$.
Second, our algorithm creates a long-range error propagation, which is highly nontrivial to analyze.
To see this, the quality of $\hGam_k$ is affected by the design objective function $\max_{w,w'} f(w,w',\hGam_{k-1},\lam)$, which depends on the error of the estimate $\hGam_{k-1}$ from the previous phase.
This error is, in turn, affected by the error of $\hGam_{k-2}$ by the same mechanism.
This is repeated all the way back to the first phase.
Thus, any abnormal behavior from the first iteration will have a cumulative impact to even the end.
In our analysis, we successfully analyze how the error is propagated from the previous iterations, which forms a complicated recursion.
Resolving this recursion is our key novelty in the analysis.

Our algorithm requires knowledge of a lower bound $\gammin$ of $\lammin(\Gam)$. The knowledge of $\gammin$ is for simplicity only as one can obtain such a lower bound that is at least half of the true value $\lammin(\Gam)$ via an efficient sampling procedure that we describe in <ref>.

Experiments. We provide experiments for the instance of Section <ref> in the Appendix <ref>. The experiments show that our approach is more sample efficient than natural passive baselines (e.g. A/B testing), or naively applying existing Pure-Exploration linear bandit methods and performs similarly to the oracle complexity.

§ CONCLUSION

This work introduces the  problem in which the learning protocol is characterized by a linear structural equation model governed by parameters $\Gamma$ and $\theta$. We provide a general solution that simultaneously estimates the structural model while optimally designing to learn the best-arm.
The key ideas behind our approach are based on linear experimental design techniques, an instrumental variable estimator whose variance can be controlled by the design, and novel finite-time confidence intervals on this estimator. This paper presents a number of directions for future work including considering situations where the $d_z\neq d_x$, analysis to improve the dependence on the underlying noise variance, and the pursuit of a tight information-theoretic instance-dependent lower-bound.
We hope that this line of work motivates increased discussion of the real impact of confounding on applicability of adaptive experimentation.

§ ACKNOWLEDGEMENTS

Kwang-Sung Jun and Yao Zhao were supported in part by the National Science Foundation under grant CCF-2327013.

[Abbasi-Yadkori et al., 2011]
Abbasi-Yadkori, Y., Pál, D., and Szepesvári, C.
Improved algorithms for linear stochastic bandits.
Advances in neural information processing systems, 24, 2011.

[Allen-Zhu et al., 2021]
Allen-Zhu, Z., Li, Y., Singh, A., and Wang, Y.
Near-optimal discrete optimization for experimental design: A regret
minimization approach.
Mathematical Programming, 186:0 439–478, 2021.

[Angrist et al., 1996]
Angrist, J. D., Imbens, G. W., and Rubin, D. B.
Identification of causal effects using instrumental variables.
Journal of the American Statistical Association, 910
(434):0 444–455, 1996.

[Bakshy et al., 2014]
Bakshy, E., Eckles, D., and Bernstein, M. S.
Designing and deploying online field experiments.
In International Conference on World Wide Web, pages 283–292,

[Bareinboim et al., 2015]
Bareinboim, E., Forney, A., and Pearl, J.
Bandits with unobserved confounders: A causal approach.
Advances in Neural Information Processing Systems, 28, 2015.

[Bilodeau et al., 2022]
Bilodeau, B., Wang, L., and Roy, D.
Adaptively exploiting d-separators with causal bandits.
Advances in Neural Information Processing Systems,
35:0 20381–20392, 2022.

[Bradlow, 1998]
Bradlow, E.
Encouragement designs: an approach to self-selected samples in an
experimental design.
Marketing Letters, 9:0 383–391, 1998.

[Chandak et al., 2023]
Chandak, Y., Shankar, S., Syrgkanis, V., and Brunskill, E.
Adaptive instrument design for indirect experiments.
arXiv preprint arXiv:2312.02438, 2023.

[Chernozhukov et al., 2018]
Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey,
W., and Robins, J.
Double/debiased machine learning for treatment and structural
parameters, 2018.

[Degenne et al., 2020]
Degenne, R., Ménard, P., Shang, X., and Valko, M.
Gamification of pure exploration for linear bandits.
In International Conference on Machine Learning, pages
2432–2442. PMLR, 2020.

[Della Vecchia and Basu, 2023]
Della Vecchia, R. and Basu, D.
Online instrumental variable regression: Regret analysis and bandit
arXiv preprint arXiv:2302.09357, 2023.

[Elbers, 2023]
Elbers, B.
Encouragement designs and instrumental variables for a/b testing,
August 2023.
Accessed: 2024-05-19.

[Engineering, 2023]
Engineering, S.
Encouragement designs and instrumental variables for a/b testing,
Accessed: December 26, 2023.

[Even-Dar et al., 2006]
Even-Dar, E., Mannor, S., Mansour, Y., and Mahadevan, S.
Action elimination and stopping conditions for the multi-armed bandit
and reinforcement learning problems.
Journal of machine learning research, 70 (6), 2006.

[Fiez et al., 2019]
Fiez, T., Jain, L., Jamieson, K. G., and Ratliff, L.
Sequential experimental design for transductive linear bandits.
Advances in Neural Information Processing Systems, 32, 2019.

[Gales et al., 2022]
Gales, S. B., Sethuraman, S., and Jun, K.-S.
Norm-agnostic linear bandits.
In International Conference on Artificial Intelligence and
Statistics, pages 73–91. PMLR, 2022.

[Gong and Zhang, 2022]
Gong, X. and Zhang, J.
Dual instrumental method for confounded kernelized bandits.
arXiv preprint arXiv:2209.03224, 2022.

[Greene, 2003]
Greene, W. H.
Econometric analysis.
Pearson Education India, 2003.

[Inoue and Solon, 2010]
Inoue, A. and Solon, G.
Two-sample instrumental variables estimators.
The Review of Economics and Statistics, 920
(3):0 557–561, 2010.

[Jun et al., 2021]
Jun, K.-S., Jain, L., Mason, B., and Nassif, H.
Improved confidence bounds for the linear logistic model and
applications to bandits.
In International Conference on Machine Learning, pages
5148–5157. PMLR, 2021.

[Kallus, 2018]
Kallus, N.
Instrument-armed bandits.
In Algorithmic Learning Theory, pages 529–546. PMLR, 2018.

[Katz-Samuels et al., 2020]
Katz-Samuels, J., Jain, L., Karnin, Z., and Jamieson, K. G.
An empirical process approach to the union bound: Practical
algorithms for combinatorial and linear bandits.
In Advances in Neural Information Processing Systems,
volume 33, pages 10371–10382, 2020.

[Kohavi et al., 2020]
Kohavi, R., Tang, D., and Xu, Y.
Trustworthy online controlled experiments: A practical guide to
a/b testing.
Cambridge University Press, 2020.

[Kveton et al., 2023]
Kveton, B., Liu, Y., Kruijssen, J. M., and Nie, Y.
Non-compliant bandits.
In Proceedings of the 32nd ACM International Conference on
Information and Knowledge Management, pages 1138–1147, 2023.

[Lattimore et al., 2016]
Lattimore, F., Lattimore, T., and Reid, M. D.
Causal bandits: Learning good interventions via causal inference.
Advances in Neural Information Processing Systems, 29, 2016.

[Lattimore and Szepesvári, 2020]
Lattimore, T. and Szepesvári, C.
Bandit algorithms.
Cambridge University Press, 2020.

[Li et al., 2023]
Li, Z., Jamieson, K., and Jain, L.
Optimal exploration is no harder than thompson sampling.
arXiv preprint arXiv:2310.06069, 2023.

[Lu et al., 2020]
Lu, Y., Meisami, A., Tewari, A., and Yan, W.
Regret analysis of bandit problems with causal background knowledge.
In Conference on Uncertainty in Artificial Intelligence, pages
141–150. PMLR, 2020.

[Malek et al., 2023]
Malek, A., Aglietti, V., and Chiappa, S.
Additive causal bandits with unknown graph.
In International Conference on Machine Learning, 2023.

[Mummalaneni et al., 2022]
Mummalaneni, S., Yoganarasimhan, H., and Pathak, V. V.
Producer and consumer engagement on social media platforms.
SSRN 4173537, 2022.

[Pearl, 2009]
Pearl, J.
Cambridge university press, 2009.

[Sen et al., 2017]
Sen, R., Shanmugam, K., Dimakis, A. G., and Shakkottai, S.
Identifying best interventions through online importance sampling.
In International Conference on Machine Learning, pages
3057–3066. PMLR, 2017.

[Soare et al., 2014]
Soare, M., Lazaric, A., and Munos, R.
Best-arm identification in linear bandits.
Advances in Neural Information Processing Systems, 27, 2014.

[Vershynin, 2018]
Vershynin, R.
High-dimensional probability: An introduction with applications
in data science, volume 47.
Cambridge university press, 2018.

[Weltz et al., 2023]
Weltz, J., Fiez, T., Volfovsky, A., Laber, E., Mason, B., Nassif, H., and Jain,
Experimental designs for heteroskedastic variance.
arXiv preprint arXiv:2310.04390, 2023.

[Zhang et al., 2022]
Zhang, J., Chen, Y., and Singh, A.
Causal bandits: Online decision-making in endogenous settings.
arXiv preprint arXiv:2211.08649, 2022.

§ BROADER IMPACTS

This work is algorithmic and not tied to a particular application that would have immediate negative impact.

§ RELATED WORKS

Our work is at the intersection of several parallel tracks of literature, pure exploration linear bandits, causal bandits, and econometrics.
The most relevant work on pure exploration in linear bandits is the  algorithm of [Fiez et al., 2019] which builds upon the early work of [Soare et al., 2014] which proposed experimental design combined with a phased elimination algorithm.  is nearly instance optimal for linear bandits in the non-confounded setting - that is it has a sample complexity which is within constant factors of the instance optimal lower bounds. Extensions of  to various noise models including logistic and heteroskedastic noise have been considered [Weltz et al., 2023, Jun et al., 2021]. Other algorithms for pure exploration linear bandits have been proposed - and we leave it for future work to extend the ideas of this paper to those settings [Li et al., 2023, Degenne et al., 2020].

Confounding in bandits was first considered in the regret minimization setting by the works of [Bareinboim et al., 2015]. The former work introduces the Multi-armed bandit with unobserved confounders (MABUC) problem. Similar to this work, they demonstrate that MAB techniques may be sub-optimal when latent unobserved confounders are present. Under a structural causal model formulation, they assume that they have additional to additional information in the form of “intuition” which can be observed at each time. The intuition captures a known distribution over the reward given latent factors and rewards. Based on this notion, they measure regret relative to the reward observed relative to playing in favor or against intuition.  They empirically demonstrate traditional bandit algorithms can have linear regret in this setting and provide an algorithm that effectively employs observed intuition.

A few works have considered the specific problem of compliance. The early work of [Kallus, 2018] assumes there is an additional unobserved latent class at each time that determines confounding. They provide novel notions of regret, relative to the instrument with the highest reward ($\arg\max Z^{\top}\Gamma^{\top}\theta$ in our notation), the highest treatment ($\arg\max_w w^{\top}\theta$), regret relative to the best latent class at each time, and regret on the set of “compliers”. They discuss the suitability of these various notions of regret, and discuss when sublinear regret is possible.
We remark that their approach is similar to ours in the sense that they assume a form of homogeneous effects across the population, and use an estimate of $\Gamma$.
Recently [Kveton et al., 2023] also consider the problem of compliance, however they don't take explicit non-confounding into account and assume an explicit parametric model that determines the non-compliance. This is analogous to the Heckman selection model considered in econometrics [Greene, 2003].

The recent works of [Della Vecchia and Basu, 2023, Zhang et al., 2022, Gong and Zhang, 2022] considered an online setting where at each time they observe a set $\{(x_t, z_t)\}$ where $x_t$ is the action of interest and $z_t$ is an associated instrument. If action $I_t$ is selected, the reward observed is $y_t = x_{I_t}^{\top}\theta_t + \epsilon_t$, where $x_t$ may be endogeneous. Similar to the standard linear bandit setting [Abbasi-Yadkori et al., 2011, Lattimore and Szepesvári, 2020], the goal is to minimize regret relative to the best action at each time. We remark that this setting is very different from ours. Effectively, we are choosing which instrument to select at each time to learn the best-performing treatment - in particular we can't choose a particular intervention. In their setting, they are choosing an intervention at each time and using the instrument purely for de-confounding the result. Experimental design for instruments to have more effective estimation has been considered by [Chandak et al., 2023].

In the causal bandit problem, an underlying causal graph

between a set of interventions and a reward value is assumed.
Actions correspond to intervening (i.e. a “do” operation [Pearl, 2009]) at one or more specific nodes in the causal graph and then observing the corresponding value at the reward node.
Causal bandits have been studied extensively in the regret setting [Lattimore et al., 2016, Lu et al., 2020, Bilodeau et al., 2022] and the pure exploration setting [Sen et al., 2017]. Though past works have allowed for unobserved confounders in the graph e.g. [Malek et al., 2023], their goal is to learn the best performing intervention, which in our setting would be  $\arg\max_{z\in \mathcal{Z}} z^{\top}\Gamma\theta$ instead of $w^{\ast}$.

Encouragement designs have been considered in many applications in online and offline settings. One of the earliest works on encouragement designs is [Bradlow, 1998], which considers the problem of using encouragements to determine the impact of coupons at a grocery store. More recent applications include  [Bakshy et al., 2014, Mummalaneni et al., 2022, Engineering, 2023] all in the context of online services and treatments that are required to be served to all users. Most of these works consider a small number of treatments and a heterogeneous treatment effect - hence are interested in LATE estimator.
As far as we are aware, we are the only work that considers adaptive encouragement design in the context of the model given in Equation <ref> and for multiple treatments.

§ ILLUSTRATIVE EXAMPLE

We now present an illustrative experiment that highlights the challenges of endogenous noise and the insufficiency of standard experimentation approaches used in the absence of confounding.

Instance Definition.
Toward connecting back to membership example in <ref>, consider that a service has $d$ membership options given by the set $\mathcal{A}=\{1,\dots, d\}$. Let the set $\mc{Z}=\cbr{e_1,\cdots,e_d}$ represent encouragements (incentives or advertisements) for the corresponding membership options given by $\mc{W}=\mc{X}=\cbr{e_1,\cdots,e_d}$. We consider a location model that assumes each user $t\in \mathbb{N}$ arriving online has an underlying unobserved one-dimensional preference $u_t \sim \mathcal{N}(0,\sigma^2_u)$. If an algorithm presents the user with encouragement $z_{I_t} = e_{I_t}$ for $I_t\in \mathcal{A}$, then the user selects into the membership level given by $J_t = \min_{j\in \mathcal{A}}|I_t + u_t - j|$ so that $x_t=e_{J_t}$. For a visual depiction, see Figure <ref>. This process captures a user being more likely to opt-in to membership levels that are closer to the encouragement that they were presented. The outcome is then given by $y_t = x_t^{\top}\theta + u_t$. This problem instance is a specific compliance instance. For this experiment, we take $d=6$, let $\theta=\begin{bmatrix}1& -0.95&0&0.45 &0.95& 0.99\end{bmatrix}$ and $\sigma_u^2=0.35$. Observe that the optimal evaluation vector is $w^{\ast}=e_1=\argmax_{w\in \mathcal{W}}w^{\top}\theta$.

We simulate a UCB strategy which maintains estimates of the average reward of each of the possible $d$ incentives, namely $\widehat{\mu}_{i,t} = \sum_{s=1}^t \mathbf{1}\{z_t = e_i\} y_t$ and then pulls the one with the highest upper confidence bound. This models current practice of using a bandit algorithm to select which incentive to show a user. Our results averaged over 100 simulations are in Figure <ref>. At each round we estimate the average reward of each level using an OLS estimator, i.e. $\widehat{\theta}^{OLS}_{i,t} = \sum_{s=1}^t \mathbf{1}\{x_t = e_i\} y_t / \sum_{s=1}^t \mathbf{1}\{x_t = e_i\}$, and check whether it matches the true value (denoted as UCB-OLS). We also consider an instrumental variable-estimator (see the next Section) which incorporates knowledge of $\Gamma$ similar to 2SLS to deconfound our estimate (UCB-IV).
As the plot demonstrates, UCB-OLS completely fails to identify $\theta_1 = \arg\max_{i\in d} \theta_i$ (this line is hard to see it is at 0) due to a biased estimate, whereas UCB-IV does better. However, UCB-IV methods seem to have a constant probability of error. To see why, note that the expected reward from pulling $z = e_{i}$ is $e_i^{\top}\Gamma\theta$. These values are plotted in orange in Figure <ref>. In particular, with some constant probability, UCB runs on the empirical rewards from pulling $z$'s zeroes in on arm 6, and as a result fails to give enough samples to learn that arm 1 is indeed the best. In contrast, our proposed method CPEG, Algorithm <ref> manages to find the best arm with significantly higher probability (the algorithm was run with $\delta = .1$) in the given time horizon.

[Depiction of the problem instance.]

[->] (0,0) – (7,0) node[right] ;

in 1,...,6
(,-0.3) – (,-1) node[below] $\x$;

[red] (4,0) circle (3pt) node[above right] $I_t$;

[black] (2,0) circle (3pt);
[black] (3,0) circle (3pt);
[black] (5,0) circle (3pt);
[black] (6,0) circle (3pt);

[orange] (1,0) circle (3pt) node[above left] $J_t$;
[blue] (1.49,0) circle (3pt) node[above right] $I_t+u_t$;
[blue, dashed] (4,0) to[out=90,in=90] (1.4,0);
[orange] (1.4,0) to[out=0,in=0] (1,0);

no markers,
axis lines=none,
axis line style=<->,
[thick, blue] exp(-x^2);
[above, blue, font=] at (axis cs: 0, 0.4) $u_t$;

[Heat-map of $\Gamma$.]

[$w^{\top}\theta$ and $w^{\top}\Gamma\theta$ for $w\in \{e_1,\dots, e_d\}$ ]
[$\mathbb{P}(\widehat{w}_t=w^{\ast})$ for the instance]
(a) A visual depiction of the problem instance from Section <ref>. The user is presented with encouragement $I_t\in \mathcal{A}$ and the user choice is given by $J_t$ where $J_t = \min_{j\in \mathcal{A}}|I_t + u_t - j|$ and $u_t \sim \mathcal{N}(0,\sigma^2_u)$. b) A heat-map showing the structural parameter $\Gamma$ for the problem instance from Section <ref>. (c) A bar chart showing $\mathbb{E}[y|x=w]=w^{\top}\theta$ and $\mathbb{E}[y|z=w]=z^{\top}\Gamma\theta$ for all $w\in \mathcal{W}$. This chart shows that the optimal evaluation vector is $w^{\ast}=e_1=\argmax_{w\in \mathcal{W}}\mathbb{E}[y|x=w]$, while $e_6=\argmax_{w\in \mathcal{W}}\mathbb{E}[y|z=w]$ and consequently estimation based on this quantity is problematic. (d)  The probability of identifying $w^{\ast}=e_1$ for a collection of algorithms on the  instance described in Section <ref>. Standard optimistic sampling approaches in combination with an ordinary least squares estimator leads to faulty inferences. Given an instrumental variable estimator, these experimental designs eventually give high probability identification but do so inefficiently compared to our proposed approach (see Section <ref>).

§ STANDARD 2SLS ESTIMATOR

Consider $\Psi = \widehat{\Gamma}_{2SLS} = (Z_T^{\top}Z_T)^{-1}Z_T^{\top}X_T$. In this setting, we recover the standard two-stage-least-squares () estimator,

\begin{align*}
\widehat{\theta}_{\texttt{2SLS}} = (X_T^{\top}Z_T(Z_T^{\top}Z_T)^{-1}Z_T^{\top}X_T)^{-1}X_T^{\top}(Z_T^{\top}Z_T)^{-1}Z_T^{\top}Y_T
\end{align*}

Note that the   estimator is a biased, but consistent estimator of the parameter $\theta$ [Angrist et al., 1996, Greene, 2003].

Note that in particular, the asymptotic variance of  is known to be $\sigma_{\epsilon}^2\|w\|_{(\widehat{\Gamma}_{\texttt{2SLS}}^{\top}Z_{T}^{\top}Z_{T}\widehat{\Gamma}_{\texttt{2SLS}} )^{-1}}$[Greene, 2003]. Recent work by [Della Vecchia and Basu, 2023] provides a confidence interval of the form $|w^{\top}(\widehat{\theta}_{\texttt{2SLS}} - \theta)| \leq O(d\sigma_{\epsilon}^2\|w\|_{(\widehat{\Gamma}_{\texttt{2SLS}}^{\top}Z_{T}^{\top}Z_{T}\widehat{\Gamma}_{\texttt{2SLS}} )^{-1}}\sqrt{\log(T/\delta)})$. However,  it's unclear how to use the form of their confidence interval directly for experimental design due to the dependence of $\hat{\Gamma}_{2SLS}$ on the random quantity $X$. In addition, their work is not sufficiently general to handle the general forms of noise that we consider in <ref>.

§ A NON-INTERACTIVE LOWER BOUND

Due to the noise model from confounding and the dependence of the noise $\theta^{\top}\eta + \epsilon$, the instance-dependent lower bounds of [Fiez et al., 2019] do not immediately apply. In this section, we develop a lower bound tailored for the confounding setting.

Toward characterizing the optimal sample complexity, we develop a lower bound for a specific non-adaptive algorithm $\mathcal{A}$ that has access to the matrix $\Gamma$ governing the structural equation model. In particular, suppose that the non-adaptive algorithm $\mathcal{A}$ is allowed to select a sequence of $T$ measurements $\{z_{I_1},\dots z_{I_t} \dots, z_{I_T}\}$ to query prior to collecting any observations, where $I_t$ represents the index of the vector $z\in \mathcal{Z}$ chosen at time $t\in \{1, \dots, T\}$. Then, given the observations $\{y_{1}, \dots, \dots y_{t}, \dots, y_{T}\}$ generated by the environment, a candidate optimal vector $\widehat{w}\in \mathcal{W}$ is returned by the algorithm. We are interested in the necessary number of observations $T$ that must be collected in order to ensure $\mathbb{P}(\widehat{w}\neq w^{\ast})\leq \delta$ for some $\delta\in (0, 1)$. Thus, it is natural that the optimal non-adaptive algorithm $\mathcal{A}$ using the estimator $\widehat{\theta}_{\texttt{oracle}}$ forms a recommendation rule such that $\widehat{w}=\argmax_{w\in \mathcal{W}}w^{\top}\widehat{\theta}_{\texttt{oracle}}$.
We now state our lower bound result with respect to the non-adaptive oracle algorithm.

Consider a problem instance characterized by $\mathcal{W}\subset \mathbb{R}^d, \mathcal{Z}\subset \mathbb{R}^d$, $\Gamma \in \mathbb{R}^{d\times d}$, and $\theta \in \mathbb{R}^d$. Assume $\Gamma$ is known, $\theta$ is unknown, and the noise process is jointly Gaussian and defined by $\gamma := \begin{bmatrix}\eta & \epsilon \end{bmatrix}\sim \mathcal{N}(0, \Sigma)$ where $\Sigma\in \mathbb{R}^{(d+1)\times (d+1)}$ is an arbitrary correlation matrix. For $\delta \in (0, 0.05]$, if the non-adaptive oracle algorithm acquires $T\leq \sigma^2\rho^{\ast}\log(1/\delta)/2$ samples on the problem instance where $\sigma^2:=v^{\top}\Sigma v$ and $v:=\begin{bmatrix}\theta & 1\end{bmatrix}\in \mathbb{R}^{d+1}$, then $\mathbb{P}(\widehat{w}\neq w^{\ast})\geq \delta$.

There exists a problem instance characterized by $\mathcal{W}\subset \mathbb{R}^d, \mathcal{Z}\subset \mathbb{R}^d$, $\Gamma \in \mathbb{R}^{d\times d}$, and $\theta \in \mathbb{R}^d$ with a noise process satisfying Assumption <ref> such that if the non-adaptive oracle algorithm acquires $T\leq \max\{d\|\theta\|_2^2, \sqrt{d}\|\theta\|_2\}\rho^{\ast}\log(1/\delta)/2$ samples, then $\mathbb{P}(\widehat{w}\neq w^{\ast})\geq \delta$ for $\delta\in (0, 0.05]$.

The proof of Theorem <ref> is in Appendix <ref>. Notably, the result is reminiscent of lower bounds for the standard pure exploration transductive linear bandit problem without confounding [Fiez et al., 2019, Katz-Samuels et al., 2020] when given the measurement set $\{\Gamma^{\top}z\}_{z\in \mathcal{W}}$, evaluation set $\mathcal{Z}$, and parameter $\theta$.

Notably, the upper bounds for our algorithms nearly match the lower bound of Theorem <ref>. However, it is interesting to observe that the sample complexity incurs an additional factor of $d^2$ relative to the standard transductive linear bandit problem in the most general setting where entries of $\eta$ are sub-Gaussian, but not necessarily independent nor bounded. This illustrates the additional difficulty that arises from confounding. We point out that this is not likely to be a tight lower bound. In particular, it is a lower bound with respect to a non-adaptive algorithm that uses the particular choice of estimator. We leave improved lower bounds to future work.

§ EXPERIMENTS

We now present experiments a collection of experiments on  problem instances. The experiments demonstrate that our approach produces efficient designs for inference and estimation.

§.§ Comparison Algorithms

The baselines that our approaches are compared with are discussed below. We run experiments both when $\Gamma$ is known and when $\Gamma$ is fully unknown.

§.§.§ Known $\Gamma$.

To standardize the experiments, the baselines considered run in rounds mirroring the structure of Algorithm <ref>. Specifically, in round $k\in \mathbb{N}$ a sampling algorithm selects a design $\lambda_k\in \Delta(\mathcal{Z})$, collects  $N_k$ samples from the design, and forms a $\Psi-\texttt{IV}$ estimate of $\theta$ with $\Psi=\Gamma$ that is combined with a confidence interval (Lemma <ref> to either eliminate evaluation vectors or validate a stopping condition. The number of samples $N_k$ taken in round $k\in \mathbb{N}$ by any of the algorithms is given by
$N_k =\lceil 2(1+\omega)\zeta_k^{-2}\rho(\mc{W}_k)L_\nu\log(4k^2|\mc{W}|/\delta)\rceil \vee r(\omega)$
where $\rho(\mc{W}_k) =\max_{w,w'\in\mc{W}_k}f(w,w',\Gamma,\lambda)$ for design $\lambda_k\in \Delta(\mathcal{Z})$ and an active set of evaluation vectors $\mc{W}_k$.
The round sample count guarantees that given any experimental design, all vectors $w\in \mathcal{W}$ such that  $(w^{\ast}-w)^{\top}\theta> 2\cdot 2^{-k}$ can be determined to be suboptimal by the end of round $k$.
The sampling methods we consider are now described.

* Static Oracle. This design selects $\lambda_k=\argmin_{\lambda \in \Delta(\mc{Z})}\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}}f(w^{\ast},w',\Gamma,\lambda)$.

* Static XY-Optimal. This design selects $\lambda_k=\argmin_{\lambda \in \Delta(\mc{Z})}\max_{w, w'\in \mathcal{W}}f(w,w',\Gamma,\lambda)$.

* Static Uniform. This design selects $\lambda_{k, z} =1/|\mathcal{Z}|\ \forall \ z\in \mathcal{Z}$.

* Adaptive Uniform (SE). This design selects $\lambda_{k, w} =1/|\mathcal{W}_k| \ \forall \ w\in \mathcal{W}_k$. Note that this algorithm is effectively an adaption of action-elimination [Even-Dar et al., 2006] .

The static designs are independent of the round and simply terminate when all evaluation vectors can be eliminated except for a recommended optimal vector $\widehat{w}$.

§.§.§ Unknown $\Gamma$.
For this set of experiments, we compare Algorithm <ref> against a collection of variations of the sampling procedures. Specifically, we compare against methods that either replace only the experimental design for estimating $\Gamma$, or only the experimental design for estimating $\theta$, or both with uniform sampling. We label the approaches as $N-N$, where $N$ represents the sampling approaches (XY or uniform) for $\Gamma$ and $\theta$ respectively. Moreover, to make our approach more practical, we modify the algorithm so that $\overline{\log}\rbr{Z_{T},\delta}=4d+\log(1/\delta)$. We find that even with this modification to the algorithm, correctness is maintained empirically.

[Experiment 1: Known $\Gamma$]
[Experiment 2: Known $\Gamma$]

[Experiment 1: Unknown $\Gamma$]
[Experiment 2: Unknown $\Gamma$]
Sample complexity for algorithms on  problems. Our approach is consistently competitive across the experiments.

§.§ Experiment 1: Jump-Around Instance

We first return back to the location model of Section <ref>. Recall that $\mc{Z}=\mc{W}=\mc{X}=\cbr{e_1,\cdots,e_d}$. For this experiment, we take $d=6$, let $\theta=\begin{bmatrix}1& -0.95&0.45&0.45 &0.95& 0.45\end{bmatrix}$ and $\sigma_u^2=0.275$. The results of the experiment are shown in Figure <ref> for the case of known $\Gamma$. We see that Algorithm <ref> performs much better than the baselines and nearly matches the oracle design. Delving into the approach, it is able to quickly eliminate all but $w_1$ and $w_{d-1}$ and then puts more mass on $z_1$ and $z_{d-1}$ to reduce the uncertainty on $w_1$ and $w_{d-1}$. For the case of unknown $\Gamma$, the results are shown in Figure <ref>, where $\theta_5$ is reduced to $0.9$ so that all approaches could finish.

§.§ Experiment 2: Interpolation Instance

Let $\mc{Z}=\mc{W}=\mc{X}=\cbr{e_1,\cdots,e_d}$ define the measurement, evaluation, and observation sets. Consider that $\Gamma:=\frac{(1-\varepsilon)}{d}1_d1_d^{\top} + \varepsilon I_{d}$ for a parameter $\varepsilon\in (0, 1)$ where $1_d$ is a $d$-dimensional vector of 1's and $I_d$ is the $d$-dimensional identity matrix. For this experiment, we take $d=4$ and let $\theta=\begin{bmatrix}0.5& 0.583&0.67&0.75\end{bmatrix}$. As in all compliance instances, $\eta_t=x_t-\Gamma^{\top}z_t$, and in this simulation $\eta_t=0.4 \eta_t^{\top}v_t$, where $v_t=\bar{v}_t/\|\bar{v}_t\|_2$ and $\bar{v}_t\sim \mathcal{N}(0, I_d)$. The results of the experiment are shown in Figure <ref> for $\varepsilon\in \{1, 0.9, 0.8, 0.7\}$ with $\Gamma$ known. Note that Static-XY and Uniform overlap, and SE and CPEG overlap. We see that Algorithm <ref> and the adaptive uniform strategy perform similarly and near optimally. This is to be expected since the most efficient way to gather observations for treatments is to encourage that treatment, given that if the encouragement is not followed each of the alternatives is equally likely and provides no additional information of interest. Moreover, as discussed earlier, the problem gets more challenging as $\Gamma \rightarrow 1_d1_d^{\top}/d$. Note that the identity matrix could be replaced with a permutation matrix, in which case uniform sampling with elimination becomes highly suboptimal. The results for the case of $\Gamma$ unknown are shown in Figure <ref> with $\varepsilon=0.99$. This shows the value that comes from the experimental design for estimating both $\Gamma$ and $\theta$.

§ PROOFS OF THE LOWER BOUND

§.§ Proof of Theorem <ref>

*app-thm:lower_bndTheorem <ref>
\beginapp-thm:lower_bnd
Consider a problem instance characterized by $\mathcal{W}\subset \mathbb{R}^d, \mathcal{Z}\subset \mathbb{R}^d$, $\Gamma \in \mathbb{R}^{d\times d}$, and $\theta \in \mathbb{R}^d$. Assume $\Gamma$ is known, $\theta$ is unknown, and the noise process is jointly Gaussian and defined by $\gamma := \begin{bmatrix}\eta & \epsilon \end{bmatrix}\sim \mathcal{N}(0, \Sigma)$ where $\Sigma\in \mathbb{R}^{(d+1)\times (d+1)}$ is an arbitrary correlation matrix. For $\delta \in (0, 0.05]$, if the non-adaptive oracle algorithm acquires $T\leq \sigma^2\rho^{\ast}\log(1/\delta)/2$ samples on the problem instance where $\sigma^2:=v^{\top}\Sigma v$ and $v:=\begin{bmatrix}\theta & 1\end{bmatrix}\in \mathbb{R}^{d+1}$, then $\mathbb{P}(\widehat{w}\neq w^{\ast})\geq \delta$.
\endapp-thm:lower_bnd

We begin by recalling the framework of the non-adaptive oracle algorithm and discussing the properties of its estimator for the noise structure described in the statement of the result.

Non-Adaptive Oracle and Instance Definition.
The non-adaptive oracle algorithm $\mathcal{A}$ selects $T$ measurements to query prior to collecting any data. Let $I_t$ represent the index of the vector $z\in \mathcal{Z}$ chosen at time $t\in \{1, \dots, T\}$.
The noise process for the instance under consideration is assumed to be jointly Gaussian and defined by $\gamma_t := \begin{bmatrix}\eta_t & \epsilon_t \end{bmatrix}\sim \mathcal{N}(0, \Sigma)$ where $\Sigma\in \mathbb{R}^{(d+1)\times (d+1)}$ is an arbitrary positive semidefinite matrix. Defining $\bar{x}_{I_t}:=\Gamma^{\top}z_{I_t}$, $v:=\begin{bmatrix}\theta & 1\end{bmatrix}\in \mathbb{R}^{d+1}$ and $\nu_t:=v^{\top}\gamma_t$, the feedback model can be described as follows:

\begin{align*}
y_t &=x_t\theta + \epsilon_t \\
&=(\Gamma^{\top}z_{I_t})^{\top}\theta+\eta_t^{\top}\theta + \epsilon_t \\
&=:(\Gamma^{\top}z_{I_t})^{\top}\theta+v^{\top}\gamma_t \\
\end{align*}

Observe that the noise is independent and identically distributed as $\nu_t \sim \mathcal{N}(0, \sigma^2)$ where $\sigma^2:=v^{\top}\Sigma v$ since $\gamma_t \sim \mathcal{N}(0, \Sigma)$. Moreover, the noise process is exogeneous with $\mathbb{E}[\nu_t|\bar{x}_{I_t}] =0$ since $\bar{x}_{I_t}$ is deterministic given the index choice $I_t$.

Let  $\{z_{I_t}\}_{t=1}^T, \{\bar{x}_{I_t}\}_{t=1}^T$, and $\{y_t\}_{t=1}^T$ denote the observations collected by the non-adaptive oracle algorithm $\mathcal{A}$ and define $Z_{T}\in \mathbb{R}^{T\times d}$, $\bar{X}_{T}\in \mathbb{R}^{T\times d}$, and $Y_{T}\in \mathbb{R}^{T}$ to contain the respective stacked observations. Algorithm $\mathcal{A}$ obtains an estimate
$\widehat{\theta}_{\texttt{oracle}}$ by minimizing the sum of squares as follows:

\begin{equation*}
\textstyle \widehat{\theta}_{\texttt{oracle}}:=\argmin_{\widehat{\theta}\in \mathbb{R}^d}\sum_{t=1}^{T}(y_t-\bar{x}_{I_t}\widehat{\theta})^2=( \bar{X}_T^{\top} \bar{X}_T)^{-1} \bar{X}_T^\top Y_T=( Z_T^{\top} Z_T\Gamma)^{-1} Z_T^\top Y_T.
\end{equation*}

Given $\widehat{\theta}_{\texttt{oracle}}$, the non-adaptive oracle algorithms $\mathcal{A}$ returns a recommendation defined by $\widehat{w}=\argmax_{w\in \mathcal{W}}w^{\top}\widehat{\theta}_{\texttt{oracle}}$.
Note that since $\widehat{\theta}_{\texttt{oracle}}$ is obtained by least squares with exogeneous, independent and identically distributed mean-zero Gaussian noise, it is straightforward to verify the estimator is distributed as

\begin{equation}
\widehat{\theta}_{\texttt{Oracle}}-\theta \sim \mathcal{N}\Big(0, \sigma^2\cdot \bar{A}(Z_T, \Gamma)^{-1}\Big),
\label{eq:dist_normal}
\end{equation}

\begin{equation*}
\bar{A}(Z_T, \Gamma):=\Big(\sum_{t=1}^T\Gamma^{\top}z_{I_t}z_{I_t}\Gamma\Big)=\bar{X}_T^{\top} \bar{X}_T= \Gamma^{\top}Z_T^{\top} Z_T\Gamma.
\end{equation*}

Proof by Contradiction. To begin, recall that

\begin{equation*}
\rho^{\ast} := \min_{\lambda \in \Delta(\mc{Z})}\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}} \frac{\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2}.
\end{equation*}

Suppose for the sake of contradiction that the number of samples collected by the non-adaptive oracle algorithm $\mathcal{A}$ is $T\leq \sigma^2\rho^{\ast}\log(1/\delta)/2$ and $\mathbb{P}(\widehat{w}\neq w^{\ast})<\delta$ for $\delta \in (0, 0.05]$. To reach a contradiction, we analyze the distribution of $(w-w^{\ast})^{\top}\widehat{\theta}$ for some $w\neq w^{\ast}$ and show that with probability at least $\delta$ it is positive. We remark that this proof follows similar techniques to that of the proof of Theorem 3 of [Katz-Samuels et al., 2020].

Let $\overline{\lambda} \in \Delta(\mathcal{Z})$ represent the empirical sampling distribution of the algorithm $\mathcal{A}$, which is defined such that $\overline{\lambda}_z=\frac{1}{T}\sum_{t=1}^T\mathbbm{1}\{z_{I_t} = z\}$ for each $z\in \mathcal{Z}$. Moreover, define

\begin{equation*}
\rho^{\ast}(\overline{\lambda}) := \max_{w\in \mathcal{W}\setminus \{w^{\ast}\}} \frac{\|w^{\ast}-w\|^2_{A(\overline{\lambda}, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2}\quad \text{and} \quad \widetilde{w}\in\argmax_{w\in \mathcal{W}\setminus \{w^{\ast}\}} \frac{\|w^{\ast}-w\|^2_{A(\overline{\lambda}, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2}.
\end{equation*}

Note that $\rho^{\ast}(\overline{\lambda}) \geq \rho^{\ast}$ and observe by definition,

\begin{equation*}
\frac{A(\overline{\lambda}, \Gamma)^{-1}}{T}:=\frac{(\sum_{z\in \mathcal{Z}}\overline{\lambda}_z \Gamma^{\top} zz^{\top}\Gamma)^{-1}}{T}=\Big(\sum_{t=1}^T\Gamma^{\top}z_{I_t}z_{I_t}\Gamma\Big)^{-1}:=\bar{A}(Z_T, \Gamma)^{-1}.
\end{equation*}

Thus, by Eq. (<ref>),

\begin{equation*}
\widehat{\theta}_{\texttt{Oracle}}-\theta \sim \mathcal{N}\Big(0, \sigma^2\cdot \frac{A(\overline{\lambda}, \Gamma)^{-1}}{T}\Big),
\end{equation*}

\begin{equation*}
\frac{(\widetilde{w}-w^{\ast})^{\top}(\widehat{\theta}_{\texttt{Oracle}}-\theta)}{\langle w^{\ast}-\widetilde{w}, \theta\rangle}\sim \mathcal{N}\Big(0, \sigma^2\cdot \frac{\|w^{\ast}-\widetilde{w}\|_{A(\overline{\lambda}, \Gamma)^{-1}}^2}{T\cdot \langle w^{\ast}-\widetilde{w}, \theta\rangle^2}\Big).
\end{equation*}

Furthermore, by the definition of $\rho^{\ast}(\overline{\lambda})$,  the assumption $T\leq \sigma^2\rho^{\ast}\log(1/\delta)/2$, and the fact $\rho^{\ast}(\overline{\lambda}) \geq \rho^{\ast}$, we obtain

\begin{equation}
\mathbb{V}\Big(\frac{(\widetilde{w}-w^{\ast})^{\top}(\widehat{\theta}_{\texttt{Oracle}}-\theta)}{\langle w^{\ast}-\widetilde{w}, \theta\rangle}\Big)= \sigma^2\cdot \frac{\|w^{\ast}-\widetilde{w}\|_{A(\overline{\lambda}, \Gamma)^{-1}}^2}{T\cdot \langle w^{\ast}-\widetilde{w}, \theta\rangle^2}:=\sigma^2\cdot \frac{\rho^{\ast}(\overline{\lambda})}{T}\geq  \frac{2}{\log(1/\delta)}.
\label{eq:var_ordering}
\end{equation}

Now, consider a random variable $W\sim \mc{N}(0, 1)$. Proposition 2.1.2 of Vershynin, 2018 gives an anti-concentration result showing that for all $\zeta >0$,

\begin{equation}
\mathbb{P}(W \geq \zeta)\geq \Big(\frac{1}{\zeta}-\frac{1}{\zeta^3}\Big)\frac{1}{\sqrt{2\pi}}e^{-\zeta^2/2}.
\label{eq:anti_concentration}
\end{equation}

We apply this result to the quantity $(\widetilde{w}-w^{\ast})^{\top}(\widehat{\theta}_{\texttt{Oracle}}-\theta)/\langle w^{\ast}-\widetilde{w}, \theta\rangle$ to conclude that $(\widetilde{w}-w^{\ast})^{\top}\widehat{\theta}_{\texttt{Oracle}}>0$ with probability at least $\delta$. Toward doing so,
let $c\in (1, 1.15]$ be a constant and define

\begin{equation*}
\widetilde{W}\sim \mc{N}\Big(0, \frac{2}{\log(1/\delta)}\Big)\quad \text{and} \quad W:=\frac{\widetilde{W}}{\sqrt{2/\log(1/\delta)}} \quad \text{and} \quad \gamma:=\frac{c}{\sqrt{2/\log(1/\delta)}}.
\end{equation*}

Observe that $W\sim \mathcal{N}(0, 1)$. The following analysis holds for $\delta \in (0, 0.05]$ given that $c\in (1, 1.15]$ as assumed:

\begin{align*}
\mathbb{P}\Big(\frac{(\widetilde{w}-w^{\ast})^{\top}(\widehat{\theta}_{\texttt{Oracle}}-\theta)}{\langle w^{\ast}-\widetilde{w}, \theta\rangle}\geq c\Big)
&\geq  \mathbb{P}(\widetilde{W}\geq c) \tag{By Eq.~\ref{eq:var_ordering}}\\
&=  \mathbb{P}\Big(\frac{\widetilde{W}}{\sqrt{2/\log(1/\delta)}}\geq \frac{c}{\sqrt{2/\log(1/\delta)}}\Big) \\
&:=  \mathbb{P}(W\geq \gamma) \\
&\geq \Big(\frac{1}{\gamma}-\frac{1}{\gamma^3}\Big)\frac{1}{\sqrt{2\pi}}e^{-\gamma^2/2} \tag{Proposition 2.1.2~Vershynin, 2018} \\
&= \Big(\frac{\sqrt{2}}{c\sqrt{\log(1/\delta)}}-\frac{\sqrt{2}^3}{c^3\sqrt{\log(1/\delta)}^3}\Big)\frac{1}{\sqrt{2\pi}}\delta^{c^2/4}\\
&\geq  \delta.
\end{align*}

The final inequality can be verified computationally.
Thus, with probability at least $\delta$ for $\delta \in (0, 0.05]$, we obtain

\begin{align*}
&\geq c(w^{\ast}-\widetilde{w})^{\top}\theta + (\widetilde{w}-w^{\ast})^{\top}\theta \\
&= c(w^{\ast}-\widetilde{w})^{\top}\theta - (w^{\ast}-\widetilde{w})^{\top}\theta \\
&= (c-1)(w^{\ast}-\widetilde{w})^{\top}\theta \\
&> 0.
\end{align*}

Observe that the final inequality holds since $c>1$ and $(w^{\ast}-\widetilde{w})^{\top}\theta> 0$ by definition.

This result directly implies that with probability at least $\delta$ for $\delta \in (0, 0.05]$, the vector $\widehat{w}$ returned by algorithm $\mathcal{A}$ is not $w^{\ast}$. This is a contradiction,  so we conclude that if the non-adaptive oracle algorithm $\mathcal{A}$ acquires $T\leq \sigma^2\rho^{\ast}\log(1/\delta)/2$ samples, then $\mathbb{P}(\widehat{w}\neq w^{\ast})\geq \delta$ for $\delta \in (0, 0.05]$.

§.§ Proof of Corollary <ref>

*app-cor:lower_bnd_corCorollary <ref>
\beginapp-cor:lower_bnd_cor
There exists a problem instance characterized by $\mathcal{W}\subset \mathbb{R}^d, \mathcal{Z}\subset \mathbb{R}^d$, $\Gamma \in \mathbb{R}^{d\times d}$, and $\theta \in \mathbb{R}^d$ with a noise process satisfying Assumption <ref> such that if the non-adaptive oracle algorithm acquires $T\leq \max\{d\|\theta\|_2^2, \sqrt{d}\|\theta\|_2\}\rho^{\ast}\log(1/\delta)/2$ samples, then $\mathbb{P}(\widehat{w}\neq w^{\ast})\geq \delta$ for $\delta\in (0, 0.05]$.
\endapp-cor:lower_bnd_cor

To begin, consider the specifications of Theorem <ref> and its result. That is, a problem an arbitrary instance characterized by $\mathcal{W}\subset \mathbb{R}^d, \mathcal{Z}\subset \mathbb{R}^d$, $\Gamma \in \mathbb{R}^{d\times d}$, and $\theta \in \mathbb{R}^d$ where the noise process is jointly Gaussian and defined by $\gamma := \begin{bmatrix}\eta & \epsilon \end{bmatrix}\sim \mathcal{N}(0, \Sigma)$ where $\Sigma\in \mathbb{R}^{(d+1)\times (d+1)}$ is an arbitrary correlation matrix. Observe that the noise process defined by $\gamma$ satisfies Assumption <ref>.
The result states that if the non-adaptive oracle algorithm acquires $T\leq \sigma^2\rho^{\ast}\log(1/\delta)/2$ samples on the problem instance where $\sigma^2:=v^{\top}\Sigma v$ and $v:=\begin{bmatrix}\theta & 1\end{bmatrix}\in \mathbb{R}^{d+1}$, then $\mathbb{P}(\widehat{w}\neq w^{\ast})\geq \delta$ for $\delta \in (0, 0.05]$. From this point, we show that there exists a parameter $\theta$ and correlation matrix $\Sigma$ such that $\sigma^2:=v^{\top}\Sigma v\geq \max\{d\|\theta\|_2^2, \sqrt{d}\|\theta\|_2\}$ in order to reach the stated conclusion.

Let $\zeta_{\eta_i, \epsilon}\in [-1, 1]$ denote the correlation between $\eta_i$ and $\epsilon$ for $i\in \{1, \dots, d\}$. Similarly, let $\zeta_{\eta_i, \eta_j}=\zeta_{\eta_j, \eta_i}\in [-1, 1]$ denote the correlation between $\eta_i$ and $\eta_j$ for $i\neq j\in \{1, \dots, d\}$. Note that the correlation of $\eta_i$ with itself for $i\in \{1, \dots, d\}$ is $\sigma_{\eta_i}^2=\zeta_{\eta_i, \eta_i}=1$ and similarly the correlation of $\epsilon$ with itself is $\sigma_{\epsilon}^2=\zeta_{\epsilon, \epsilon}=1$. The correlation matrix $\Sigma$ is then given by

\begin{equation*}
\Sigma =
\begin{bmatrix}
1& \zeta_{\eta_2, \eta_1} & \cdots & \zeta_{\eta_d, \eta_1} &\zeta_{\epsilon, \eta_{1}}  \\
\zeta_{\eta_1, \eta_2} & 1 & \cdots & \zeta_{\eta_{d}, \eta_{2}} & \zeta_{\epsilon, \eta_{2}} \\
\vdots & \vdots &\ddots & \vdots & \vdots \\
\zeta_{\eta_1, \eta_d} &  \zeta_{\eta_{2}, \eta_{d}} & \cdots & 1 & \zeta_{\epsilon, \eta_{d}}   \\
\zeta_{\eta_1, \epsilon} &  \zeta_{\eta_2, \epsilon} & \cdots & \zeta_{\eta_d, \epsilon} & 1
\end{bmatrix} :=
\begin{bmatrix}
\Sigma_{\theta} & \zeta_{\eta, \epsilon} \\
\zeta_{\eta, \epsilon}^{\top} & 1
\end{bmatrix},
\end{equation*}

\begin{equation*}
\Sigma_{\theta} =
\begin{bmatrix}
1& \zeta_{\eta_2, \eta_1} & \cdots & \zeta_{\eta_d, \eta_1}   \\
\zeta_{\eta_1, \eta_2} & 1 & \cdots & \zeta_{\eta_{d}, \eta_{2}}  \\
\vdots & \vdots &\ddots & \vdots  \\
\zeta_{\eta_1, \eta_d} &  \zeta_{\eta_{2}, \eta_{d}} & \cdots & 1
\end{bmatrix} \in [-1, 1]^{d\times d}
\quad \text{and} \quad  \zeta_{\eta, \epsilon}=
\begin{bmatrix}
\zeta_{\eta_1, \epsilon}  \\
\zeta_{\eta_2, \epsilon}   \\
\vdots   \\
\zeta_{\eta_{d}, \epsilon}
\end{bmatrix} \in [-1, 1]^{d}.
\end{equation*}

Moreover, we use the notation $\Sigma_{\theta, i}^{\top}\in \mathbb{R}^{d}$ to denote the $i$–th row of $\Sigma_{\theta}^{\top}$, or equivalently the $i$–th column of $\Sigma_{\theta}$, for any $i\in \{1, \dots, d\}$.

Lower Bounding Noise Variance. Given the above notation, we now work toward lower bounding $\sigma^2:=v^{\top}\Sigma v$.
Observe that by algebraic manipulations,

\begin{align*}
v^{\top}\Sigma v &:=
\begin{bmatrix}
\theta \\ 1
\end{bmatrix}^{\top}
\begin{bmatrix}
\Sigma_{\theta} & \zeta_{\eta, \epsilon} \\
\zeta_{\eta, \epsilon}^{\top} & 1
\end{bmatrix}
\begin{bmatrix}
\theta \\  1
\end{bmatrix} \\
\begin{bmatrix}
\theta^{\top}\Sigma_{\theta, 1}^{\top} +\zeta_{\eta_1, \epsilon} & \theta^{\top}\Sigma_{\theta, 2}^{\top} +\zeta_{\eta_2, \epsilon} & \cdots & \theta^{\top}\Sigma_{\theta, d}^{\top} +\zeta_{\eta_d, \epsilon} & \theta^{\top}\zeta_{\eta, \epsilon}+1
\end{bmatrix}
\begin{bmatrix}
\theta \\  1
\end{bmatrix} \\
&= \theta^{\top}\sum_{i=1}^d\Sigma_{\theta, i}^{\top}\theta_i + \sum_{i=1}^d\zeta_{\eta_{i}, \epsilon}\theta_i+\theta^{\top}\zeta_{\eta, \epsilon} + 1 \\
&= \theta^{\top}\Sigma_{\theta}\theta + 2\theta^{\top}\zeta_{\eta, \epsilon} + 1.
\end{align*}

Since $\Sigma_{\theta}$ is a real symmetric matrix, an eigendecomposition exists such that $\Sigma_{\theta} = Q\Lambda Q^{\top}$ where $\Lambda=\text{diag}(\lambda_1,\dots, \lambda_d)\in \mathbb{R}^{d\times d}$ is a diagonal matrix containing the eigenvalues of $\Sigma_{\theta}$ and $Q\in \mathbb{R}^{d\times d}$ is an orthogonal matrix with columns corresponding to the eigenvectors of $\Sigma_{\theta}$. Let $q_i:=Q_i^{\top}$ denote column $i$ of the matrix $Q$ for $i=\{1, \dots, d\}$, which is equivalently eigenvector $i$ of $\Sigma_{\theta}$ for $i=\{1, \dots, d\}$. Without loss of generality, assume that the eigenvectors are of unit length so that $\|q_i\|_2=1$ for all $i=\{1, \dots, d\}$.

Given this information, suppose that the parameter $\theta$ in the instance is equal to a scalar multiple of the eigenvector of $\Sigma_{\theta}$ corresponding to the maximum eigenvalue. Note this is equivalent to the statement that $\theta$ is equal to some scalar multiple of the column $q_{\ast}\in \mathbb{R}^d$ of the matrix $Q$ where $q_{\ast}$ is the eigenvector of $\Sigma_{\theta}$ corresponding to the maximum eigenvalue $\lambda^{\ast}$. Thus, we take $\theta = c\cdot q_{\ast}$ for some $c\in \mathbb{R}$ and observe that $\|\theta\|_2=c$. Toward quantifying the value of $v^{\top}\Sigma v$ for the problem instance, we begin by characterizing $\theta^{\top}\Sigma_{\theta}\theta$ for the choice of $\theta$. Consider the following analysis:

\begin{align*}
\theta^{\top}\Sigma_{\theta}\theta
&= \theta^{\top}Q\Lambda Q^{\top}\theta \\
&:= (cq_{\ast})^{\top}Q\Lambda Q^{\top}(cq_{\ast}) \tag{$\theta:=cq_{\ast}$} \\
\begin{bmatrix}
q_1&\cdots &q_{\ast}&\cdots & q_d
\end{bmatrix}
\Lambda
\begin{bmatrix}
q_1&\cdots &q_{\ast}&\cdots & q_d
\end{bmatrix}^{\top}q_{\ast} \\
\begin{bmatrix}
0&\cdots &\|q_{\ast}\|_2^2&\cdots & 0
\end{bmatrix}
\text{diag}(\lambda_1, \dots, \lambda^{\ast}, \dots, \lambda_1)
\begin{bmatrix} \tag{$q_{i}^{\top}q_j=0 \ \forall i\neq j$}
0&\cdots &\|q_{\ast}\|_2^2&\cdots & 0
\end{bmatrix}^{\top} \\
&= \|\theta\|_2\lambda^{\ast}.\tag{$\|q_{\ast}\|_2=1$}
\end{align*}

Thus, in general for this choice of $\theta$,

\begin{equation*}
v^{\top} \Sigma_{\theta}v=\|\theta\|_2^2\lambda^{\ast} + 2\theta^{\top}\zeta_{\eta, \epsilon} + 1.
\end{equation*}

To conclude, take $\Sigma_{\theta}:=\mathbf{1}_d\mathbf{1}_d^{\top}$ where $\mathbf{1}_d$ represents the $d$-dimensional vector of all ones. Since this is a rank-1 matrix, the maximum eigenvalue is $\lambda^{\ast}=\mathbf{1}_d^{\top}\mathbf{1}_d=d$ and the remainder of the eigenvalues are zero. Observe that $q_{\ast}=\mathbf{1}_d/\sqrt{d}$ is an eigenvector corresponding to the maximum eigenvalue since $\mathbf{1}_d\mathbf{1}_d^{\top} q_{\ast}=dq_{\ast}$. Thus,

\begin{align*}
v^{\top} \Sigma_{\theta}v &=\|\theta\|_2^2\lambda^{\ast} + 2\theta^{\top}\zeta_{\eta, \epsilon} + 1\\
&:= \|\theta\|_2^2\lambda^{\ast} + 2\|\theta\|_2\mathbf{1}_d^{\top}\mathbf{1}_d/\sqrt{d} + 1 \tag{$\theta:=cq_{\ast}:=\|\theta\|_2\mathbf{1}_d/\sqrt{d}=$ and $\zeta_{\eta, \epsilon}:=\mathbf{1}_d$} \\
&= d\|\theta\|_2^2 + 2\sqrt{d}\|\theta\|_2 + 1 \\
&\geq \max\{d\|\theta\|_2^2, \sqrt{d}\|\theta\|_2\}.
\end{align*}

This completes the proof since we have shown that there exists a parameter $\theta$ and correlation matrix $\Sigma$ such that $\sigma^2:=v^{\top}\Sigma v\geq \max\{d\|\theta\|_2^2, \sqrt{d}\|\theta\|_2\}$, which by Theorem <ref> allows us to make the stated conclusion.

§ PROOFS OF THE CONFIDENCE INTERVAL

§.§ Proof of <ref>

The first statement is an immediate consequence of <ref> and the second statement is proven in <ref>.

In the compliance model, the noise $\eta^\top\theta+\eps$ follows a $\rbr{8\vcn{\theta}_2^2+2}$-sub-Gaussian distribution.

In compliance, we have $z,x\in\cbr{e_1,\cdots,e_d}$, and

\begin{align*}
\eta=x-\rbr{\pp{e_1\mid z},\cdots, \pp{e_d\mid z}}^\top.
\end{align*}

Let us figure out the sub-Gaussian parameter of the random vector $\eta$.
Fix any unit vector $a$.
First, we have $\EE[\la \eta, a\ra] = 0$.
Second, we have

\begin{align*}
\abr{\eta^\top a}
&\le \vcn{\eta}_2\tag{Cauchy–Schwarz inequality}\\
&\le\vcn{  x-\rbr{\pp{e_1\mid  z},\cdots, \pp{e_d\mid  z}}^\top}_2\\
&\le\rbr{\vcn{x}_2+\vcn{\rbr{\pp{e_1\mid  z},\cdots, \pp{e_d\mid  z}}^\top}_2}\\
&\le\rbr{1+\vcn{\rbr{\pp{e_1\mid  z},\cdots, \pp{e_d\mid  z}}^\top}_1}\tag{$x\in\cbr{e_1,\cdots,e_d}$ and $\vcn{x}_2\le\vcn{x}_1, \forall x$}\\
\end{align*}

Thus, $\eta^\T a$ is bounded and zero-mean and thus $2^2$-sub-Gaussian.
This implies that

\begin{align*}
\forall \beta, \max_{a: \|a\|\le 1}\EE[\exp(\beta \la  \eta,a \ra)] \le \exp(\frac{\beta^2 2^2 }{2})~.
\end{align*}

and thus $\eta$ is a $2^2$-sub-Gaussan random vector.
Then, $\eta^\T\th$ is $(2 \|\th\|)^2$-sub-Gaussian.

Using <ref>, we have that $\eta^\T\th + \eps$ is $2(4\|\th\|^2 + 1) $-sub-Gaussian.

Let $A$ and $B$ random variables that are each $\sigma^2_A$- and $\sigma^2_B$-sub-Gaussian but are correlated.
Then, $A+B$ is $2(\sigma_A^2 + \sigma_B^2)$-sub-Gaussian.

By definition of sub-Gaussian, we have for any $\gamma\in\mathbb{R}$,

\begin{align*}
\ee{\ep{\gamma(A+B)}}
&= \ee{\ep{\gamma A} \ep{\gamma B}}
\\&\le \sqrt{\ee{\ep{2\gamma A}}} \sqrt{\ee{\ep{2\gamma B}}} \tag{Cauchy-Schwarz}
\\&\le \sqrt{\exp\left(2\gamma^2 \sigma_A^2\right)} \sqrt{\exp\left(2\gamma^2 \sigma_B^2\right)}
\\&\le \exp\left(2\gamma^2 (\sigma_A^2+\sigma_B^2)\right).
\end{align*}

§.§ Proof of <ref>

*app-thm:oracle 2SLSLemma <ref>
\beginapp-thm:oracle 2SLS
Suppose that $T$ observations are collected non-adaptively from the structural equation model in Eqs. (<ref>) and $\Gamma\in \mathbb{R}^{d\times d}$ is known. Then, with probability at least $1-\delta$ for $\delta\in (0, 1)$ and $w \in \mathbb{R}^d$,

\begin{align*}
&|w^{\top}(\widehat{\theta}_{\texttt{oracle}} - \theta)| \leq \sqrt{2\sigma_{\nu}^2\|w\|_{\bar{A}(Z_T, \Gamma)^{-1}} \log\left(2/\delta\right)}.
\end{align*}

where $\sigma_{\nu}^2$ is the sub-Gaussian parameter of the noise process $\nu:=\eta^{\top}\theta+\epsilon$ as characterized in Lemma <ref>.
\endapp-thm:oracle 2SLS

Given the knowledge of $\Gamma$, we have the oracle 2SLS estimator

\begin{align*}
\widehat{\theta}_{\texttt{oracle}}=\rbr{\sum_{t=1}^{T} z_s\rbr{\Gamma^\top z_s}^\top}^{-1}\sum_{t=1}^{T} z_sy_t=\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}\sum_{t=1}^{T} z_sy_t.
\end{align*}

Note that

\begin{align*}
y_t=x_t^\top\theta+\eps_t=\rbr{\Gamma^\top z_t}^\top\theta+\eta_t^\top\theta+\eps_t.
\end{align*}

Denote $\nu_t:=\eta_t^\top\theta+\eps_t$. For any $w\in\mc{W}$, we have

\begin{align*}
\ibr{\widehat{\theta}_{\texttt{oracle}}-\theta}{w}&=\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}\sum_{t=1}^{T} z_sy_t-\theta}{w}\\
&=\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}\sum_{t=1}^{T} z_s\rbr{ z_s^\top\Gamma\theta+\nu_t}-\theta}{w}\\
&=\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma\theta+\sum_{t=1}^{T} z_s\nu_t}-\theta}{w}\\
&=\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}\sum_{t=1}^{T} z_s\nu_t}{w}\\
&=\sum_{q=1}^{t}\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1} z_q}{w}\nu_q.
\end{align*}

By Lemma <ref>, the noise $\nu_t$ is $\sigma_{\nu}^2$-sub-Gaussian, we have

\begin{align*}
\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1} z_q}{w}\nu_q
\end{align*}

is $\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1} z_q}{w}^2\sigma_{\nu}^2$-sub-Gaussian. Thus

\begin{align*}
\pp{\ibr{\widehat{\theta}_{\texttt{oracle}}-\theta}{w}\ge\sqrt{2\sum_{q=1}^{t}\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1} z_q}{w}^2\sigma_{\nu}^2\log\rbr{\fr{1}{\delta}}}}\le\delta.
\end{align*}

\begin{align}
\sum_{q=1}^{t}\ibr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1} z_q}{w}^2&=w^\top\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}\rbr{\sum_{q=1}^{t} z_q z_q^\top}\rbr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}}^\top w \nn\\
&=w^\top\Gamma^{-1}\rbr{\rbr{\sum_{t=1}^{T} z_s z_s^\top\Gamma}^{-1}}^\top w \nn\\
&=w^\top\rbr{\Gamma^\top\rbr{\sum_{t=1}^{T} z_s z_s^\top}\Gamma}^{-1}w. \label{ieq:concise norm}
\end{align}

\begin{align*}
\pp{\ibr{\widehat{\theta}_{\texttt{oracle}}-\theta}{w}\ge\sqrt{2w^\top\rbr{\Gamma^\top\rbr{\sum_{t=1}^{T} z_s z_s^\top}\Gamma}^{-1}w\sigma_{\nu}^2\log\rbr{\fr{1}{\delta}}}}\le\delta.
\end{align*}

We can further write it as

\begin{align*}
\pp{\ibr{\widehat{\theta}_{\texttt{oracle}}-\theta}{w}\ge\sqrt{2\vcn{w}^2_{\rbr{\Gamma^\top\rbr{\sum_{t=1}^{T} z_s z_s^\top}\Gamma}^{-1}}\sigma_{\nu}^2\log\rbr{\fr{1}{\delta}}}}\le\delta.
\end{align*}

By taking a union bound over, we have the confidence interval for the absolute value as in the statement of the lemma.

§.§ Proof of <ref>

*app-thm:pseudo 2SLSTheorem <ref>
\beginapp-thm:pseudo 2SLS
Suppose that $\widehat{\Gamma}$ is estimated through a design matrix $Z_{T_1}\in \mathbb{R}^{T_1\times d}$ and $\widehat{\theta}_{\texttt{P-2SLS}}$ is estimated through a design matrix $Z_{T_2}\in \mathbb{R}^{T_2\times d}$. Then, for any $w \in \mathcal{W}$, with probability at least $1-\delta$,

\begin{align*}
|w^{\top}(\widehat{\theta}_{\texttt{P-2SLS}} - \theta)| \leq\vcn{w}_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4}{\delta}}}+\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}_2\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta/4}}.
\end{align*}

where $\sigma_{\nu}^2$ is the sub-Gaussian parameter of the noise $\nu:=\eta^{\top}\theta+\epsilon$, $\sigma_{\eta}^2$ is the sub-Gaussian parameter of the noise $\eta$, and

\begin{align*}
\overline{\log}\rbr{Z_{T},\delta}:=8d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(Z_{T_1}^\top Z_{T_1}))} } + 16\ln\del{\fr{2\cd6^d}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(Z_{T_1}^\top Z_{T_1})} }  }
\end{align*}

\endapp-thm:pseudo 2SLS

For the pseudo 2SLS estimator, we have

\begin{align*}
\widehat{\theta}_{\texttt{P-2SLS}} - \theta&=\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}\sum_{t=1}^{T_2}z_{I_t}y_t-\theta\\
\end{align*}

For any $w\in\mc{W}$, we have

\begin{align}
\ibr{\widehat{\theta}_{\texttt{P-2SLS}}-\theta}{w}&=\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}\sum_{t=1}^{T_2}z_{I_t}\nu_t+\rbr{\widehat{\Gamma}^{-1}\Gamma-I}\theta}{w} \nn\\
&=\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}\sum_{t=1}^{T_2}z_{I_t}\nu_t}{w}+\ibr{\rbr{\widehat{\Gamma}^{-1}\Gamma-I}\theta}{w} \nn\\
&=\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}\nu_q+\ibr{\rbr{\widehat{\Gamma}^{-1}\Gamma-I}\theta}{w} \nn\\
&=\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}\nu_q+\ibr{\rbr{\widehat{\Gamma}^{-1}-\Gamma^{-1}}\Gamma\theta}{w}. \label{ieq:confidence interval breakdown}
\end{align}

We upper bound the first term and the second term separately. For the first term, by <ref>, we know $\nu_t$ is $\sigma_{\nu}^2$-subGaussian, we have $\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}\nu_q$ is $\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}^2\sigma_{\nu}^2$-subGaussian. Thus by the concentration inequality of subGaussian random variables, we have

\begin{align*}
\pp{\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}\nu_q\ge\sqrt{2\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}^2\sigma_{\nu}^2\log\rbr{\fr{2}{\delta}}}}\le\fr{\delta}{2}.
\end{align*}

By similar calculation as (<ref>), we have with probability at least $1-\fr{\delta}{2}$,

\begin{align}
\sum_{q=1}^{T_2}\ibr{\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top\widehat{\Gamma}}^{-1}z_{I_q}}{w}\nu_q\le&\sqrt{2w^\top\rbr{\widehat{\Gamma}^\top\rbr{\sum_{t=1}^{T_2}z_{I_t}z_{I_t}^\top}\widehat{\Gamma}}^{-1}w\sigma_{\nu}^2\log\rbr{\fr{2}{\delta}}}\nn\\
\le&\sqrt{2\vcn{w}^2_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\sigma_{\nu}^2\log\rbr{\fr{2}{\delta}}}. \label{ieq:confidence interval p2sls 1}
\end{align}

Thus with probability at least $1-\fr{\delta}{2}$,

\begin{align*}
\ibr{\widehat{\theta}_{\texttt{P-2SLS}}-\theta}{w}\le\vcn{w}_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{2}{\delta}}}+\ibr{\rbr{\widehat{\Gamma}^{-1}-\Gamma^{-1}}\Gamma\theta}{w}
\end{align*}

By Theorem <ref>, we have with probability at least $1-\fr{\delta}{2}$,

\begin{align}
\ibr{\rbr{\widehat{\Gamma}^{-1}-\Gamma^{-1}}\Gamma\theta}{w}\le\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta/2}}\label{ieq:confidence interval p2sls 2}
\end{align}

Combining (<ref>) and (<ref>), we have with probability at least $1-\delta$, for any $w\in\mc{W}$,

\begin{align*}
\ibr{\widehat{\theta}_{\texttt{P-2SLS}}-\theta}{w}\le\vcn{w}_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{2}{\delta}}}+\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta/2}}.
\end{align*}

By a union bound, we have the confidence interval for the absolute value of the inner product,

\begin{align*}
|w^{\top}(\widehat{\theta}_{\texttt{P-2SLS}} - \theta)| \leq\vcn{w}_{\bar{A}(Z_{T_2}, \widehat{\Gamma})^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4}{\delta}}}+\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta/4}}.
\end{align*}

Suppose that the least square estimator $\widehat{\Gamma}$ is estimated through a design matrix $Z_{T_1}\in \mathbb{R}^{T_1\times d}$, then it satisfies, with probability at least $1-\delta$,

\begin{align*}
\ibr{\rbr{\widehat{\Gamma}^{-1}-\Gamma^{-1}}\Gamma\theta}{w}\le\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta}}
\end{align*}

for any $w\in\mc{W}$.

Define $V=Z_{T_1}^\top Z_{T_1}$, $S\in\mathbb{R}^{T_1\times d}$ as the matrix with $i$-th row being $\eta_i^\top$, the stacked noise, i.e., the data collection process of the design matrix $Z_{T_1}$ is $X=Z_{T_1}\Gamma+S$. Then we have

\begin{align*}
\ibr{\rbr{\widehat{\Gamma}^{-1}-\Gamma^{-1}}\Gamma\theta}{w}=&w^\top\rbr{\widehat{\Gamma}^{-1}-\Gamma^{-1}}\Gamma\theta\\
=&w^\top\rbr{\Gamma+V^{-1} Z_{T_1}^\top S}^{-1}V^{-1} Z_{T_1}^\top S\Gamma^{-1}\Gamma\theta \tag{\cref{lem:estimateGamma inverse}}\\
=&w^\top\widehat{\Gamma}^{-1}V^{-1/2}V^{-1/2}Z_{T_1}^\top S\theta\\
\le&\vcn{w^\top\widehat{\Gamma}^{-1}V^{-1/2}}\vcn{V^{-1/2}Z_{T_1}^\top S\theta}\\
\le&\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{V^{-1/2}Z_{T_1}^\top S}_{\mathrm{op}}\vcn{\theta}\\
\le&\vcn{w}_{\bar{A}(Z_{T_1}, \widehat{\Gamma})^{-1}}\vcn{\theta}\sqrt{\sigma_{\eta}^2\overline{\log}\rbr{Z_{T_1},\delta}},
\end{align*}

where the last inequality is due to Lemma <ref>.

Suppose we have $z_1,\ldots, z_T\in\RR^d$ and $\eta_1,\ldots,\eta_T\in\RR^d$ such that $\eta_T \mid z_1,\eta_1,\ldots,z_{T-1},\eta_{T-1},z_T$ is $\sigma_{\eta}^2$-sub-Gaussian vector (defined in <ref>).
Let $Z, S\in\mathbb{R}^{T\times d}$ be matrices whose $t$-th row is $z_t^\T$ and $\eta_t^\T$ respectively.
Suppose $\|z_t\| \le L_z,\forall t$.
Let $V= Z^\T Z$.
Then, $\forall \delta\in(0,1)$, we have, with probability at least $1-\delta$,

\begin{align*}
\vcn{V^{-1/2}Z_{T}^\top S}_{\text{op}}\le\sigma_{\eta}\sqrt{8d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(V))} } + 16\ln\del{\fr{2\cd6^d}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V)} }  } }.
\end{align*}

We abbreviate $\overline{\log}\rbr{Z_{T},\delta}:=8d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(V))} } + 16\ln\del{\fr{2\cd6^d}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V)} }  } $.

By the definition of operator norm, we have

\begin{align*}
\vcn{V^{-1/2} Z_{T}^\top  S}_{\text{op}} &= \sup_{\cbr{x\mid\vcn{x}_2=1}} \vcn{V^{-1/2} Z_{T}^\top Sx}_2\\
&= \sup_{\cbr{x\mid\vcn{x}_2=1}} \sqrt{x^{\top}S^{\top}Z_{T} V^{-1}Z_{T}^{\top}S x} \\
&= \sup_{\cbr{x\mid\vcn{x}_2=1}} \vcn{Z_{T}^{\top}S x}_{V^{-1}}.
\end{align*}

Considering a fixed $w\in\mathbb{R}^d$, by <ref>, we have with probability at least $1-\delta$,

\begin{align*}
\vcn{Z_{T}^{\top}S x}_{V^{-1}}&\le\sqrt{2}\sigma_{\eta}\sqrt{d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(V))} } + 2\ln\del{\fr{2}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V)} }  } }.
\end{align*}

By Lemma <ref> and a union bound, for the $\eps$-covering $\cC_\eps$, we have the following event happens with probability no more than $\delta$:

\begin{align*}
\mc{E}:=\cbr{\exists x\in \cC_\eps, \vcn{Z_{T}^{\top}S x}_{V^{-1}}\geq \sqrt{2}\sigma_{\eta}\sqrt{d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(V))} } + 2\ln\del{\fr{2\abr{\cC_\eps}}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V)} }  } }}.
\end{align*}

We abbreviate $\widehat{\log}\rbr{Z_{T},\delta}:=\sqrt{2}\sqrt{d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(V))} } + 2\ln\del{\fr{2\abr{\cC_\eps}}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V)} }  } }$
When $\mc{E}$ does not happen, we have

\begin{align*}
\vcn{V^{-1/2} Z^\top  S}_{\text{op}}&= \sup_{\cbr{x\mid\vcn{x}=1}} \vcn{Z_{T}^{\top}S x}_{V^{-1}} \\
&= \sup_{\cbr{x\mid\vcn{x}=1}} \min_{\cbr{y\mid y\in C_\eps}}\vcn{Z_{T}^{\top}S (x-y+y)}_{V^{-1}} \\
&\le \sup_{\cbr{x\mid\vcn{x}=1}} \min_{\cbr{y\mid y\in C_\eps}}\rbr{\vcn{Z_{T}^{\top}S (x-y)}_{V^{-1}}+\vcn{Z_{T}^{\top}S y}_{V^{-1}}} \\
&\le \sup_{\cbr{x\mid\vcn{x}=1}} \min_{\cbr{y\mid y\in C_\eps}}\rbr{\vcn{V^{-1/2}Z^\top S}_{\text{op}}\vcn{x-y}_2+\sigma_{\eta}\widehat{\log}\rbr{Z_{T},\delta}} \\
&\le \eps\vcn{V^{-1/2} Z_{T}^\top  S}_{\text{op}}+\sigma_{\eta}\widehat{\log}\rbr{Z_{T},\delta}.
\end{align*}

\begin{align*}
\vcn{V^{-1/2}Z_{T}^\top S}_{\text{op}}\le\fr{\sigma_{\eta}}{1-\eps}\widehat{\log}\rbr{Z_{T},\delta}.
\end{align*}

By choosing $\eps=\fr{1}{2}$ and <ref>, we have

\begin{align*}
\vcn{V^{-1/2}Z_{T}^\top S}_{\text{op}}\le\sigma_{\eta}\sqrt{8d \ln\del{1 + \fr{2TL_z^2}{d(2\wed \sigmin(V))} } + 16\ln\del{\fr{2\cd6^d}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V)} }  } }.
\end{align*}

[Self-Normalized Bound for Vector-Valued Martingales]
Suppose we have $z_1,\ldots, z_t\in\RR^d$ and $\eta_1,\ldots,\eta_t\in\RR^d$ such that $\eta_t \mid z_1,\eta_1,\ldots,z_{t-1},\eta_{t-1},z_t$ is $\sigma_{\eta}^2$-sub-Gaussian vector (defined in <ref>).
Let $Z, S\in\mathbb{R}^{t\times d}$ be matrices whose $s$-th row is $z_s^\T$ and $\eta_s^\T$ respectively. Suppose $\|z_s\| \le L_z,\forall s$.
Let $V_t= Z^\T Z$.
Then, $\forall b\in\mathbb{R}^t, \delta\in(0,1)$, we have, with probability at least $1-\delta$,

\begin{align*}
\forall t\ge1, \vcn{ Z^\top  Sb}_{V_t^{-1}}\le\sqrt{2}\vcn{b}_2\sigma_{\eta}\sqrt{d \ln\del{1 + \fr{2tL_z^2}{d(2\wed \sigmin(V_t))} } + 2\ln\del{\fr{2}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V_t)} }  } }.
\end{align*}

Since each row of $S$ is a $\sigma_{\eta}^2$-subGaussian vector, we have that $(Sb)_i / \|b\|$ is $\sigma_{\eta}^2$-subGaussian.
Using <ref> with $\eps_s = (Sb)_i / \|b\|$ completes the proof.

[Lemma 20.1]

There exists a set $\cC_\eps\subset \mathbb{R}^d$ with $\abr{\cC_\eps}\leq \rbr{\fr{3}{\eps}}^d$ such that for any $x\in\cbr{x\mid x\in\mathbb{R}^d, \vcn{x}_2=1}$, there exists a $y\in\cC_\eps$ such that $\vcn{x-y}_2\leq \eps$.

Let $z_1,z_2,\ldots \in \{z\in\RR^d: \|z\|_2 \le L_z\}$ and $\eps_1,\eps_2,\ldots \in \RR$ be random variables such that $\eps_k \mid z_1,\eps_1,\ldots,z_{t-1},\eps_{t-1},z_{t}$ is $\sig^2_\eps$-sub-Gaussian.
Let $V_t = \sum_{s=1}^t z_s z_s^\T$.

\begin{align*}
1-\dt \le \PP\del{\forall t\ge1, \norm{\sum_{s=1}^t z_s \eps_s}_{V_t^{-1}}
\le \sqrt{2} \sig_\eps\sqrt{d \ln\del{1 + \fr{2tL_z^2}{d(2\wed \sigmin(V_t))} } + 2\ln\del{\fr{2}{\dt}\cd \log^2_2\del{\fr{4}{2\wed \sigmin(V_t)} }  } } }
\end{align*}

Let $Z_t \in \RR^{t\times d}$ be the design matrix and define $\eps_t := (\eps_1,\ldots, \eps_t)^\T \in \RR^t$.
Let us omit the subscript $t$ from $Z_t$, $\eps_t$ and $V_t$.
Note that

\begin{align*}
\normz{Z^\T\eps}_{V^{-1}}
=\normz{Z^\T\eps}_{(\fr12 V + \fr12 V)^{-1}}
~\le\normz{Z^\T\eps}_{(\fr12 V + \fr12 \sigmin(V) I)^{-1}}
~\le\sqrt{2}\normz{Z^\T\eps}_{( V +  \sigmin(V) I)^{-1}}
\end{align*}

It remains to bound $\normz{Z^\T\eps}_{( V +  \sigmin(V) I)^{-1}}$.
We use union bound with the standard self-normalized inequality of Lattimore and Szepesvári, 2020[Theorem 20.4 and Note 20.2].
Specifically, let $\lam_k = 2^{-k+1}$ for $k\ge1$.

\begin{align*}
1-\dt \le \PP\del{\cE := \cbr{\forall k \in \NN_+, \forall t\ge1, \normz{Z^\T\eps}_{(V + \lam_k I)^{-1}} \le \sig_\eps\sqrt{d \ln\del{1 + \fr{tL^2}{d \lam_k} } + 2\ln\del{\fr{(\pi^2/6) k^2}{\dt} } } }  }
\end{align*}

Under the event $\cE$, we have two cases:

* $\sigmin(V) \ge 1$: choose $k=1$. Then,

\begin{align*}
\normz{Z^\T\eps}_{( V +  \sigmin(V) I)^{-1}}
&\le \normz{Z^\T\eps}_{( V + I)^{-1}}
\\&\le \sig_\eps\sqrt{d \ln\del{1 + \fr{tL^2}{d} } + 2\ln\del{\fr{(\pi^2/6)}{\dt} } }
\end{align*}

* $\sigmin(V) < 1$: choose $k = \lcl  \log_2(2\sigmin(V)^{-1}) \rcl \ge 2$, which is the $k$ that satisfies $\lam_k \le \sigmin(V) < \lam_{k-1}$. Then,

\begin{align*}
\normz{Z^\T\eps}_{( V +  \sigmin(V) I)^{-1}}
&\le \normz{Z^\T\eps}_{( V +  \lam_k I)^{-1}}
\\&\le \sig_\eps\sqrt{d \ln\del{1 + \fr{tL^2}{d \lam_k} } + 2\ln\del{\fr{(\pi^2/6) k^2}{\dt} } }
\\&\le \sig_\eps\sqrt{d \ln\del{1 + \fr{tL^2}{d  \sigmin(V)/2} } + 2\ln\del{\fr{(\pi^2/6) \lcl  \log_2(2\sigmin(V)^{-1}) \rcl^2}{\dt} } } \tag{$\lam_k = \fr12 \lam_{k-1} > \fr12 \sigmin(V)$; def'n of $k$ }~.
\end{align*}

Altogether, we have

\begin{align*}
\normz{Z^\T\eps}_{( V +  \sigmin(V) I)^{-1}}
\le \sig_\eps\sqrt{d \ln\del{1 + \fr{tL^2(1 \vee 2\sigmin^{-1}(V))}{d} } + 2\ln\del{\fr{(\pi^2/6) (1 \vee \lcl  \log_2(2\sigmin^{-1}(V)) \rcl^2)}{\dt} } } ~.
\end{align*}

We conclude the proof by simplifying the RHS.

§ PROOFS OF SAMPLE COMPLEXITY WHEN GIVEN

§.§ Proof of <ref>

Confounded pure exploration with known $\Gamma$

Input $\mc{Z}, \mc{W}, \Gamma, \delta, \epsilon, L_{\nu}\ge\sigma_{\nu}^2$
Initialize: $k=1, \hat{\mc{W}}_k=\mc{W}$
$f(w,w',\Gamma,\lambda):=\|w-w'\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}$
$N_k=\ur{2(1+\omega)\zeta_k^{-2}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}}\vee r(\omega)$
Pull arms in $Z_{N_k}$ and observe $Y_{N_k}$
Compute $\widehat{\theta}_{\texttt{oracle}}^k=\rbr{Z_{N_k}^\top Z_{N_k}\Gamma}^{-1}Z_{N_k}^\top Y_{N_k}$
$\mc{W}_{k+1}=\mc{W}_k\backslash\cbr{w\in\mc{W}_k\mid \exists w'\in\mc{W}_k,\text{s.t.},\ibr{w'-w}{\widehat{\theta}_{\texttt{oracle}}^k}>\zeta_k}$
Output: $\mc{W}_k$

*app-thm:complexity knownGammaTheorem <ref>
\beginapp-thm:complexity knownGamma
Algorithm <ref> is $\delta$-PAC for the  problem and terminates in at most $c(1+\omega)L_{\nu}\rho^{\ast}\log(1/\delta)+cr(\omega)$ samples,
where $c$ hides logarithmic factors of $\Delta$ and $\abr{\mc{W}}$.
\endapp-thm:complexity knownGamma

Part 1 $\delta$-PAC

By the confidence interval in <ref>, we have, with probability at least $1-\fr{\delta}{2 k^2\abr{\mc{W}}}$,

\begin{align*}
\abr{\ibr{\widehat{\theta}_{\texttt{oracle}}^k-\theta}{w-w^*}}&\le\vcn{w-w^*}_{\rbr{\Gamma^\top\rbr{\sum_{s=1}^{N_k }z_{I_t}z_{I_t}^\top}\Gamma}^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4 k^2\abr{\mc{W}}}{\delta}}}\\
&\le\fr{\sqrt{1+\omega}\vcn{w-w^*}_{\rbr{\Gamma^\top\lambda_z zz^{\top}\Gamma}^{-1}}}{\sqrt{N_k }}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4 k^2\abr{\mc{W}}}{\delta}}}\\
&\le\fr{\sqrt{1+\omega}\vcn{w-w^*}_{\rbr{\Gamma^\top\lambda_z zz^{\top}\Gamma}^{-1}}}{\sqrt{\ur{2(1+\omega)\zeta_k^{-2}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}}\vee r(\omega)}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4 k^2\abr{\mc{W}}}{\delta}}}\\
\end{align*}

Define a good event $\mc{E}_{k, w }$ for each $k$ and $\mc{W}_k$ as

\begin{align*}
\mc{E}_{k, w }=\cbr{\abr{\ibr{\widehat{\theta}_{\texttt{oracle}}^k-\theta}{w-w^*}}\le\zeta_k }.
\end{align*}

We claim that with probability at least $1-\delta$, the event $\bigcap_{k=1}^{\infty}\bigcap_{w\in\mc{W}_k}\mc{E}_{k,w}$ holds. It can be proved by a union bound.

\begin{align}
\pp{\rbr{\bigcap_{k=1}^{\infty}\bigcap_{w\in\mc{W}_k}\mc{E}_{w,k}}^c}\le&\sum_{k=1}^{\infty}\sum_{w\in\mc{W}_k}\pp{\mc{E}_{k,w}^c}\nn\\
\le&\sum_{k=1}^{\infty}\sum_{w\in\mc{W}_k}\fr{\delta}{2k^2\abr{\mc{W}}}\nn\\
\le&\fr{\delta}{2}\sum_{k=1}^{\infty}\fr{1}{k^2}\nn\\
\le&\delta, \label{ieq:good event}
\end{align}

where the last step is by the fact that $\sum_{k=1}^{\infty}\fr{1}{k^2}=\fr{\pi^2}{6}<2$. Under the event $\bigcap_{k=1}^{\infty}\bigcap_{w\in\mc{W}_k}\mc{E}_{w,k}$, to show that the best arm is never eliminated, it suffices to show that for any sub-optimal arm $w\in\mc{W}_k$,

\begin{align*}
\ibr{w-w^*}{\widehat{\theta}_{\texttt{oracle}}^k}&=\ibr{w-w^*}{\widehat{\theta}_{\texttt{oracle}}^k-\theta}+\ibr{w-w^*}{\theta}\\
\end{align*}

Thus the best arm $w^*$ never satisfies the elimination condition. Next we show that at the end of stage $k$, any suboptimal arm $w$ that satisfies

\begin{align*}
\ibr{\theta}{w^*-w}\ge&2\zeta_k
\end{align*}

is eliminated. To show this, we need to show that $w$ satisfies the elimination condition,

\begin{align*}
\max_{w'\in \mc{W}_k}\ibr{\widehat{\theta}_{\texttt{oracle}}^k}{w'-w}\ge&\ibr{\widehat{\theta}_{\texttt{oracle}}^k}{w^*-w}\\
\ge&-\zeta_k+2\zeta_k\\
\end{align*}

This implies that with probability at least $1-\delta$, $ w ^*$ always survives.

Part 2 sample complexity

Define $S_k =\cbr{ w \in\mc{W}\mid \ibr{ w ^*- w }{\theta}\le4\zeta_k }$. Thus with probability at least $1-\delta$, we have $\bigcap_{k}\cbr{\mc{W}_k \subseteq S_k }$. This implies the following is true with probability at least $1-\delta$ for all $k$

\begin{align*}
\rho(\mc{W}_k )&=\min_{\lambda\in\Delta(\mc{Z})}\max_{ w, w'\in \mc{W}_k }\vcn{w-w'}^2_{\rbr{\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma}^{-1}}\\
&\le\min_{\lambda\in\Delta(\mc{Z})}\max_{ w, w'\in S_k }\vcn{w-w'}^2_{\rbr{\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma}^{-1}}\\
&=\rho(S_k ).
\end{align*}

Define $\Delta$ to be the minimum gap between $ w ^*$ and any other $ w \in\mc{W}$, i.e., $\Delta:=\min_{ w \in\mc{W}\setminus\cbr{ w ^*}}\ibr{ w ^*- w }{\theta}$. Then for $k\ge\ur{\log\rbr{4\Delta^{-1}}}$, we have $S_k =\cbr{ w ^*}$ with probability at least $1-\delta$. The total sample complexity is the summation of the number of samples in each round,

\begin{align}
&\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}N_k \nn\\
=&\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}\rbr{\ur{4(1+\omega)\zeta_k^{-2}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}}\vee r(\omega)} \nn\\
\le&\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}\rbr{8(1+\omega)2^{2k}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}\vee r(\omega)} \nn\\
\le&\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}\rbr{8(1+\omega)2^{2k}\rho(\mc{S}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}\vee r(\omega)}. \label{ieq:sum of Nk}
\end{align}

On the other hand, we have

\begin{align}
\rho^{\ast} =& \min_{\lambda \in \Delta(\mc{Z})}\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}} \frac{\|w^{\ast}-w\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2} \nn\\
=&\min_{\lambda\in\Delta(\mc{Z})}\max_{k}\max_{ w \in S_k }\frac{\|w^{\ast}-w\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2} \nn\\
\ge&\fr{1}{\ur{\log\rbr{4\Delta^{-1}}}}\min_{\lambda\in\Delta(\mc{Z})}\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}\max_{ w \in S_k }\frac{\|w^{\ast}-w\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2} \nn\\
\ge&\fr{1}{16\ur{\log\rbr{4\Delta^{-1}}}}\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}2^{2k }\min_{\lambda\in\Delta(\mc{Z})}\max_{ w \in S_k }\|w^{\ast}-w\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}} \nn\\
\ge&\fr{1}{64\ur{\log\rbr{4\Delta^{-1}}}}\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}2^{2k }\min_{\lambda\in\Delta(\mc{Z})}\max_{ w , w '\in S_k }\|w-w'\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}} \nn\\
=&\fr{1}{64\ur{\log\rbr{4\Delta^{-1}}}}\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}2^{2k }\rho(S_k ), \label{ieq:rho*}
\end{align}

where the last inequality is by the triangle inequality, i.e.,

\begin{align*}
\max_{w,w'\in S_k}\vcn{w-w'}^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}\le&\max_{w,w'\in S_k}\rbr{\vcn{w-w^*}_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}+\vcn{w^*-w'}_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}}^2\\
\le&4\max_{w\in S_k}\vcn{w-w^*}^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}.
\end{align*}

Combining (<ref>) and (<ref>), we have

\begin{align*}
\sum_{k=1}^{\ur{\log\rbr{4\Delta^{-1}}}}N_k\le c(1+\omega)L_{\nu}\log\rbr{4\Delta^{-1}}\log\rbr{\fr{\log\rbr{4\Delta^{-1}}^2\abr{\mc{W}}}{\delta}}\rho^*+\log\rbr{4\Delta^{-1}}r(\omega).
\end{align*}

§ PROOFS OF SAMPLE COMPLEXITY WHEN GIVEN

Confounded pure exploration with known $\widehat{\Gamma}$

Input $\mc{Z}, \mc{W}, \widehat{\Gamma}, \delta, \epsilon, L_{\nu}\ge\sigma_{\nu}^2, L_{\eta}\ge\sigma_{\eta}^2$
Initialize: $k=1, \hat{\mc{W}}_k=\mc{W}$, calculate $\gamma$ using (<ref>)
$f(w,w',\Gamma,\lambda):=\|w-w'\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}$
$\exists w, w'\in\mc{W}_k,\text{s.t.},\ibr{w'-w}{\widehat{\theta}_{\texttt{oracle}}^k}>4\gamma$ or $\zeta_k>\gamma$
$N_k=\ur{2(1+\omega)\rbr{\zeta_k^{-2}\wedge\fr{1}{\gamma^2}}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}}\vee r(\omega)$
Pull arms in $Z_{N_k}$ and observe $Y_{N_k}$
Compute $\widehat{\theta}_{\texttt{P-2SLS}}^k=\rbr{Z_{N_k}^\top Z_{N_k}\widehat{\Gamma}}^{-1}Z_{N_k}^\top Y_{N_k}$
$\mc{W}_{k+1}=\mc{W}_k\backslash\cbr{w\in\mc{W}_k\mid \exists w'\in\mc{W}_k,\text{s.t.},\ibr{w'-w}{\widehat{\theta}_{\texttt{P-2SLS}}^k}>\zeta_k+\gamma}$
Output: any $w\in\mc{W}_k$

Suppose that we have $\widehat{\Gamma}$ that is an OLS estimate from an offline dataset $\cbr{Z_T, X_T}$ collected non-adaptively through a fixed design $\xi$ and the efficient rounding procedure $\text{ROUND}$,as well as $T\ge\fr{4\sigma_\eta^2}{\sigma_{\min}\rbr{A(\xi, \Gamma)}}\overline{\log}\rbr{Z_{T},\delta/2}$. Algorithm <ref> guarantees that with probability at least $1-\delta$, a $6\gamma$-good arm is returned, where

\begin{align}
\gamma:= \max_{w,w'\in\mc{W}}\vcn{w-w'}_{\bar{A}(Z_T, \widehat{\Gamma})^{-1}}\vcn{\theta}_2\sqrt{L_\eta\overline{\log}\rbr{Z_T,\delta}}. \label{eq:gamma}
\end{align}

Also, the algorithm terminates in at most

\begin{align*}
c(1+\omega)L_{\nu}\log(1/\delta)\rho^{\ast}(\gamma)+c(1+\omega)^2L_{\nu}\log(1/\delta)\fr{\sigma_\eta^2\overline{\log}\rbr{Z_{T},\delta}}{\sigma_{\min}\rbr{A(\xi, \Gamma)}}\fr{\rho(\xi,\gamma)}{T}\vee cr(\omega)
\end{align*}

samples, where $c$ hides logarithmic factors of $\Delta$ and $\abr{\mc{W}}$ and $\gamma$.

\begin{align*}
\rho(\lambda,\gamma):=\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}}\frac{\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2\vee\gamma^2}
\end{align*}

\begin{align*}
\rho^{\ast}(\gamma)=\min_{\lambda\in\Delta(\mc{Z})}\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}}\frac{\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2\vee\gamma^2}.
\end{align*}

The proof can be divided into four steps:

* The best arm $w^*$ is never eliminated.

* At the end of stage $k$, any suboptimal arm $w$ that satisfies $\ibr{\theta}{w^*-w}\ge2\zeta_k+2\gamma$ is eliminated.

* The stopping condition is met in finite time.

* When the stopping condition is met, there are only $6\gamma$-good arms left.

* The upper bound of the sample complexity.

Step 1: The best arm $w^*$ is never eliminated.

By the confidence interval in <ref>, we have, with probability at least $1-\fr{\delta}{2 k^2\abr{\mc{W}}}$, for any $k$ and $w\in\mc{W}_k$,

\begin{align*}
\abr{\ibr{\widehat{\theta}_{\texttt{P-2SLS}}^k-\theta}{w-w^*}}&\le\vcn{w-w^*}_{\bar{A}(Z_{N_k}, \widehat{\Gamma})^{-1}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4 k^2\abr{\mc{W}}}{\delta}}}+\gamma\\
&\le\fr{\sqrt{1+\omega}\vcn{w-w^*}_{A(\lambda_k, \widehat{\Gamma})^{-1}}}{\sqrt{N_k }}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4 k^2\abr{\mc{W}}}{\delta}}}+\gamma\\
&\le\fr{\sqrt{1+\omega}\vcn{w-w^*}_{A(\lambda_k, \widehat{\Gamma})^{-1}}}{\sqrt{\ur{2(1+\omega)\zeta_k^{-2}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}}\vee r(\omega)}}\sqrt{2\sigma_{\nu}^2\log\rbr{\fr{4 k^2\abr{\mc{W}}}{\delta}}}+\gamma\\
\end{align*}

Define a good event $\mc{E}_{k, w }$ for each $k$ and $ w \in\mc{W}_k$ as

\begin{align*}
\mc{E}_{k, w }=\cbr{\abr{\ibr{\widehat{\theta}_{\texttt{P-2SLS}}^k-\theta}{w-w^*}}\le\zeta_k+\gamma }.
\end{align*}

By the same calculation as (<ref>), we claim that with probability at least $1-\delta$, the event $\bigcap_{k=1}^{\infty}\bigcap_{w\in\mc{W}_k}\mc{E}_{k,w}$ holds. Under the event $\bigcap_{k=1}^{\infty}\bigcap_{w\in\mc{W}_k}\mc{E}_{w,k}$, to show that the best arm is never eliminated, it suffices to show that for any sub-optimal arm $w\in\mc{W}_k$,

\begin{align*}
\ibr{w-w^*}{\widehat{\theta}_{\texttt{P-2SLS}}^k}&=\ibr{w-w^*}{\widehat{\theta}_{\texttt{P-2SLS}}^k-\theta}+\ibr{w-w^*}{\theta}\\
\end{align*}

Thus the best arm $w^*$ never satisfies the elimination condition.

Step 2: At the end of stage $k$, any suboptimal arm $w$ that satisfies $\ibr{\theta}{w^*-w}\ge2\zeta_k+2\gamma$ is eliminated.

To prove this, we show that such arm $w$ must satisfy the elimination condition,

\begin{align*}
\max_{w'\in \mc{W}_k}\ibr{\widehat{\theta}_{\texttt{P-2SLS}}^k}{w'-w}\ge&\ibr{\widehat{\theta}_{\texttt{P-2SLS}}^k}{w^*-w}\\
\ge&-\zeta_k-\gamma+2\zeta_k+2\gamma\\
\end{align*}

Thus the arm $w$ is eliminated.

Step 3: The stopping condition is met in finite time.

Given the result in Step 2 and the fact that $\zeta_k$ is an exponentially decreasing sequence, we know that all of the arms $w$ satisfying $\ibr{\theta}{w^*-w}>2\gamma$ will be eliminated in finite time. This means that only the arms $w$ satisfying $\ibr{\theta}{w^*-w}\le2\gamma$ will remain. We need to show that $\forall w, w'\in\mc{W}_k,\ibr{w'-w}{\widehat{\theta}_{\texttt{P-2SLS}}^k}\le4\gamma$ can be achieved in finite time. When $\zeta_k\le\gamma$, (which will happen in finite time),

\begin{align*}
\ibr{w'-w}{\widehat{\theta}_{\texttt{P-2SLS}}^k}=&\ibr{w'-w+w^*-w^*}{\widehat{\theta}_{\texttt{P-2SLS}}^k}\\
\le&\zeta_k+\gamma+2\gamma+\zeta_k+\gamma\\
\end{align*}

Step 4: When the stopping condition is met, there are only $6\gamma$-good arms left.

For any $w\in\mc{W}_k$, we have

\begin{align*}
\ibr{w^*-w}{\theta}=&\ibr{w^*-w}{\theta-\widehat{\theta}_{\texttt{P-2SLS}}^k}+\ibr{w^*-w}{\widehat{\theta}_{\texttt{P-2SLS}}^k}\\
\le&\zeta_k+\gamma+4\gamma\\
\end{align*}

Step 5: The upper bound of the sample complexity.

Define $\mc{W}\rbr{2\gamma}$ as the set of $2\gamma$-good arms, i.e., $\mc{W}\rbr{2\gamma}:=\cbr{w\in\mc{W}\mid \ibr{\theta}{w^*-w}\le2\gamma}$. Then the best arm in the set $\mc{W}\backslash\mc{W}\rbr{2\gamma}$ has a suboptimality gap $\min_{w\in\mc{W}\backslash\mc{W}\rbr{2\gamma}}\ibr{\theta}{w^*-w}$. We define $\Delta_{\min}\rbr{2\gamma}:=\min_{w\in\mc{W}\backslash\mc{W}\rbr{2\gamma}}\ibr{\theta}{w^*-w}-2\gamma$. By the result in Step 2, we know that after $k^*:=\ur{\log\rbr{4\Delta_{\min}\rbr{2\gamma}^{-1}}}$ stages, all of the arms in $\mc{W}\backslash\mc{W}\rbr{2\gamma}$ are eliminated. Define $S_k =\cbr{ w \in\mc{W}\mid \ibr{ w ^*- w }{\theta}\le4\zeta_k+2\gamma}$. Thus with probability at least $1-\delta$, we have $\bigcap_{k}\cbr{\mc{W}_k \subseteq S_k }$.

The sample complexity of the algorithm is the total number of samples pulled, which is

\begin{align}
&\sum_{k=1}^{k^*}N_k \nn\\
=&\sum_{k=1}^{k^*}\rbr{\ur{2(1+\omega)\rbr{\zeta_k^{-2}\wedge\fr{1}{\gamma^2}}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}}\vee r(\omega)} \nn\\
\le&\sum_{k=1}^{k^*}\rbr{3(1+\omega)\rbr{2^{2k }\wedge\fr{1}{\gamma^2}}\rho(\mc{W}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}\vee r(\omega)} \nn\\
\le&\sum_{k=1}^{k^*}\rbr{3(1+\omega)\rbr{2^{2k }\wedge\fr{1}{\gamma^2}}\rho(\mc{S}_k)L_{\nu}\log\rbr{\fr{4k^2\abr{\mc{W}}}{\delta}}\vee r(\omega)}. \label{ieq:sum of rho}
\end{align}

Note that the factor of $\min_{\lambda\in\Delta(\mc{Z})}\max_{ w , w '\in S_k }\|w-w'\|^2_{A(\lambda, \Gamma)^{-1}}$ has the underlying true $\Gamma$ in it, while the algorithm uses the plugged-in $\widehat{\Gamma}$. We need to relate it to $\rho(S_k):=\min_{\lambda\in\Delta(\mc{Z})}\max_{w,w'\in S_k}\|w-w'\|^2_{A(\lambda, \widehat{\Gamma})^{-1}}$. By <ref>, and defining $\lambda_z^*(S_k)$ as the optimal design for $S_k$, i.e.,

\begin{align*}
\lambda_k^*=\argmin_{\lambda\in\Delta(\mc{Z})}\max_{w,w'\in S_k}\|w-w'\|^2_{A(\lambda, \Gamma)^{-1}}.
\end{align*}

Define $\lambda'_k=\fr{1}{2}\lambda_k^*+\fr{1}{2}\xi$, we have

\begin{align}
&\max_{w,w'\in W_k}\vcn{w-w'}^2_{A(\lambda_k, \widehat{\Gamma})^{-1}} \nn\\
=&\min_{\lambda\in\Delta(\mc{Z})}\max_{w,w'\in W_k}\vcn{w-w'}^2_{A(\lambda, \widehat{\Gamma})^{-1}} \nn\\
\le&\min_{\lambda\in\Delta(\mc{Z})}\max_{w,w'\in S_k}\vcn{w-w'}^2_{A(\lambda, \widehat{\Gamma})^{-1}} \nn\\
\le&\min_{\lambda\in\Delta(\mc{Z})}\max_{w,w'\in S_k}3\vcn{w-w'}^2_{A(\lambda, \Gamma)^{-1}}+2\max_{w,w'\in S_k}\vcn{(\Gamma^{-\top} - \widehat{\Gamma}^{-\top})\rbr{w-w'}}^2_{A(\lambda, I)^{-1}} \nn\\
\le&\max_{w,w'\in S_k}3\vcn{w-w'}^2_{A(\lambda'_k, \Gamma)^{-1}}+2\max_{w,w'\in S_k}\vcn{(\Gamma^{-\top} - \widehat{\Gamma}^{-\top})\rbr{w-w'}}^2_{A(\lambda'_k, I)^{-1}} \nn\\
\le&\max_{w,w'\in S_k}6\vcn{w-w'}^2_{A(\lambda_k^*, \Gamma)^{-1}}+\fr{16\sigma_\eta^2\overline{\log}\rbr{Z_{T},\delta/2}}{\sigma_{\min}\rbr{A(\xi, \Gamma)}}\max_{w,w'\in S_k}\vcn{w-w'}^2_{A(\xi, \Gamma)^{-1}}\fr{1}{T}, \label{ieq:lmabdaK2lamdbda*}
\end{align}

where the last inequality is due to

\begin{align*}
\vcn{w-w'}^2_{A(\lambda'_k, \Gamma)^{-1}}\le&\vcn{w-w'}^2_{A\rbr{\fr{1}{2}\lambda_k^*+\fr{1}{2}\xi, \Gamma}^{-1}}\\
\le&\vcn{w-w'}^2_{A\rbr{\fr{1}{2}\lambda_k^*, \Gamma}^{-1}}\\
=&2\vcn{w-w'}^2_{A(\lambda_k^*, \Gamma)^{-1}},
\end{align*}

as well as, by <ref>, when $T\ge\fr{4\sigma_\eta^2}{\sigma_{\min}\rbr{A(\xi, \Gamma)}}\overline{\log}\rbr{Z_{T},\delta/2}$, we have, with probability at least $1-\delta/2$,
we have

\begin{align*}
&\max_{w,w'\in S_k}\vcn{(\Gamma^{-\top} - \widehat{\Gamma}^{-\top})\rbr{w-w'}}^2_{A(\lambda'_k, I)^{-1}}\\
\le&\fr{4\sigma_\eta^2\overline{\log}\rbr{Z_{T},\delta/2}}{\sigma_{\min}\rbr{A(\lambda'_k, \Gamma)}}\max_{w,w'\in S_k}\vcn{w-w'}^2_{A(\xi, \Gamma)^{-1}}\fr{1}{T}\\
\le&\fr{4\sigma_\eta^2\overline{\log}\rbr{Z_{T},\delta/2}}{\sigma_{\min}\rbr{A(\alpha\lambda_k^*+(1-\alpha)\xi, \Gamma)}}\max_{w,w'\in S_k}\vcn{w-w'}^2_{A(\xi, \Gamma)^{-1}}\fr{1}{T}\\
\le&\fr{8\sigma_\eta^2\overline{\log}\rbr{Z_{T},\delta/2}}{\sigma_{\min}\rbr{A(\xi, \Gamma)}}\max_{w,w'\in S_k}\vcn{w-w'}^2_{A(\xi, \Gamma)^{-1}}\fr{1}{T}.
\end{align*}

We can lower bound $\rho^{\ast}(\gamma)$ by

\begin{align}
\rho^{\ast}(\gamma)=& \min_{\lambda \in \Delta(\mc{Z})}\max_{w\in \mathcal{W}\setminus \{w^{\ast}\}} \frac{\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2 \vee \gamma^2} \nn\\
=&\min_{\lambda\in\Delta(\mc{Z})}\max_{k}\max_{ w \in S_k\setminus \{w^{\ast}\} }\frac{\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2 \vee \gamma^2} \nn\\
\ge&\fr{1}{k^*}\min_{\lambda\in\Delta(\mc{Z})}\sum_{k=1}^{k^*}\max_{ w \in S_k\setminus \{w^{\ast}\} }\frac{\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}}}{\langle w^{\ast}-w, \theta\rangle ^2 \vee \gamma^2} \nn\\
\ge&\fr{1}{16k^*}\sum_{k=1}^{k^*}\rbr{2^{2k }\wedge\fr{1}{\gamma^2}}\min_{\lambda\in\Delta(\mc{Z})}\max_{ w \in S_k\setminus \{w^{\ast}\} }\|w^{\ast}-w\|^2_{A(\lambda, \Gamma)^{-1}} \nn\\
\ge&\fr{1}{64k^*}\sum_{k=1}^{k^*}\rbr{2^{2k }\wedge\fr{1}{\gamma^2}}\min_{\lambda\in\Delta(\mc{Z})}\max_{ w , w '\in S_k }\|w-w'\|^2_{A(\lambda, \Gamma)^{-1}}. \label{ieq:lowerbound rho 1}
\end{align}

Given (<ref>), (<ref>) and (<ref>), we have

\begin{align*}
&\sum_{k=1}^{k^*}N_k\le c(1+\omega)L_{\nu}\log(1/\delta)\rho^{\ast}(\gamma)+c(1+\omega)^2L_{\nu}\log(1/\delta)\fr{\sigma_\eta^2\overline{\log}\rbr{Z_{T},\delta}}{\sigma_{\min}\rbr{A(\xi, \Gamma)}}\fr{\rho(\xi,\gamma)}{T}+cr(\omega),
\end{align*}

where $c$ hides logarithmic factors of $\Delta$ and $\abr{\mc{W}}$ and $\gamma$.

§ PROOFS OF SAMPLE COMPLEXITY WITH UNKNOWN

Optimal design with unknown $\Gamma$

Input $\mc{Z}, \mc{W}, \delta, L_{\nu}\ge\sigma_{\nu}^2, L_\eta\ge\sigma_{\eta}^2, \omega, \gammin\le\lammin(\Gam), \lambda_E,\kappa_0$
Initialize: $k=1, \mc{W}_1=\mc{W},\widehat{\Gamma}_0=\perp, \zeta_1=1$
Define $f(w,w',\Gamma,\lambda):=\|w-w'\|^2_{(\sum_{z\in \mathcal{Z}}\Gamma^{\top}\lambda_z zz^{\top}\Gamma)^{-1}}$,
$M:=\fr{32L_\eta}{\gammin^2\sigma_{\min}\rbr{A(\lambda_E, I)}}\vee1$, $\dt_{k,\ell}:=\fr{\delta}{4k^2\ell^2}$
$\widehat{\Gamma}_k=\Gamma-\mathtt{estimator}\rbr{\mc{W}_k, \widehat{\Gamma}_{k-1}, \zeta_k, \delta, k, \omega,\lambda_E, M, L_\eta}$ Step 1: update $\widehat{\Gamma}$
$\widehat{\theta}_{\texttt{P-2SLS}}^k=\theta-\mathtt{estimator}\rbr{\mc{W}_k, \delta, \zeta_k, \widehat{\Gamma}_{k}, \omega,L_{\nu},k}$ Step 2: update $\widehat{\theta}$
$\mc{W}_{k+1}=\mc{W}_k\backslash\cbr{w\in\mc{W}_k\mid \exists w'\in\mc{W}_k,\text{s.t.},\ibr{w'-w}{\widehat{\theta}_{\texttt{P-2SLS}}^k}>\zeta_k}$ Step 3: elimination
$k\leftarrow k+1$, $\zeta_k=2^{-k}$
Output: $\mc{W}_k$

The algorithms of $\Gamma-\mathtt{estimator}$ and $\theta-\mathtt{estimator}$ we present below are slightly different from the one in the main text. In the main text, we omit the phase index $k$ in the algorithm for simplicity.

Input $\mc{W}_k, \widehat{\Gamma}_{k-1}, \zeta_k, \delta, k, \omega,\lambda_E, M, L_\eta$
Define $\mathtt{Stop}\rbr{\mc{W},Z,\Gamma,\dt}:=\max_{w,w'\in\mc{W}}\vcn{w-w'}_{\bar{A}(Z, \Gamma)^{-1}}\vcn{\theta}_2\sqrt{L_{\eta}\overline{\log}\rbr{Z,\delta}}$
Initialize $\ell=1$, $N_{k,0,0}=0$ doubling trick initialization
$\ell=1$ or $\mathtt{Stop}\rbr{\mc{W}_k,Z_{k,0,\ell},\widehat{\Gamma}_k,\dt_{k,\ell}}>1$
get $2^{\ell-1}\rbr{r(\omega)\vee\fr{2}{\kappa_0}}$ samples denoted as $\cbr{Z_{k,0,\ell},X_{k,0,\ell},Y_{k,0,\ell}}$ per design $\lambda_E$ via ROUND
Update $\widehat{\Gamma}_k$ by  on $\cbr{Z_{k,0,\ell},X_{k,0,\ell}}$, $\ell\leftarrow \ell+1$
$N'=\drd{4gdM\ln\del{1+2M\rbr{d+L_z^2}+2M2gdM}+8M\ln\del{\fr{2\cd6^d}{\dt_{k,1}}}\vee r(\omega)}$
$\ell=1$ or $\mathtt{Stop}\rbr{\mc{W}_k,Z_{k,0,\ell}\cup Z_{k,1,\ell},\widehat{\Gamma}_k,\dt_{k,\ell}}>\zeta_k$
$N_{k,1,\ell}=2^{\ell}N'$ doubling trick update
get $N_{k,1,\ell}$ samples per $\tilde{\lambda}_k$ denoted as $\cbr{Z_{k,1,\ell},X_{k,1,\ell},Y_{k,1,\ell}}$ via ROUND
$N_{k,0,\ell}=\urd{2gdM\ln\del{M\rbr{d+N_{k,1,\ell}+L_z^2}}+4M\ln\del{\fr{2\cd6^d}{\dt_{k,\ell}}}\vee r(\omega)\vee\fr{2}{\kappa_0}}$
get $(N_{k,0,\ell}-N_{k,0,\ell-1})$ samples per $\lambda_E$ augmented to $\cbr{Z_{k,0,\ell-1},X_{k,0,\ell-1}}$ and get $\cbr{Z_{k,0,\ell},X_{k,0,\ell}}$