Finite Population Regression Adjustment and Non-asymptotic Guarantees for Treatment Effect Estimation
Authors: Mehrdad Ghadiri
David Arbour
Tung Mai
Cameron N Musco
Anup Rao
Venue: NeurIPS 2023 poster (2023)
Keywords: regression adjustment; treatment effect estimation; average treatment effect
TL;DR: We present regression adjustment methods based on numerical linear algebra techniques for experimental design, average treatment effect, and individual treatment effect estimation.

Abstract: The design and analysis of randomized experiments is fundamental to many areas, from the physical and social sciences to industrial settings. 
Regression adjustment is a popular technique to reduce the variance of estimates obtained from experiments, by utilizing information contained in auxiliary covariates. 
While there is a large literature within the statistics community studying various approaches to regression adjustment and their asymptotic properties, little focus has been given to approaches in the finite population setting with non-asymptotic accuracy bounds. 
Further, prior work typically assumes that an entire population is exposed to an experiment, whereas practitioners often seek to minimize the number of subjects exposed to an experiment, for ethical and pragmatic reasons.
In this work, we study the problems of estimating the sample mean, individual treatment effects, and average treatment effect with regression adjustment. 
We propose approaches that use techniques from randomized numerical linear algebra to sample a subset of the population on which to perform an experiment. We give non-asymptotic accuracy bounds for our methods and demonstrate that they compare favorably with prior approaches.

Finite Population Regression Adjustment and Non-asymptotic Guarantees for Treatment Effect Estimation
Mehrdad Ghadiri MIT mehrdadg@mit.cm edu
David Arbour Adobe Research arbour@adobe.cm com
Tung Mai Adobe Research tumai @adobe. com
Cameron Musco UMass Amherst cmuscoθcs umass. edu
Anup B. Rao
Adobe Research
anuprao@adobe. com
Abstract
The design and analysis of randomized experiments is fundamental to many from the physical and social sciences to industrial settings. Regression adjustment is popular technique to reduce the variance of estimates obtained from experiments, by utilizing information contained in auxiliary covariates. While there is a large literature within the statisties community studying various approaches to regression adjustment and their asymptotic properties, little focus has been given to approaches in the finite population setting with non-asymptotic accuracy bounds. Further, prior work typically assumes that an entire population exposed to an experiment, whereas practitioners often seek to minimize the number of subjects exposed to experiment, for ethical and pragmatic reasons. In this work, we study the problems of estimating the sample mean, individual treatment effects, and average treatment effect with regression adjustment. We propose approaches that use techniques from randomized numerical linear algebra to sample subset of the population on which to perform an experiment. We give non-asympiotic accuracy bounds for our methods and demonstrate that they compare favorably with prior approaches.
Introduction
Randomized experiments play crucial role in estimating counterfactual outcomes. The strength of randomization lies in its ability to ensure independence between the treatment group and pre- treatment covariates [14] (21). This independence, in turn, results in unbiased estimation of treatment effects through the comparison sample averages. However, in finite sample settings, the use of simple randomization techniques result in high variance of the estimated quantities. There have been two broad approaches to address this issue by reducing the variance the estimated effects when subject level covariates
available: (1) design- based approaches, which focus the mechanism for assigning individuals to treatment and control groups in order to minimize imbalance (c.f. 21 18 27 34 5 (5)), and (2) regression adjustment based approaches, which correct for imbalances post hoc by incorporating a regression of the covariates on the outcome (c.f.
Much of the prior work examining variance reduction for experimentation has focused asymptotic properties of approaches [23]24] 13 23] Recently, there has been work focusing on basign-basce variance reduction, which ties the treatment assignment problem to discrepancy minimization and gives finite sample error bounds for estimating causal estimands such as the average treatment effect (ATE) 18 and individual treatment effect (TTE) In particular, Harshaw et al. propose
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

design based on the Gram-Schmide Schmidd random walk (GSW) and give finite sample guarantees on the mean squared error for average treatment effect estimation when using this design with the classic Horvizz-Thompsnn estimator. Addanki, et al. [ use GSW and leverage score sampling to provide algorithms and finite sample guarantees for both average and individual treatment effect estimation in the partial observation setting where only a subset of subjects are included in the experiment.
However, little work has examined the finite sample properties for experimental designs based on regression adjustment. Mou, et al. [28] provide non-asymptotic guarantees, but the analysis is in infinite population settings. The main contribution of this work to analyze the finite population behavior of regression adjustment for variance reduction in treatment effect estimation, for ATE and ITE estimation in both the full and partial (sub-pope soab-oooplatio) settings.
1.1 Problem Statement: Finite Population Treatment Effect Estimation
In the finite population treatment effect estimation problem, we have n units (or individuals) associated with covariates X1...., X,, E Rd which are the rows of covariate matrix X E Raxd Each unit i is associated with two real numbers yi Vi control, respectively. We denote the corresponding potential outcome vectors with y()), y E IR.. The Individual Treatment Effeett (ITE) vector E IR" is defined as t := t=:y"" y( (() and the Average Treatment Effect (ATE) is defined as the mean of this vector, i.e.,
(0) called the potential outcomes under the treatment and
2i=1
In treatment effect estimation, we seek to estimate either the ITE or ATE, under the constraint that at most one of the potential outcomes, either y, or y, can be observed for each unit i. E.g., that unit may be assigned treatment group in controlled trial, and thus y is observed while y! is not.
In full observation treatment effect estimation, we observe one potential outcome for each unit 24 [8]. In partial observation treatment effect estimation, we further restrict ourselves to observing outcomes for just a small subsample of the population. This setting is important when the goal is to minimize experimental costs by e.g., limiting the size of controlled
in the population 9 13
trial (3 20 32
To reduce variance, many approaches to ATE and ITE estimation, including our own, attempt to leverage the population covariates X make inferences about the treatment outcomes. In particular, throughout this work, we will target error bounds that depend on the best linear fit of the treatment outcomes to our covariates. Such linear effects models are common in the literature on treatment and provide useful baseline for establishing theoretical bounds.
effect estimation
33 18
1.2 Our Contributions
We give new approaches to both ITE and ATE estimation in the full and partial observation settings that combine the classic statistical technique of regression adjustment with techniques from randomized numerical linear algebra. Our algorithms give natural, non-asymptotic error bounds depending on the best linear fit for either the ITE vector t == y(1) (0) or the sum of the potential outcome vectors µ := (l) + (0) our covariates. µu is a natural quantity that arises e.g., in variance bounds for classic approaches to ATE estimation, such as the Hovitz-Thomppe estimator [19]
1.2.1 Individual Treatment Effect Estimation
Our first result for individual treatment effect estimation in the full and partial observation settings. Theorem 1 (ITE Estimation). For any e, ô, a E (0, 1], there exists randomized algorithm (Algorithm that observes a potential outcome for (ddoo(d//))/ + a n) units and outputs vector t. Moreover there is event E such that Pr(E) > - (over the randomness of the algorithm) and
E [ [1B]E] < (1 -minbeRd XXb -112 _. 1.1..11100 (1+ --(0)
Theorem shows that with high probability (at least - 5). our algorithm achieves an ex- pected error bound that nearly matches the best linear fit of the ITE vector to our covariates (i.e.. minbeR" || Xbb- t) up to an additive term depending on the maximum magnitude of the sum of potential outcomes for any individual. In the full observation setting, when 1, this additive

---

term is O(d), assuming that the maximum magnitude does not scale with the population size n. In contrast, expect the best fit error (minbeR= | - term is lower order. In the partial observation setting, as long as we set a = w(1/n), we still expect this additive term to be lower order. Thus, the case when d < n, we able nearly match the best fit error while only observing a very small subset of the full population.
-t t[l2] to grow (()), and thus the additive
We sketch the proof of Theorem in Section 2 Our algorithm is based pseudo-outccoe regression, which is studied in the asymptotic setting by Kennedy et al. for conditional ATE estimation [22]. Roughly, in the full observation setting, we construct a vector where U; = 2y;
(1) with probability
1/2 and U/ = 22 (0) with probability 1/2. We can see that E[v] = t (i.e., v is equal to the ITE vector in expectation), and importantly, that constructing v only requires observing one potential outcome for each individual. By regressing v onto our covariates, we obtain our estimate È for t.
In the partial observation setting, further subsample individuals according to the leverage scores of the covariate matrix X. This is a standard technique in randomized linear algebra [26] 10 35), which allows us to approximately solve the above regression problem while only observing a subset of Oddee²) entries of v (and thus only observing outcomes for a subset of the population). Addanki et al. [1] similarly propose algorithm for the partial observation setting based on leverage score sampling. However, instead of pseudo-outcome regression, they learn separate approximations to both outcome vectors y(l) and (0)), Thus, their final error depends on the best linear fit error for both these vectors, rather than the more natural best linear fit error for the ITE vector t.
1.2.2 Average Treatment Effect Estimation
We can extend our general approach to give bounds for average treatment effect estimation. In the full observation setting we show:
Theorem 2 (ATE Estimation - Full Observation). For any e, E (0, 1J, there exists a randomized algorithm (Algorithm D that computes an unbiased estimate of the ATE t, and for which there is an event E with Pr(E) 1 - 8 (over the randomness of the algorithm) such that
E [(f (( 1 8(1+c) mimbER" (IIXb
11112 1001og(n/5)) 5² 1b(12)
32d n²
(1)
(0))
where Ç == max ||x:|l and μ == y0)) + y() is the total outcome vector.
As with our result in ITE estimate, the error bound of Theorem 2 matches the error of the best linear fit of the total outcome vector μs to our covariates up to a constant factor plus additive terms depending on 1) the maximum magnitude of the sum of potential outcomes for any individual, and 2) the norm of the coefficient vector used to approximately reconstruct from the covariates. Again, we expect these additive terms to scale as O(1/n²)) while we expect the term depending on the best fit error to scale as (1/nn). Thus, we generally expect the additive error to be lower order.
Harshaw et al. [1] proposed algorithm based the Gram-Schmiae walk (GSW) design (see Algorithm 5 in the appendix) for balancing assignments to treatment and control groups that achieves an unbiased estimator for ATE with variance of J ming | Xb| 2 + f ||b|2] where á (0, 1) is chosen by the experiment designer. This guarantee is comparable to but stronger than ours, e.g., when one sets Φ = 1/2. If we ignore the additive error terms and focus just on the best linear fit error |Xbb 14122- then GSW is better than our guarantees by about a factor of 8. But the GSW design is computationally more expensive and, since is based on balancing the covariates, requires availability of all covariates before running the experiment and hence cannot be applied in an online setting 4 contrast, our approach can be applied directly the online setting since we place each unit independently in the treatment or control groups with equal probability. Moreover, as we show in our experiments, the empirical performance of our algorithm is much closer to GSW than the theoretical bounds we are able to show.
We sketch the ideas behind Theorem 2 Section 2 Our algorithm is based on the classic Horvitz- Thompson estimator 19]. Roughly, this estimator randomly assigns individuals to control and treatment groups. and estimates ATE as the difference in the average outcome between these groups, appropriately weighted by the assignment probabilities. It is well known that the variance of this estimator is bounded by + ||11/2 To reduce this variance, we introduce filloduceeilllooobsratttoooo regression
3

---

adjusted Horvitz-Thompson (RAHT) estimator, which estimates ATE as the following, We partition the units into two groups 5 and S (this is a different partitioning than the partitioning into control and treatment), and we regress Mı and Mg onto their corresponding covariates. As in the ITE case, since we cannot directly form μ4, we instead perform the above regressions using random vector u with [E[u] = μ., and where u can be formed by only observing one potential outcome per individual. We then use the solution vector of group S to adjust the outcomes of the group S and vice versa, and we apply the HT estimator to the adjusted outcomes. This gives the bound in Theorem 2
Again, following a similar approach as in the ITE case, we can apply leverage score based subsampling to tackle ATE estimation in the partial observation setting as well, obtaining:
Theorem 3 (ATE Estimation Partial Observation). For any 6, E (0, 1). Φ E (0, and m E [n], there exists a randomized algorithm (Algorithm 4 that observes potential outcome for (og//////))/ + m) units and outputs an unbiased estimate & of the ATE т, for which there is an event & with Pr(E) 2 1 - 8 (over the randomness of the algorithm), such that
E [(f [|] -
1 |Xb* - -1112 - m² (1-0) ||b^112 +
1 c²
100d-log(d/6) n²²²
||
mn
+ (1+c) m 11 + 1 II(X X}b - (t - t)113 + mn
A
[B1B)
where Ç := maxie(m] ||x:||2. 2: 1t == y()) + y(u), A G-logdd/0) <2 t E 12" is a vector where all the
entries are equal to I = 2/4 -i==1 Li and
b arg min BERd
Xbb m1l2 + c² [ 6 arg min [[[(X - [I(X-X-b- (t. E)IB All
n
bER"
The error bound of Theorem depends on 1) the error of best linear fits of the covariates onto μ and t - E, 2) the norm of the coefficient vector used to approximately reconstruct μa and - E from the covariates, and 3) the largest component of the total outcome vector μ.. We expect both the best linear fit and || terms to scale as O ).
). All other terms are lower order. Namely, The terms depending on the norms of the coefficient vectors scale as O( mm) and O =. For the case where md & n, by increasing the number of samples to O log (d/))/²* + md), we can make the || term to scale as ((1mm²) as well.
We deploy a different regression adjustment technique for achieving Theorem 3 which we call partial- observation regression adjusted Horvitz-Thompson (RAHT) estimator. In the partial observation setting, a simple approach is just to apply the full observation algorithm to uniformly selected subset of m units and then report the ATE estimate for these m units. However, this leads to a variance
bound depending on |t E E To reduce this variance, the partial observation RAHT adjusts the estimate obtained from the mi units using regression adjustment techniques. We give the details behind the full algorithm in Section 2.3 and provide the full analysis in the appendix.
1.2.3 Experimental Evaluation
We compliment the above theoretical bounds with experiments using both synthetic and real data, focusing on the full observation ATE and ITE estimation problems. For ATE, we compare our method with classical regression adjustment (24, the GSW method [18), and the algorithm of (28). We observe that our estimator (which unbiased) performs as well or better than all other estimators except the classic regression adjustment, which is a biased estimator. For ITE, we compare our approach with the method of [I] and observe superior performance. A discussion about the running time of our algorithm and comparison to GSW design is included in the appendix.
1.3 Roadmap
We present the high-level ideas and techniques for ITE and ATE estimation for the full and partial observation settings in Section [ We then present notation and preliminaries required for a more in-depth discussion of our result in Section 3 Section 4 we give a detailed non-asymptotic analysis of the random vector (i.e., pstudo- regression technique, which forms the basis for all of our algorithms. We present our algorithm for the full observation ATE estimation and sketch the proof of its accuracy (Theorem 2). Finally, we present our experimental results in Section 5

---

We relegate the details of other algorithms and proofs to the appendix, along with additional details and results for the experimental results. In the appendix. we also include warm-up result for the mean estimation problem in the partial observation setting (where there is only one outcome y and we seek estimate + behind leverage score sampling, which is used for our results on both ITE and ATE estimation. In addition we discuss how to remove the 1 - probability from our results in Section C.1
n >i= yr from a small number of samples). This result illustrates the key ideas
2 Technical Overview
In this section, we sketch the key ideas behind our main results on ITE estimation (Theorem I and ATE estimation (Theorems) 2 aadd] 3
2.1 ITE Estimation
We first consider ITE estimation in the full observation setting. The key challenge here is that we must infer the full treatment effect vector = y(()) - -y(0), while only observing one outcome, ( or (0) for any unit i. To do so, we use the idea of pxeudo-outcome regression in [22] We construct random vector v E II" that is equal to t in expectation, i.c., E [v] = but only requires observing one outcome per unit. Specifically, we independently set v; 2y (1) with probability 0.5 and = -2y; (0)
yi
with probability 0.5. We then regress v onto our covariates, and use the result to estimate t. In particular, we set b == arg minbERd Xb v 2 and construct our ITE estimate as f Xb.
To give non-asymptotic error bound for this approach, we observe that we can write v = t Ja, where E {1, 11" is random sign vector, fe = =y (y is the sum of potential outcome vectors, and @ denotes the Hadamard (entrywise) product. Thus, letting xx E Raxx be the projection matrix onto the column span of X, we can write b = T = xxt + TX (2 p)) In turn, we have
E
|it - [1 = E | [1[] = E ||tx|- t112 + Tx (z 2 2 2x x
X 2 11)
The first term the righthand side is simply the best linear fit enror for t, || tl2 minpeRd |Xb tl The third term is in expectation since E[z @ u1] = Finally, since the entries of z u are independent and mean 0, the expectation of the second term can be bounded by ir(Tx) |||||| d | Putting these all together, we have
[
[12] S minbER Xbb
t12 + || 11113
(I)
which gives the bound of Theorem in the full observation setting. We note that this expected error bound can also be turned into a high probability bound using concentration inequalities, i.e., Wright- inequality (31).
ITE Estimation with Partial Observation. Our next step is extend the above the partial ob- servation setting. To do so, we apply leverage score sampling, which is a standard technique for approximate regression via subsampling [10] 26 35). It is well known that if we sample ((dlog(d/))/e) rows of X and the corresponding entries of v according to the leverage scores of X, then we can solve reweighted regression problem just these rows to find 6 satisfying Xb - v 2
(1++)) minbER" Xb - v Observe that to compute b, we only need to observe potential outcomes for the units (i.e., the entries of v) that are sampled. Thus, this almost allows to recover similar bound to (I] up to a multiplicative (1 + e) factor, in the partial observation setting. However, the leverage score of any one row may be very small, and if sampled, that row may be scaled up by large factor in our reweighted regression problem (to preserve expectation). It thus becomes difficult to control to error introduced by approximating t by v.
To handle this issue, we mix leverage score sampling, with uniform sampling at rate a. This increases our sample complexity to (d log(d/))) + an), but allows to bound the error due to
approximating t by v by (I+e)d || This yields the final sample complexity and error bound of Theorem For a full proof, see Appendix
5

---

2.2 ATE with Full Observation
We build on our techniques for ITE estimation tackle the ATE estimation problem. We start with the classic Trriiiso (HT) estimator for ATE, defined below.
Definition 4 (Horvitz-hommsonn estimator). Given two outcome vectors y y yy0) E R", let Z
and Z [n] \ Z + be random partitioning of the units to two groups under distribution P. Then
y (0)
the Morvitz-Thompsse
(HT) estimator is 7 := 1,2 ((vEz+
FplEz where
n
P
Pp j E - denotes the probability of placing i in Z when the partitioning is performed according to distribution P.
Note that if units are assigned independently with equal probability to treatment and control, the HT estimator is equivalent to outputting the average of the random vector used in our regression-base estimate for ITE. In this case, the HT estimator is unbiased and well known to have variance bounded by n.a ll 14 where μ := y()) y()). Our goal is to improve this variance to instead depend on the error of the best linear fit to μa, minDER= |Xb classic technique of regression adjustment.
1412. as in Theorem To do so, we will employ the
Horvitz-Thompson Thompson on Regression Adjusted Outcomes. Suppose we were given vector b R determined independently of the randomness in the HT estimator. Then applying the HT estimator
directly on the regression adjusted vectors g((+) := (())) - Xb (which have the same ATE as the original vectors), would yield an unbiased estimate for the ATE with variance at most n² | - Naturally, we would like to find such b making this variance bound as small as possible.
To do so, we will apply the same approach as in ITE estimation. We cannot directly regress μ onto our covariates to find an optimal b, as we cannot observe any entry of ja since is the sum of the two potential outcomes. Instead, we can construct vector u with B(u] μe, and use it as surrogate in
our regression problem. In particular, we let u; = 2y/ (1) with probability 1/2 and u, = 2y, with probability 1/2. As in the ITE setting, we can bound the error incurred by approximating l by u in
(0)
our regression problem, and thus bound the variance of our regression adjusted HT estimator. We
note that to the best of our knowledge, this idea has not been leveraged in ATE estimation.
Unfortunately, a difficulty arises in this approach. We need to assign each unit a treatment or control group both when computing u to solve our regression problem, and later when applying the HT estimator. If we use the same assignment for both steps (as is required since we can only observe one of y; (1) or Ui
(0) for each unit), we introduce bias in our estimator.
Avoiding Estimator Bias. To alleviate this issue, use ridge leverage score sampling to first partition the units into two groups. Ridge leverage score sampling is a technique similar to lever- age score sampling that gives similar approximation guarantees to linear regression but for ridge regression problems (i.e., linear regression problems with an Ee regularization term). We pick a ridge (regularization) parameter that guarantees that each unit only needs to be in the first second group with probability of at most 0.5 while guaranteeing that the solution obtained from this sampling is (1 + e) approximation of the optimal ridge regression solution. Then we solve ridge regression problems on the groups separately to obtain two vectors ((2) and ((2), for the first and the second group, respectively. We then use (11) to adjust the outcomes of the second group and use ((2) to adjust the outcomes of the first group. Since the vector of each group is not used adjust the outcomes of itself, this gives an unbiased estimator.
2.3
ATE with Partial Observation
We now consider ATE estimation where we desire to only observe y, or (0) for small subset of the population. If draw a subset of size m uniformly at random and observe the outcomes of the samples in the subset independently and with equal probability and apply the Horvitz-Thompson
(1)
estimator, the variance is E [f - ?)? <
[
+ |t We would like to replace terms || and
||t - with min, |/ Xbll2 and min, ||t i x)bbi2 respectively. To deal with the ||aa||2 term, use the HT estimator applied to the GSW design [18] the uniformly selected subset of size m of the data. Then for & E (0, that the experiment designer selects, The variance reduces
to m minbberd
Xb-
[al] 1 + I- ||b|
+
||t-tl1ā To replace the t-Ellā term with the error mn
mm

---

of the best linear fit, we use the partial observation RAHT estimator that given solution vector b adjusts the estimate obtained by the GSW design on the sample of size m. Similar to our approach for partial observation ITE estimation, we compute by sampling about d units according to leverage scores of the matrix and performing the random vector regression gpseudo-outcom regression) on a vector that is equal to t/2 in expectation. Let S be the set of units selected by leverage score sampling and 5 = [n]] \ S. We denote the estimate obtained by sampling m units from (we denote this set with 7) and applying GSW design with Fg. Then, the RAHT estimator estimates the ATE as Pg. Eer - }).. of being included in S and, therefore, would contribute to the estimate with lower probability. To resolve this issue, we also estimate the ATE on the set S with the same random vector of outcomes that is used in regression to learn b. We use HT estimator for this estimate, and we denote it by îs. Then, our final estimate is convex combination of the two estimates as the following.
This estimate is biased because some units have a higher probability
f5s m (x x3 bb)
Tg.
This estimator is unbiased and achieves the variance bounds stated in Theorem 3
3 Preliminaries
Notation. Vectors and matrices are denoted with bold lowercase and capital letters, respectively. the Hadamard (entrywise) product of x, is denoted by x @ y. We denote the vectors of all ones and zeros in I by I, and O,, respectively. We drop the subscript the dimension is clear. We denote the identity matrix of dimension n by I, For X E RRadd we denote the projection matrix of X by π == X(X X), where + denotes the pseudoinverse. For x E IR" and S n, 2, denotes its 1th entry and X1: denotes a vector in RC that is equal to its first entries. We use X and x; interchangeably to denote row of X. Similarly for S C [m]. Xs: E IR ()xx denotes the submatrix with row indices in S and for y E R", Ys E Rlil denotes the vector restricted to indices in S.
Definition 5 (Ridge leverage scores (2).. Let E Raxd, X 2 O. and X; E IR be the ith row of Then the X-ridge leverage score of row i of X is defined as e!²) := x; (X X TXTx
AI) x;. We denote
the leverage scores (i.e., X = 0) with &, If corresponding matrix X is not clear from the context, denote the X-ridge leverage scores with e!A1 (X).
Note that & is the i'th diagonal entry of the projection matrix - = X(XTX)--XX, Thuss, since tr(TT) rank(XX) Σi= 4 rXn. Since in our approach, we need to sample different sets of population for different purposes, we require a bound on the ridge leverage scores that can be achieved by taking the regularization factor À to be large enough.
Theorem 6. Let Ç == max,e(n] ||x|||2- Thenn for X c. S2 e 2 1, G. (X) ] for all i E [m]]
The following theorem provides a sampling procedure based on X-ridge leverage scores to solve linear ridge regression problem approximately.
Theorem 7 12). Let 5 E 1), and u E R be a vector that upper bounds the X-ridge leverage scores of a matrix X E Raxdd i.e., e(A) S U₁ Let Raxn be a random diagonal matrix in which SA = 1/ with probability p. min u; log (d/8)fc}, and S = 0,, otherwise. Then for
any y E IR", letting b arg minbeR* ||SXb- Syl A [[b][] with probably at least 1 8.
xxb yll2 + {|lbi| & (1 + <) min
| - y12 + A [b]]
4 Average Treatment Effect Estimation
In this section, we first analyze the random vector regression (i.e., pseudo-outeome regression) approach which is used both for our ITE and ATE estimation. We then discuss our full observation ATE estimation using this approach and provide a sketch of the proof of Theorem 2
4.1 Random Vector Regression
For both ATE and ITE estimation, regression adjustment allows us to reduce the variance of estimation to the error of the best linear fit on vectors μ and t, respectively. Since we do not have access to the
7

---

entries of these vectors, we cannot compute solution vector directly. Therefore, for ATE (or ITE), our approach is to instead compute solution vector b by performing regression on random vector that, in expectation, is equal to ja (or t). The following theorem characterizes the error of the solution vector obtained in this way.
The following is more general result that combines this random assignment technique with leverage score sampling allow us to observe only small subset of the population. Algorithm 3 and the proof of the following theorem are presented in the appendix.
Theorem (Random vector regression). Let yr))) y(()) E R and y E ER" be a random vector such that for each i E [n]. y: is independently and with equal probability is either equal to (0) (1)
y;
Moreover, let b* arg min, ||Xb y11 Let μ = -y(1) + y() Then
E 1|2b - _113] 2 d|ly - -y + minb [ - -
Proof. Let 2 E 1-1, +13", where 24 = +1 if y/1 = yi (1) anddd 2; = -1, otherwise. First, note that
2y = + z @ t, where μ = ( + +y(0) and t = y(")) -y Therefore
μ 2Xb*
x = μ x + o t) = (μ TXH) - xx)
Moreover since I 7x and x are orthogonal to each other,
|1 2Xb b0bl1 1/2 (1 xx ||
1 x min, X - 11 So we only need to bound E [ の 11 Now note that 1144-1100 -
We have (x @ t)l? (z ② t) (2 t) z T x Tz, where T is the diagonal matrix associated with the vector t. Since 2; and 2j are independent for i + j. we have
2 T Tzj = Ei= (Txx) sa S ||t||20 (Tx)
Then the result follows by noting that (Tx):: sss is equal to the leverage score of row i, and the sum of leverage scores is less than or equal d.
We can use the same theorem to characterize the regression error for t by changing the sign of (0).
4.2 Full Observation ATE
We now sketch proof of our result for ATE estimation with full observation. Proof sketch of Theorem Here, we only sketch our proof for the variance bound. For S S [n]] and
2 E { (11, +1 (where 7 {i:
+1} and 7 (i 2 -1)), we define
= [s| Eies( (1) yi (0)), and is 2
(EEzztn 5. (2.1) EiEz-ns v/200)
[s|
Similarly, define 73 and Pz. Then for 7 as defined in A Algorithm 1
we have,
E, [Es [[ -] Ex
Es
[S] 151-7s+131-7
1S [3
2
EEs
1 n Ez (7s 21
1ss nd
Ez [(fs rs)
+1s 2.
EE [(Ts
Ts) - Ts)
n
By Sauchy-Schward inequality, Ez Pf Tz) (Fs [s)
S E2 [ 13) )2 E2 [ Ts)
Therefore we only need to bound E. [ }3 and E. [
)=) We can bound
E, [(fs rs) by Theorem 7 (ridge leverage score sampling) and using the fact that is is ob-
tained by (2.1) and (22.0) which are adjusted by (2)), vector that is learned from the units in the set S. Since S and are disjoint, ((2) is independent of 2, and this allows us bound E, [(Ts Ts) )s)²
)²]
The bound on Ez [ - ?] follows from similar argument.
Experiments
In this section, we compare our method for estimating ATE and ITE with full observation with prior art and baselines on synthetic and real-world datasets. For ATE estimation, our experiments
demonstrate that our approach either outperforms other unbiased estimation approaches has comparable performance terms of the achieved variance.
8

---

Algorithm 1: ATE estimation with leverage score sampling and cross adjustment
Input: X E Rnxd, y(0), y E IR", r, E (0, 1]. 2 Set y E IR" random vector such that for each E [n]. yi independently and with equal
probability either equal to y y (t)
3 Set - to the set of entries i such that ys is equal to yj (1) and z = [x]] \ Z+
4 Set A = - < lo)(n/ where Ç = max ||x: 12-
s Set 5 C [n] to random set of indices where each i E [n] is included in S independently and
with probability 0.5. Set S [n|] SS.
6 Set to an n X n diagonal matrix corresponding to S where each entry S. = 2 if E S, and
S. = 2, otherwise. Similarly, set S.
7 Set (1) = arg min |SXb- Syll + A | 1 60 = arg min SSb - 5y + X |b]] 2
8 Set (1100) = y0) Xb(1) ((11)) -y(1) - (1) -((200 =y(0) (0) xb -y(1) -x(2))
Return == 2 n
(Eiezin of'1)) Eieztns Ui vf Eiez nB U EiEz ns Dj of2001).
ATE methods. Our ATE estimation method is compared against five other approaches. (i) HT Uniform: This employs the Hurvitz Thompson estimator, using uniform assignment of units to treatment and control groups. (ii) GSW: This design follows the methodology in [8]. (iii) Classic Regression Adjustment: As described in (24), this estimator is inherently biased. (iv) Leverage Score-Based: In this method, we learn two distinct vectors for y()) and (() through leverage score sampling. The difference between these vectors predicts the treatment effect. Though this method has its advantages. it is biased and was previously employed by Addanki et al. M. (v) 4 Vectors: Obtained by specializing Mou et al. [28] to linear regressors, this approach adopts the same cross-adjustment mechanism to mitigate bias as in this paper. However, instead of random vector regression, it determines separate vectors for (1) and y!) for each group, leading to four vectors obtained by linear regressions.
ITE methods. We discuss two distinct approaches for ITE methods. (i) Baseline: The ideal scenario or "oracle" baseline involves utilizing a linear function, trained on vector t. Naturally, this is not feasible in practical applications. (ii) leverage score based: This is based on learning two different vectors by ridge leverage score sampling and using the difference of the two linear functions for estimating t. This is an adaptation of the algorithm in D to the full observation setting - we make this adaptation because the exact algorithm of [1 cannot be applied to the full-observation setting due to sampling with replacement. We do not compare to 28 because their guarantees are only for an infinite population.
Synthetic Datasets. For synthetic datasets, we will construct y(),, .y(()) such that (for the ATE case) or t (for the ITE case) is a linear function of the covariate matrix.
ATE Dataset. There are 50 covariates (i.e., features). Each entry of X E Raad is a uniform random number in J0, 0.01]. Each entry of b E 55 is uniform random number in [0, 1]. The individual treatment effect vector f1 is Xb+r, where each entry of E IR" is randomly picked from a mean zero Gaussian distribution with standard deviation of 0.2. Each entry of (0) is picked + μ and t t=y (0) The results for this dataset with different number of samples is shown in the right plot of Figure [ Our experiments illustrate that on the synthetic dataset, our approach outperforms all other approaches except the GSW design. However, we note that our approach is computationally more efficient than GSW design (see T and gives much simpler design that can be used in different settings.
uniformly at random from [0, 5). Then y = -
2. ITE Dataset. The covariate matrix X and the vectors b, r, y are picked similar to the ATE dataset. Then individual treatment effect vector is Xb + r and y (1) = y((0) + t. The results for this dataset with different number of samples is shown in the left plot of Figure I ITE error for our method is consistently smaller than leverage score based method.
WoaldWorld Datasets We analyze the following three distinct real-world datasets. A comprehensive breakdown, including variance and other measurements for each dataset, is available in Table in the appendix.

---

leverge besed
Ht u Kor
MM5 base ine
dasac regression mojur
loverage score besed pecion
0.06
11 1500
2/22
3020
3.05
Population 17500
2230
Population Size
Figure 1: Results for synthetic ITE and ATE datasets are shown on the left and right, respectively. For different population sizes (i.e., n), estimation is performed for 1000 trials. Then the average of ||t-ū|a vn-I[t]2 over these trials is shown with solid lines for ITE and the shades around these lines denote the standard deviation. For ATE, the solid lines represent the average of | Tl/|T| over 1000 trials.
1. Boston Dataset 17 This is dataset of housing prices in the Boston area, consisting of 506 samples and 13 features. Since it has only one label, we set y = =y(") i.c., ATE is zero.
2. IHDP Dataset [16, [I]: Derived from the characteristics of children and their mothers, this
dataset comprises 747 samples and 25 features. ATE for is -4.016.
3. Twins Dataset (3): This dataset is constructed based on the characteristics and mortality rates of twin births in the US. We specifically selected samples that have complete feature values, resulting in dataset of 32,120 samples with 50 features. The ATE of this dataset is 0.0064.
Table 1: Results of ATE estimation. For each result, the first number the average of |r 7 over
1000 trials and the second number is the standard deviation of this quantity.
Dataset
Uniform HT
GSW
Classic Reg Adj
Lev Score
4 vecs
Ours
1.736
0.663
0.333
0.658
1.677
0.628
Boston
1.339
+0.510
+0.255
#0.504
+1.256
+0.459
0.272
0.042
0.012
0.536
0.264
0.040
IHDP
+0.206
+0.031
#0.009
+0.050
+0.203
+0.030
1.351e-3
1.231e-3
1.201e-3
1.226e-3
1.369e
1.218e-3
Twins
+1.025e-3
+0.937e-3
+0.899e-3
10.91le-3
+1.015e-3
+0.936e-3
Table 2: Results of ITE estimation. For each result, the first number is the average of J
1||2
over 1000 trials and the second number is the standard deviation of this quantity. The relative error
columns report <
2 |t|l2-
Dataset
Baseline
Leverage Score
Ours
Leverage Score Relative
Ours Relative
Boston
0.0
0.827 0.469
7.654
1.678
IHDP
0.548
2.449 0.071
1.765
0.223
0.021 0.637e-3
0.015 1.9e-
Twins
0.155
0.156 + 2.66e-5
0.156 2.21ee4
5.57le 9.5e-7
5.57 7.9e-
6 Conclusion
In this paper, considered the problems of mean estimation, ATE and ITE estimations in the presence of covariates. We considered the finite population setting and provided non-as con-asymptotic variance bounds for several novel variants of the classical regression adjustment method-based estimators. Our guarantees are model-free, even the covariates are arbitrary (are not informative of the estimand and may be chosen adversarially), the variance bounds still meaningful and almost match the variance bounds of widely used estimators not based on covariates.
Our algorithms are simple and efficient, and we believe they readily extend to many related settings like arbitrary assignment probabilities and online treatment assignments. We also believe the results can be stated for the kernel setting. These extensions will be part of future work.
10

---

References
[1] Raghavendra Addanki, David Arbour, Tung Mai, Cameron Musco, and Anup Rao. Sample constrained treatment effect estimation. In Advances in Neural Information Processing Systems, 2022.
[2] Ahmed Alaoui and Michael W Mahoney. Fast randomized kernel ridge regression with statistical guarantees. Advances in neural information processing systems, 28, 2015.
[3] Douglas Almond, Kenneth Y Chay, and David S Lee. The costs of low birth weight. The Quarterly Journal of Economics, 120(3):1011-1083 2005.
[4] David Arbour, Drew Dimmery, Tung Mai, and Anup Rao. Online balanced experimental design. In International Conference on Machine Learning, pages 844-864. PMLR, 2022.
[5] David Arbour, Drew Dimmery, and Anup Rao. Efficient balanced treatment assignments for experimentation. In International Conference on Artificial Intelligence and Statistics, pages 3070-3078. PMLR, 2021.
[6] Haim Avron, Kenneth L Clarkson, and David P Woodruff. Sharper bounds for regularized data fitting, Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, 2017.
[7] Adam Bloniarz, Hanzhong Liu, Cun-Hui Zhang, Jasjeet S Sekhon, and Bin Yu. Lasso adjust- ments treatment effect estimates in randomized experiments. Proceedings of the National Academy of Sciences, 1 113(27):7833-739) 2016.
[8] Xue Chen and Eric Price. Active regression via linear-sample sparsification. Conference on Learning Theory, pages 663-695. PMLR, 2019.
[9] William G Cochran. Sampling techniques. John Wiley & Sons, 1977.
[10] Michael B Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, and Aaron Sidford. Uniform sampling for matrix approximation. Proceedings of the 2015 Conference on Innovations in Theoretical Computer Science, pages 18 181-190, 2015.
[11] Vincent Dorie. NPCI: non-parametrics for causal inference. https: https:///github..cm//vvoo//
npci 2016.
[12] Matthew Fahrbach, Gang Fu, and Mehrdad Ghadiri. Subquadratic kronecker regression with applications to tensor decomposition. In Advances in Neural Information Processing Systems, 2022.
[13] David A Freedman. On regression adjustments to experimental data. Advances in Applied Mathematics, 40(2):180-193, 2008.
[14] Alan S Gerber and Donald Green. Field experiments and natural experiments, 2008
[15] Robert Greevy, Bo Lu, Jeffrey H Silber, and Paul Rosenbaum. Optimal multivariate matching hefore randomization. Biostatistics, 5(2):663-275 2004.
[16] Ruth T Gross al. Infant health and development program (ihdp): Enhancing the outcomes of low birth weight, premature infants in the united states, 1985-1988, 1993.
[17] David Harrison Jr and Daniel L Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of environmental economics and management, 5(1)811-10,, 1978.
[18] Christopher Harshaw, Fredrik Sävje, Daniel Spielman, and Peng Zhang. Balancing covariates in randomized experiments with the gram-schmidt walk design. arXiv preprint arXiv:1911.03071,
2019.
[19] Daniel G Horvitz and Donovan Thompson. generalization of sampling without replacement from a finite universe. Journal the American statistical Association, 47(260);663-68) 1952.
11

---

[20] Andrew Jesson, Panagiotis Tigas, Joost van Amersfoort, Andreas Kirsch, Uri Shalit, and Yarin Gal. Causal-bald: Deep bayesian active learning of outcomes to infer treatment- from observational data. Advances in Neural Information Processing Systems, 34:30465-30478, 2021.
[21] Nathan Kallus. Optimal priori balance the design of controlled experiments. Journal of the Royal Statistical Society. Series B (Statistical Methodology), (0(1):85-112, 2018.
[22] Edward H Kennedy. Towards optimal doubly robust estimation of heterogeneous causal effects. arXiv preprini arXiv:2004. 14497, 2020.
[23] Xinran Li and Peng Ding. Rerandomization and regression adjustment. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82(1):241-268, 2020.
[24] Winston Lin. Agnostic notes on regression adjustments to experimental data: Reexamining freedman's ss critique. THE ANNALS of APPLIED STATISTICS, pages 295-318, 2013.
[25] Ping Ma. Michael Mahoney, and Bin Yu. A statistical perspective algorithmic leveraging. In International Conference on Machine Learning, pages 91-99. PMLR, 2014.
[26] Michael W Mahoney et al. Randomized algorithms for matrices and data. Foundations and Trendy in Machine Learning, 3(2): 2011.
[27] Kari Lock Morgan and Donald Rubin. Rerandomization to improve covariate balance in experiments. The Annals of Statistics, pages 1263-1282, 2012.
[28] Wenlong Mou, Martin J Wainwright, and Peter I. Bartlett. Off-policy estimation of linear func- tionals: Non-asymptotic theory for semi- cemi-parametric efficieec. arXiv preprint arXiv:2209. 13075, 2022.
[29] Akanksha Negi and Jeffrey M Wooldridge. Revisiting regression adjustment in experiments with heterogeneous treatment effects. Econometric Reviews, 40(5):504-533 -344, 2021.
[30] Tian Qin, Tian-Zuo Wang, and Zhi-Hua Zhou. Budgeted heterogeneous treatment effect estimation. In International Conference Machine Learning, pages 8693-8702. PMLR, 2021.
[31] Mark Rudelson and Roman Vershynin. Hanson-wright inequality and sub-gaussian concentra- tion. Electronic Communications in Probability, 18:1-9, 2013.
[32] liris Sundin, Peter Schulam, Eero Siivola, Aki Vehtari, Suchi Saria, and Samuel Kaski. Active learning for meaision from imbalanced observational data. International conference on machine learning, pages 6046-6055. PMLR, 2019.
[33] Stefan Wager, Wenfei Du Jonathan Taylor, and Robert Tibshirani. High- -ighh-immnsional regres- sion adjustments randomized experiments. Proceedings of the National Academy of Sciences, 1 113(45):12633-124) 2016.
[34] Martin B Wilk. The randomization analysis of generalized randomized block design. Biometrika, 42(1/2):70-79, 1955.
[35] David P Woodruff et al. Sketching as a tool for numerical linear algebra. Foundations and Trends@ in Theoretical Computer Science, 10(1-2):1-157 2014.
12

---

A Warm-up: Mean Estimation
In this section, we consider the problem of estimating the mean of a single vector y by observing only a few entries (i.e., y == Ei=1 Vi). The most classic approach is to simply uniformly sample m units from the population and report their average. This results in variance of at most _ lly v1 where y is vector of the same size y with all of its entries equal to y-
mn
The main theme of regression adjustment is to replace the dependence of variance on ly y 2 with a dependence on the error the best linear fit, i.e. millber ||y y (X X)b where X is vector of the same size X with all of its rows equal to X 1 Σ i=1 X, We wish to do so without making any modeling assumptions - a common modeling assumption in statistics is that the response vector is a linear function of the covariates plus Gaussian noise. However, our goal is to achieve significant variance reduction if y - y can be well approximated by a vector in the span of X- - X without any such modeling assumptions. In the worst case, when y - y orthogonal to the lly yl which recovers the result
column span of (X - X), mimbeR" y y - (X - X)ble
(Xbl 2
of the classic approach.
For the mean estimation problem, if an oracle provided b* = arg minDER | - y x)b|l the estimator ÿ [S] Eies Di (x; x) x)=b* is unbiased with its variance bounded by
mimbERa ly y - (X X)b| where S is a uniformly sampled subset of the population of size This is classic result (see Chapter 7 of (9]) that also follows from the following Lemma.
Lemma (9- Chapter 2). Let y E R2 and y be the average of uniformly sampled set of entries of y of size m. Then
ly - y 2
E [0
m n
m(n
where y is a vector of size equal to y for which all the entries are equal to y.
One can easily see that by applying Lemma y, where y: y, (x; x) b, a variance of Ily 1 (1 can be obtained. m(n-1)
Key challenge. The main challenge for the mean estimation problem to compute by (or a vector close to it) without observing all the entries of vector y. Note that do not have access to any entry of y - y because we cannot compute ÿ without observing all the entries of y.
We resolve the issue of not being able to compute y by performing the regression a modified matrix. We add a column of all ones to the matrix X X and perform a regression on using the modified matrix see the proof of Theorem 10 More specifically, we prove the following.
min be
X - b (y
min bER-H1
(X - X) 1] b xy1
To solve the regression problem on the right-hand side without observing all the entries y, we use leverage score sampling 25 [TO]. Leverage score sampling is numerical linear algebraic tool that
allows for solving a linear regression problem approximately by subsampling about d units from the
population and solving regression problem over these units, where d is the number of features [ see Section
3 for details. Combining these techniques gives the following results.
Theorem 10 (Leverage score sampling based regression adjustment). Let e, E (0,1) and m E N and X E R be a covariate matrix. Then, Algorithm 2 computes an unbiased estimator of the mean of an outcome vector (on a population of size n.) with ((dlog(d/8)fe + m) samples. Moreover, there an event E with Pr(E) 2 1-8 8 (over the randomness of the algorithm) such that the variance of the estimator conditioned to E is bounded by
(1 + 6) m(n - 1) (1 m n
min ||X
X)b (y -
Note that since the dependence of number of samples on d is logarithmic, we can pick d polynomially small, e.g., 8 = 1/n?. The variance upper bound achieved by our method is compared to other methods in Table 3
13

---

Table Variance of different methods for estimating the mean population with few samples. denotes the maximum norm over rows of X, 0 < & < 1 is an input parameter of GS walk design, and y is vector where all entries are equal the mean y. Variance
Method
Reference
Uniform sampling
m 72
v ly
((m) samples
m (n 1)
GS walk design
(by adjusting
[ ||(x
(y y)l1 - σ² " |bbl
probabilities)
m
min beR
X)bb
1-6"
18
((m) samples
Two-phase
O(m 0(mlog(mm88//4
(1 +) m
1
min bERd
m -
||(x - X)b (y --)ll <² [|b]
This paper
samples
Two-phase
O(m + dllo(d/8)/e{) (1+c)
min beR
X)bb (y - x>>1 [X
This paper
m(n 1)
samples
Regression adjustment. Our estimators are based classic technique in statisties called regres- sion adjustment [9] Given vector b (either given by an oracle or learned from the data), the classical technique uses the sample covariates to "adjust" the estimate of the sample mean |8t| Σies and 1
outputs the following estimate
1
Σ -x -x-xx)-.
S
If S' is a random sample independent of b, the adjusted estimator is unbiased. The adjusted estimate will have lower variance than the sample mean estimator if |ly y x lyy 2
Lemma 11. Lety R and X Raxd, Let b be a fixedhreassigned vector. Then the variance of b, where S' a uniformly sampled subset of [ni] of size m,
1 Eies Us- (x; x)
m(n 1) (1 (1-m) m))
1
- (X X)b
Proof. Since b is fixed, by linearity of expectation
E
E
Σ
[x
x;) T
= + (x
Σ
x;) Th
T = yy
S
S
'ES'
"ES'
$ES"
Therefore by Lemma (applied to the quantities Vi + (x - xx) (b)
7 - (xi - x)
Tb
E
[ 3)
1
(1
m
71
and the result follows.
Therefore our goal is to find b that minimizes |yy - - (X -X Our approach is to use leverage score sampling (Theorem 7 to find vector that approximately minimizes this error. A main technique in our approach that the samples picked for leverage score sampling are independent of the samples drawn for estimation. More specifically, our estimator is based on two sets of samples.
We of samples S to estimate b and independent of S) in the following estimator ÿ |St Eies Vi + (x xx) b. By Lemma this unbiased and achieves variance of
the another set S' of uniform samples (that
m (1 (1-m) n
| x X X- (y 8)
m(n 1)
Although it straightforward to carry this independence for the mean estimation problem, as we will see, we require more involved algorithm provide it for the treatment effect estimation problems.
14

---

Algorithm 2: (Ridge) leverage score sampling based regression adjustment
Input: X E Raxd y E R", m E N, 0 < <<< <1 For Theorem 10 set = O, and for Theorem
n/m Compute the leverage scores c!A) of X [x -) 1].
12 set
x² =
11 5 m |xel2
Set p; [1,3 EfA) log(d/0)/e* for E [n].
Let S E R be random diagonal matrix with S. 1/ /pi with probability pi, and
S = 0, otherwise.
~ Let b arg bin SXb - Sy 3y|2 + A || 12 and b E Rd be a vector consisting of
the first d entries of b.
Let S' be uniformly sampled subset of size of [n] and f = Return 7.
1 m Eies Vi (x (x-x;) Tb
Theorem (Leverage score sampling based regression adjustment). Let e, E (0, 1], and m E N and X Bnx be a covariate matrix. Then, A computes an unbiased estimator of the mean of an outcome vector (on population of size n) with Oodlog(d/ (/ + m) samples. Moreover there an event & with Pr(E) 2 1 (over the randomness of the algorithm) such that the variance of the estimator conditioned to & is bounded by
(1 m(n 1) (1 (1-m)
min ||X ([XXX))b (y -))1
Proof. Let x = [X- X)
1].
We first show that
min be
X
X)b (y - --)1
min bER
Xb $
(3)
2 y
Note that minbeRd-1
Xb
| y 1 Since 1 is orthogonal to X X.
x
x TTy.
Moreover TTy y. Since 1 is orthogonal to X. y y is orthogonal X X, and = 0. Therefore
Tx== x (y y) +y.
Subtracting both sides by y and taking the norm, we have
x v11 x (y y (y y
min beR
X - X)b (y - -
which implies our claim. Now note that by sampling with high probability, we can find b such that
entries of y. by Theorem 7
Xb v13 < (1++))
x113 = (1 + t) min
min beR
Xb
||X (IX-X)--- (y -))
beR
Let Rdd be a vector obtained by taking the first d entries of b. Then, we have
Xb - y = (x - -x6 + bat1 1.
Now consider the following regression problem:
1 . t - - -x (x-X)b) 2
min tR
-
The optimal solution is 1 1 (y (X -x)b) = 1 y. Therefore
Xb y >/2 Xx x)b y + Tost 1 11/ X x)b y k)b-y+3 +y 2
(5)
15

---

Thus by 44 and (
| X]b (y 111 (1
min bERd
| - X)b (y- -))l
and the result follows by Lemma
Note that the number of samples in Theorem T0 depends on the number of covariates (also called features). This dependence can be alleviated by using ridge leverage score sampling at the expense of higher variance bound. If we want to sample < dlog(d/8)fe² sampling, which gives a bound on the ridge regression error instead of the regression error. Using X-ridge leverage scores, for À <² where & is the maximum &g norm over rows of X. gives the following result. Before providing this result, note that the block leverage score sampling approach of [12] shows that to solve
use ridge leverage score
min BER
Xb - 11b--112-1 + /
we can use the same ridge leverage sampling as for mimbeRd | y122 + and the number of samples only needs to be increased by constant factor. We require this for ridge leverage score sampling-base baased mean estimation result.
Theorem 12 (Ridge leverage score sampling based regression adjustment). Let <<<<<<<o<
0.5, E N, and X be a covariate mairix. Then Algorithm 2 computes an unbiased estimator of the mean of an outcome vector y (on a population of size n) by using O(m log(m/8) (mmmmm/8/)))) samples. Moreover, there is an event & with Pr(E) 2 1 - 8 (over the randomness of the algorithm) such that the variance of the estimator conditioned on E is bounded by
[|[X K)b (y - --)1112 + 11b12
(1 + c) m² beRd
min
(6)
Proof. Let x = [XX X)
1]. and
[X X)
(X - X) 1 ViL
Z
Al
Matrix Z corresponds to the regression problem
y11 + b1
min
Xb
bER
In other words, defining the vector y E Rn++d as a vector that is equal to y on the first entries and is zero on other entries, we have
Xb - y 2 + +A +blbbadle
Zb
Zb. 2b-y1 151 is
for any vector b E Rd+11 Now note that the solution to arg min bER
7 ブ
2 Z)
TTy y,
by definition of y. Moreover by definition,
I 0 +
z X
X
0
Therefore since 1T(X - X)= = OT, we have
ZT X X
(X-XX
+
AL
0
0
Since this is block diagonal matrix
Z -1
X
+
ALa) -1
0
T
7
16

---

Therefore
Z Z)
Thus
(X- (x-X) (X- XX + Aaa) -1 XX-X) y+ +
y - -
1 - 8112=
x
VdIa
0g
Note that (x-X) XX - X) +l1) -1 (X xx))((x y is the solution to
[|(X-X]0 - (y - 8)1 * |b|
arg min
bERd
and because 1 in the null space of X-X,
((X - (x-x) X) + Ala) (x-X) T (y -y) = =0(=(()--) -(X-Xx (0x-x)((2 (X X) + Aaa) (X-X)] xxx
Therefore
|| - 11 |(T2= - I) y y
0₂ 1
and
min BE
Xjb (y 3)11 -y)1 [
+ A b
min bE
Xb
x|| + b 2
(7)
Now note that by ridge leverage score sampling, we find a vector b E Rd+1 such that, with high probability,
Xb 111 + A |ind (1 +e)
xb y1 + A Albbaal
min beR
Therefore
Xb- 11( + 2 bi
||= s (1+ (1++)
2
[|(x- - X)b (y - --) + X [bb]
min beRd
(8)
Now let E Rd be the vector that is equal to the first entries of b. Then we have
xx-xx6+baar - 11 6 |
Xb
Now consider a regression problem of the following form.
-(x-(x-- 1 1 - (y - (X -xxb) 2
min ER
2
The optimal solution to this problem is 11T ((y - (X - Xib)) = 1 = 7. Therefore
Xb- - > (x
X)b - y + b041
2
>
X - x -
(9)
2
Thus by 8 and (9) and because > |16/|² 2 0,
| 10-x)6-oe (y (y- 1) S (1++)) e)) E) min
[|[X X-b (x 20122 + |b|2
->>
ber
Therefore by Lemma the variance of the regression-adjiue than.
estimator with m samples is less
()+ min beRd
-bb (y
[[[X-
+ +
mn
Finally by definition 2_ Cnn/ S²
2 and the result follows.
17

---

We now explain previous approaches to mean estimation. The classical approach to regression adjustment [9] is to learn vector b from the same set of samples that are used for estimation.
= arg min lys - -ys - (Xs X b
The Gram- -chhmih walk design of Harshaw al. 18 has originally been developed for the average treatment effect estimation and achieves the following bound using the Horvitz- Thompson estimator see Algorithm 5
Theorem 13 (Gram-Schmidt walk design (8). Let O
1. Given y X E R the Horvitz-Thompsor estimator under the Gram-Schmidt walk design is umbiased and has variance of
R and
1 [
E [
7)
S
min
Xb
1112
+ 1
|
n BER"
where μ y() + y()))
special setting of ATE estimation problem reduces mean estimation. Suppose y(1) y and y()) = 0. Then the average treatment effect for this problem equals the mean of y. This approach partitions the samples into two groups. The first outcome is observed for the first group, and the second outcome is observed for the second group. The estimate is then produced by the Horvitz- Thompson estimator. Note that in our special case, only the samples in the first group amount to observation since y()) is the fixed vector of all zeros. This yields the following guarantee for mean estimation
Theorem 14 Schamdt walk based design 18).
The Horvitz- Prriit-Thoompon estimator applied to the output of Algorithm 3 gives an unbiased estimator of the mean of y with a variance of less than or equal to
min beRn
[|[X-
X-b (y --)1 1B11
(10)
m²
with (m)) samples.
In comparison, the bound of Theorem 12 is better by a factor of m/n on the |Xb yllă term.
B Individual Treatment Effect Estimation
In this section, we give an ITE estimation algorithm with partial observation (Algorithm 3). Our approach is to learn a vector b that nearly optimally fits t := y()) y()) linearly onto X, i.e., ||b - t/|2 is small. Since do not have access to t/2, we perform regression
random vector for which the expectation is t.. To this end, we construct the random vector by setting y; independently
(1)
(0) with equal probability. If we observe all samples, and use Theorem 8
to either y; -y;
vectors (1) and (0), this approach gives the following expected error bound.
E |2xx
= < d a|py
y y 1|
+ min |Xb t
We can achieve similar bound using leverage score sampling by thresholding the probability of sampling rows from below while we guarantee that the number of samples does not increase too much with high probability. The following is the main technical contribution of this part.
Theorem 1 (ITE Estimation). For any 6, a E (0, 1]. there exists randomized algorithm (Algorithm that observes potential outcome for (d (ddoogd//)// + n) units and outputs vector t.
Moreover there is an event such that Pr(E)= 2 (over the randomness of the algorithm) and
tl|21]] $ (1 c1) minDER [Xb t + _ (1+ (1+)))1100)+00000
=01011
Proof.
Note that 2y + z μ. Then
2x6 = 2X SSX) - 1x
= -X(X SSX) 'x SSt + X(X SSX) -ex')-x'' SS(2 o pt)
18

---

Algorithm 3: ITE estimation with random vector regression and leverage score sampling
Input: X E guxd , E RR", 0 < E 1,0 <
Compute the leverage scores of X.
Set p; min maxfa, 3 4 (r((d/))/=)) for E [n].
Let S E _ be random diagonal matrix with Sn 1/ /Pi with probability ps, and
S 0, otherwise.
Let b = arg min |SXb - Syll2-
Return t == 2Xb.
Moreover
X SSX) 'x SSt = arg min SXX St bER
Therefore by Theorem
7 and since | tll ming
|Xb tll
X SSX) 'X SSt tl = + <) || ||xt - tll2
Now note that
||X(X] SSX) X X C 14) [3] - (1+ c)E (z p)) SSX(X SSX) 'x SS(z + p1)
E
(SX) 2
=11+e) Σ
Pr
iSu
5 1 -1 (e+ |ull
Σ (SX)
E:S,
< (1 + +) |1
where the last inequality follows from the fact that the sum of leverage scores is less than d, and the penultimate inequality follows from the definition of p;. Finally, the result follows by noting that,
[XX(X - SSt t) T X(X SSX) - TSS( n)] 0.
C ATE with Full Observation
Our approach for ATE with full observation is similar ITE in the sense that we perform regression on a random vector to learn solution vector b. Then, for estimation, we apply regression adjustment based on b to adjust the outcomes. We first present naive version of this approach that is biased and achieves factor proportional d2 /n² in one the terms bounding the variance. We improve this to d/n² by partitioning the units into two groups and using the learned vector for each group to adjust the estimation in the other group. With this approach, our estimator unbiased.
Before discussing our results, we discuss the other approaches in the literature. The Gram-Schmer walk design of [18] that uses the Horvitz-Thine Thompson estimator gives the following variance for the ATE estimation.
1
min bER"
Xb - 11l2 +
bb
n
where 0 < is input parameter. The classic regression adjustment with interaction term gives an asymptotic variance proportional to 1 (of + 20 + σ11u), where of is the variance of
V = v1 (1) ((()) +X; ()) (1) (24, where ((1) is the best linear fit for ()) -2(2) using X when
124). Similarly, of is the variance of ((0)
y (0) 7( (X, ()) and 0
the covariance of y. (()) Note that the above mentioned bound is achieved by learning two
separate vectors for regression over (1) (1)) and y; (0) (() Moreover, that because these bounds are asymptotic, we need to assume a model data generation.
For our naive approach, first need the following technical lemma.
19

---

Lemma 15. Let X E Raxd, T E Ruxn be a diagonal matrix and E {-1 be random vector where each entry is picked independently and uniformly at random. Then
E
2
TTx
Tx Tz) - (d + 2d) l|t|l
where t E R is the vector corresponding to the diagonal matrix T.
Proof. Let M R and. Then
E 2 Mzz M
=
44 M3y
E
22 M. MA
M,
Note that for any term such that there exists one of the i, k, that is not equal to any of the other three, the expectation is equal to zero. Therefore we need to look at terms where for each i, j. k, I, at least of the other three is equal it. Therefore have four cases:
1. --=== = 1. This gives an expectation of M So total from these terms is Σ M
2. i j and k = but # k. This gives expectation of M. So total from these terms is
2 E Ek=++1 M&
3. i = j. = and i + j. This gives an expectation of M So total from these terms is
2 2 2) M
4. 1 1 k, i # j. This gives expectation of M, M, So total from these terms is
2 2' Ej M,,
Recall that we denote the leverage scores of matrix X with 4. Therefore taking the above over our matrix M Tzx, since Mij- = (TTx)ij, we have
E
Mzz M z)
(xx()))(()xx
n-
+
(TT))
|=
+2
Σ (Tx)
J=H1
$ 11t11
2.) 2 +2 1
< | (d + 2
< ||t|| (de + =d)
where the last two inequalities follow from the fact that the sum of leverage scores (which are the diagonal entries of the projection matrix σx) is at most the rank of the matrix.
are now equipped analyze the ATE estimation obtained by purely using random vector regression and regression adjustment.
20

---

Theorem 16 (ATE with full observation). Let y E E" and y E R" be a random vector such that for each i [n], Vi is independently and with equal probability either equal to yi or Vi Moreover let Z+ be the set of entries i such that yi is equal to yi (1) and Z = [n] \ Zt. Let
(0)
b arg minl |Xb -yll Let y y( Xb and 5( y(1) Xb. Let 7 be the following Horvitz-Thooe Horvitz-Thompson estimator, i.e,
Σ Di
Σ
(0))
V/
iEZ+
iE
Then 1 S d (1l
-y(0) |00 and
d2 + 2d
((()= x 12
E [f - -)
+ n²
min Xb -m112
n²
where μ y() + y!0).
Proof. Let {-1,+1}"* where 2, +1 i Z+ and 2, -1, otherwise. Moreover let
Ā = ÿ (()) + ÿ(1)) = =y(0) + y 2Xb. First, note that 2y = I + 2 о Therefore
4= 22-6- μ 2xxx (x + (μ πx,) (x
Moreover
8 - 1
Σ
(20)
Vi + +i
(0))
(2 + Vi Vi w;
(0)
(1)
(11)
VEZI
IE
(u y;
(1) (0) 2x6)
Σ
(0) + y
(1)
(0) 2x b)
n
ViEz
1
T
n 4 z.
Therefore
E
7|
EE - [ z - (z 2 t) TXZ
(12)
Thus by the linearity of expectation and since is equal to +1 and -1 with an equal probability. (1 xxp) z] Let T be diagonal matrix associated with t. Then t Tz. Therefore
1 n
E F
TTx2]
Therefore since 2; and 2j are independent for + j. and 23 with probability one.
n
n
E B
4 (TX) )i
tel
Therefore by the triangle inequality
||||
E F - r||
Σ n
t.l
Σ
Then the first part of the theorem follows by the fact that the sum of leverage scores is less than d, and noting that by definition, t. =y =y - (0)
We now prove the second part of the theorem. First, note that since 2; and zy are independent for i # j and 2, = +1 with equal probability,
E [H T zz TTxz] E [(me TxH) TTxz] 0.
21

---

Therefore by 11 and Lemma 15
E
[f - -)
= n²
E
(z) T
- -x - TxpL 22]
1
+
E
2 TTxzz TX
n²
S n? (1 (1le
1
xx/1 (dP 2d) ||tl(²2)
The next step is two use subset of samples for regression and the rest of the samples for estimating the ATE.
Since we want to have a guarantee on the regression error, we need to make sure that the sampling (i.e., partitioning) to two groups is done according leverage scores. However, if the leverage score of a row is greater than half, then it cannot have a probability larger than its leverage score for being assigned to each group. To reduce this probability, we use the ridge leverage score. The following lemma states that by picking a large enough ridge parameter. we can achieve a small enough ridge leverage score for all rows (therefore, small enough sampling probabilities). This is at the expense of a regularization term of the form A [|b|² in the regression error.
Theorem 6. Let Ç := [ ||x|| Then for X 2 c. f2 2 1. A (X) 1 for all i [n].
Proof. Recall that Ç == ma mn) ||x l22- Therefore, ||x;|| 1 5² for all i E [n]. We have
(x) = x X x + -I) S x; (AI) -x, = ||x|| 1/c.
Equipped with this, we present the main result of this section.
Theorem (ATE Estimation - Full Observation). For any e, 6 E (0, 1), there exists a randomized algorithm (Algorithm) 7p that computes an unbiased estimate i of the ATE r, and for which there is an event E with Pr(E) 2 1 - (over the randomness of the algorithm) such that
E [(f - - +[]+c) minDER" (IXXb Hal2 + 100log(m/5) 5² (1b112) + 32d 322.111y -ly y 120.
where Ç = max ||x: and sa ===; y()) + y y() is the total outcome vector.
Proof. For S C [n] and z E (-1,+1] (where Z = fi: 2 = +1} and Z = fi: = -1)), we
define
(0) and Ps =
2
Σ
(2,1)
Σ
y ((200)
Ts =
S
S
IES
NEZ ns
iEZ nn
Similarly, define T5.
Then we have,
Ez
ES [f
1)
Eg
Es
|S|
Pg + |S| -
|S| Ts +
|S|
Ts 22
n
n
|SS
= EE
Ex [TS [3
n
|S|2
+ Es
E2 [(fs - Ts)
n2
|S| |S|
Es
2.
Ex [(Fg Tg) (Ts [
n²
We now bound all of the terms above. Let
(1) := (1.0) + (1,1)) = (() +y(1) - (1)
ju (2) = ((2.0) + ((2.1) = (())) + +y - 2x6(2)
22

---

Note that
Σ
(20)
(1,1)
Vi
(1) (0)
(1.0)
+ yi y3
(1) (0)
[S|
iEZ+nS
iEZ ns
1
Σ
(Ui
(2)
+ yi
(0)
2x, (x]T
Σ
(vi 9/
(1) (0)
2x, B())
[S]
iEE
n0
1
(mgg)
Zg.
[S]
Similarly, have
(15 ( T
Gs - Ts |5|
Therefore
|||
[S]
Eg
[3) (Ts - rs)]
I Ex
(TLS
(2)
28
n
n²a
<
n²
EEz
(1)
( z3
E8
S
zs}
1/2
where the inequality follows from Caucch-Schward inequality. Moreover
S
Ex
[Tg - -g) = n?
Ez
(
(1)
25
and
n
S Eg [(Fs - rs)
E2
(iz x
(2))
n
n
Therefore we only need to bound Ez 2 (ps (11 5 are disjoint, 5 is independent of Zg. Moreover, the entries of independent from each other. Therefore, since 2y μ t
(2)
(2)
Tzs and Eg - s
(1) ( ) zg
Since S and
Ex
2
S
2s
= Ez
Σ Lie
[l(ce x(x
SSX + -TT-sSXTSXxeX + 20))s 2
Ez
where x = 100ç² lon(n). We have
E2
[1(1=
-
x--XXXXSSSX SSX + A. -
1 SS 2 tt)s
< 2E,
SSx + X 1) SSp)
2
||[X(X
sSx + A X SST2)
Note that due to guarantees of leverage score sampling.
[1(u X(x SSX - sSp)s < 1| X(X (IXTSSXX SSX + - SSp
$ (1 + e) (|Xb - اlll +
Moreover
Ey [|( TSX + +X-I--IX SST2)
< Ex ||x[X SSX + A - T SST2|12]
Since
x (1+ SSX + I
23

---

we have
Ex [|X(X SSx + A -I' SSTT2|2 3 Ez [z TSSXXXX SSX + -I- SSTz]
(Ix(XSSXX+
(SX(X SSX + + - -
SES
1 4 ||| d.
Therefore
T Ez 2 S
/3 s 8 |t| 2 (1 + c) min ( - alll + +1|lb
P
With similar argument, can show that
Ex Mislis (1g
(2))
Tzs
s 8. ||t1l3
d + 2. (1 +) (IXb rall + lbbbi²
Thus
Ez
Es (F
[f-
32d
((())
- (0)
10)
2
+
8(1 + t)
min
XXb 14/2 + 100000 (00log(n)(² [[b|13)
n²
n²
We can get high probability bound for the above result using the Hanson- Wright inequality [31, Adjusting the vectors y() and y()) using Xb" as 5(()) y(0)) Xb and 5(1) y(1)) Xb, and
using the Horvitz-Thompson estimator on the adjusted outcomes y:()) and y(1), achieve the following result.
C.1 Unconditional Bound for ATE
The only reason the 1 - å probability appears in our bound for the variance of ATE estimation (Theorem 2) is that leverage score sampling only gives a high probability guarantee for the regression error. Moreover, our variance only depends the regression error, and the estimator is always unbiased, even if the leverage score sampling does not give the (1+e) approximation guarantee on the regression. This "bad event" happens with probability most 6. Therefore, it is enough to bound the regression error when the bad event happens. Here we present such an approach. However, to keep our argument simple, we just work with a single solution vector. To actually apply this to ATE estimation approach, one needs to consider both solution vectors for groups S and 5 and consider the bad event for both. Such an argument follows similar to the following and provides the same result (up to constant factors). This approach essentially can make all of our results unconditional (i.e., remove the 1 - N probability) at the expense of increasing the variance/error bounds by a constant factor and adding
additive term proportional to 8.
X S T SX is not a spectral approximation of T i.e., (1 (1-c)XTX 3 X SX < (1+6)XTX does not hold (this can be checked since we have access to matrix X and only happens with probability at most 8), we just set = In this case, X 21/2 + X bbl = L Otherwise (i.c., (1- x_-e)XXxx X S SX ~ (1 (1+ x)X let X E R dn xd be the matrix X concatenated with factor X of the identity matrix and ça E Rn+dd be the vector fa concatenated with a zero vector and
b" = arg min
Xb
faliz = arg min
Xb- 1al12 + \|bb//
Let X E xd xd be the matrix SX concatenated with the factor À of the identity matrix and g(n+d) be the vector Su concatenated with zero vector and
b = arg min
Xb 2 = arg min SXbb - Sp 1 + A||b|||.
X X (1 c) X). Then since X is spectral approximation of X (i.c., (1 - x) X
XT Í S 18b- Ā+/l2 + || Xb*ll2 < 18b' 2ا112 +
Xb Xb" 1|,, (13)
1-c
24

---

where the first inequality follows by triangle inequality. Moreover, by triangle inequality and optimality of b
Xb Xb" 112 S Xb 2 12 + |Xb' - S 2- |bb 8 - 12122
(14)
Now let Ei be the event that S SX is not spectral approximation of X X. Let & be the event that X S SX is a spectral approximation of _ "X but |SXb - Sp 3 + K A|bb| > (1+ (I 1 [/ + b [3)) Finally let &g be the event that X S SX a spectral approximation of X X, but ||SXb Sui 22 xlbbl S (1+e) X| expectation,
pa 2 + X x||b(l)) Then, by the law of total
- - = E[ - []] P [ + B[( - 7)² إ2ع P [ + E - 7) Ex] IP [E&)
Note that ,[Ei], ,P [E2] < and therefore by the above argument and Theorem
[ - -) 1 5. ||ا| + 5 ( - [2)
(15)
+ (1 6)
[8(1 +e)
min
(IXb 1l12 + 100 log(n/)) C² (b()))
+ n²
32d
-1ly(1) - 1200
n
beRd
Now note that
B[(F - 1)2)
R
8(1 +)
(1 (Xxb-
Xb
1
2
+ 100 log(n/8)) -22
|b]])
[6]
+
32d
11y - y100 E2)
n²
n²
E
8(1 + +c)
(|xx 11) 2 + 1000 (n/6)
.c2 [
2
+
32d
-1y - 11
P
[8]]
n
2
n²
1
E
8(1 + c)
(IIbb' 11/2 100 ") x² ||b 11
+
32d
||y 12
P
[8]]
n
n²
1
E
8(1 + E) 2
(ISXX
Sp + 100 log(
n ")
||b*|
||b
[
n²
1-c
where the last inequality follows from and (1
Now since leverage score sampling preserves
the norms in expectation, we have E Xbb T0ll
Xbb A0l1.. Therefore
B[(f 112) [Ea]
8(1 + c) (1
2 - (ISXD- Sull2 - 100 Th - -|bb*
n²
32d
(1y (0) 11
n
Combining this with 115.. we have
[7 -T) 7)
64d
8(2 22))
$ n²
(1y(1)
12 +
min (IXb 1112 + 100 log(n/0) <2 -1bb11))
+ 8
||111/2
D ATE with Partial Observation
A simple approach for estimating ATE using partial observations is to first uniformly sample a subset of size m of the population and then apply any of the methods for ATE with full observation to this subset. If we observe the outcomes of the samples in the subset independently and with equal probability and apply the Horvitz- Horvitz-Thompson estimator, we achieve an squared error of less
|||l + ||t (1 )" If we use the Gram-Schmidt walk design and the Horvitz- Horvitz-Thompson Thompson than mn
estimator on the uniformly sampled subset, then we achieve an unbiased estimator with a variance of less than
E - [f- <
min BER
1 Xb 1112 + 1
52 |bbl
+
E
1/ (1
m
m
m
n
25

---

Our approach is to adjust the Horvitz-Thane Thompson estimator itself. We call this partial observation regression adjusted Shriiz--hoompoon (RAHT) estimator. The first result we prove in this section is independent of the algorithm we use for estimation over the uniformly selected subset of the population. We introduce the following novel estimator by applying regression adjustment to HT estimator.
Definition 17 (Partial Observation pigressioonAdjuste Horvitz- Thompson (RAHT) Estimator).
Given matrix of covariates X E R"xd, fixed vector b, and an integer 1 I S introduce the partial observation regression adjusted Horvitz- Thompson estimator as the following. Estimation is based on a uniformly sampled set S of size m that partition two sets Z and Z by a random process (e.g., Gram-Schmidi walk). Then the ATE estimate is the following.
Σ LIEZ+ P
Σ iEZ
1
ΣxT iES
(16)
P
The following result indicates that the partial observation regression-adjusted Horvitz-Thompson estimator is a powerful tool for estimating ATE with partial observation.
Lemma 18. Let Le) (0)))(00))))) R", b E R be fixedhreassiined vector and X E Raxd Given
unbiased estimator Ps for Ts = [S] Eies Vi - "² and a uniformly sampled set of size m, the
estimator 7 == Pg ( (x: x) b is unbiased and its variance is
X (
E [f--)) - ESnD E (ffs 11)
(17)
2
where the second expectation on the right-hand side is over the randomness os. and DD denotes the distribution of choosing a uniformly random set of size m.
Proof. We have E
[(F 7) = ESSD
E
(īs
(n m
Eiss(x; xx) - -
Moreover
(fs
1 Σ (x, x) b) - -)² -Ps
Ts + Ts
Σ (x: xx] () 7)
m
m
iES
(fs rs)² 2 + (TS
- Σ (x: - x) b) - 7) + 2( TS)
-x
1 Σ (x: [b) -7).
m
Since ig is an unbiased estimator of Tg.
2(Fs - (T)(TS - ( Σ (x; - x) -x) - -)
1
Eg-D E
m iES
= ESSD
E [2 - rs)] (Ts
Σ
(xi - x) Tb) 7)
= 0.
iES
Moreover
EgnD
E
(TS
1
Σ x) -7)2
ESND
(TS (TS- 1 Σ
(x; - b) -)2
m iES
m $ES
Since Es-D
T8 ( Eiss(x;
- x) [b] = T, by Lemma 9
- m/n
EsnD
(TS - Σ (x, x) b)= - -
X --xb +5x (x(()) - (0)) -
m - 1)
iES
X x)b (t- -t) =1
< mn
26

---

We next prove varaince bound if the Gram-Schmidt walk design is applied for estimation over the uniformly selected subset.
Theorem 19. Let is be estimate obtained by the Horvitz-Thomee Thompson estimator the Gram- Schmidt walk design. Let b be a fixed/preasigned vector and 0 observation regression-adjusted Horvitz-Thompsit estimator with m samples is an unbiased estima- tor with variance of
1. Then, the partial
E [ -)
mn XXb 1el22 m²
1
2²
|bb'1
1
[
x
(t (t-))
(1
where b* arg minDeRd
Xbb 1ul2 +
c²
|
(1-0)
Proof. By Theorem 13
and Lemma 18.
[(f - [ < BS-D
[
min DER"
1
[xs:b- - -ssl2 + 1-
2²
|b|2
m
1
X) 6 (t (t- 8)1{
X
mm
BSwD
1xs:b= - Msll2 + 2²
(x
T
(t
t)
m=
c²
ESND
1xs.b* a-b - psl12
b
m
m 1
1
|X
X) Th 15-tt
-(t - F)|(
mn
Now note that
Xs:b* - -ssl12 = (x, [x.'b= [x, -b,²²
iSS
The probability of any i E [n] appearing in a uniformly sampled set of size m is
n-1
m n
11
Therefore
ESnD
Xg.bb* - -1sl2
= m
(x] b' - 11)
m Xb* n2
-
n
EE(N)
and the result follows by combining the above.
In Lemma 18 and Theorem 19
have assumed that the vector b (which is used for adjusting the PevvitzThhomsonn estimator) is given. For our final result, we combine our approach for regression on random vectors (with partial observation according to leverage score sampling) with partial observation RAHT to obtain a complete algorithm for estimation.
Theorem 3 (ATE Estimation Partial Observation). For any c, E (0, 1), Φ E (0, and m E [n], there exists a randomized algorithm (Algorithm 4 that observes potential outcome for O(dlog(d/8)/e + m) units and outputs an unbiased estimate ē the
the ATE T, for which there is an event & with Pr(E) Z 1 - 8 (over the randomness of the algorithm), such that
E (f - [ <
1
Xbb
1112
1
| 112 + 100d-log(d/d) n²
1 3
mn á
m² (1-0)
+ + E)
m 1
X)6 (t T)l
1
mn
mm
where : maxie(n]
||, 1 1 = y (0) (1)), A = 6-log(d/))
2 E Te is a vector where all the
entries are equal to f
Σ tu and
b' arg min beRd
Xb 11 + 2 ||b|2
6
arg min [[[XX - X)b (t. E)IB +
m
bERd
27

---

Algorithm 4: Regression Adjusted Horvitz-hhompson Estimation
Input: X E guxd y(1, y IR", E EN.0<<<Ll<0<c
Set a = d/n and 6:logdd/4) <²
Compute the X-ridge leverage scores p!7)) of x= [X X) 1
Set p. min maxfa, C (og(d/6) ( For all E [na] define the random variable that is 1 with probability Pt, and is zero,
otherwise. Let S be the set of all indices with qi Let S E Rox" be random diagonal matrix with S 1/ l/VPi if 1, and S 0, otherwise.
(0) Let E R be random vector where ys is independently equal to either yf (2) or Vi with equal probability.
SXb - Srib-syliz + Alb:aall and v E IRd be a vector consisting of
Let V = arg min
the first d entries of V
Let T be uniformly sampled set of size m from [n] \
10 Obtain 3z by Horvitz-Thompsee estimator and Gram Schmiid walk design on T with
parameter é.
// See Theorem 13
11 Let Ts = 2ies
12 Return
S
(fg
-x)
Exe -x
v)
SS
Fs-
Proof. We first prove that the estimator is unbiased. By the law of total expectation, we have
E F r] Σ E -
where the summation is over all sets S that are selected from leverage score sampling. Note that given and y, the solution vector fixed. We denote this fixed solution vector with v Let S E R' be the vector obtained by taking the first d entries of v..y Let x be the average row of Xg.
Then we have
f T|S]
Σ
E
S
Ts
(x:
x5) x
S
Ts S y
PP [y]
71
m
y
S
S
+
Fs
TS S
where the summation over the assignments of the random vector y, x the average row over the of rows indexed by S, 7' uniformly random subset of 5 of size m, Fg is estimated by using GSW design and Horvitz-Thompson estimator on set T, and is is estimated by the Bernoulli design and Morvitz-Thomper is an unbiased estimator of rS, the second term is equal to zero. Moreover since Tg is an unbiased estimator of T5. over the expectation, they cancel out and therefore.
estimator (the observation is compatible with vector y). Trivially since řs
S n
E - Σ
E
(x (x)xx))
|S.y Pyll.
m
Now since &Sy is fixed given S and y, we have
S
E 1
E
Σ
(x;
x5)
FS.. .P PP [x].
mn
LiET
We have
15|-1
E
Σ X4+ x3
x3)
m-1
T
m
Σxi x
0.
(191)
|S|
iET
iES
192
iE3
28

---

Therefore E - T|S] = 0 and B 7] = 0. Thus, the estimator is unbiased.
We define the event & be the event where the regression enror is within factor of (1 e) of error and the number of nonzero entries in the matrix is less than 100d log(d/6) Note that
E
[ T)
E
[ 1) P [[,
3.9
Note that condition and y.
[S|
Σ
[x; - x
T
| + 7
-s.
m
ET
where ig and ns are unbiased estimators for Tg and rs, respectively. Then since T this is unbiased estimator. Now we have
.Tg + [8]
TT.
E [
E
S
15
1
Σ
T
|
Tg
S
Tg
|||
12
m
72
n
5
E
( - x5) T
Tz)
m IE
n
|
S
+ E
TS)
1s),
+
2E
SS]
1
(x+ x3))
x
5
[g)
S
Fs
|S|
Ss
(18)
ET
Note that condition to S. y.
[S]
Fs
|S|
Ts) is fixed number. Therefore
5
S
|S|
|S|
E
T3
[ x3) T
T
-s
SS,
n
m
n
n
ET
S
S
S
[S]
E
[x x3) TySy
75)
S
E
-ss
n
m
n
n
n
Since E [Tgl
Tg and
Lis X; = x the above quantity is equal to zero. Now note that since
[ST
we use the Horvitz- Norvitz-Thompee Thompson estimator to estimate F, we have
-
S
S
|S|
100d log(d/6)
E
TS) 1S,
I 11
1111120
(19)
IS2
ne
Now note that is fixed given S.y. Let
1
2
b
arg min be
|
13l12
|
(1 - Φ)
Then by Theorem
19
E
S
1
X
3 T88y
S
TS)
2
n
m
11
iET
S
2
1
X Xyb5
1
2
2
X
x
x5)
118
t
2
m
[S|
Φ
2
2
|
1
2
Xzh3
1'
X
Sy
mn
m2
(1
3)
m
1 1 mn
1
2
|xbb
|
X 5
v. - (t.
m
1
2
Xbb
-111
lb
2
X -
X
- (
tz²
(20)
mn
m
mn
29

---

We now use the guarantees of leverage score sampling to bound the following
(t" 5 tz
X
The matrix X denotes a matrix with a size equal to X such that all of its rows are equal to the
average row of Similarly, t. is a vector with a size equal to tz such that all of its entries are
equal to the average entry of Eg. Similar to the proof of Theorem 12 let X == -xx-x 1 and
z
X S
1
VAL
ViL
0
Moreover, let 5(1.5), 50,
R +d be vectors that are equal to y . respectively, on their
(1) (0)
first [S] entries, and the rest of their entries are zero. Also let =)(1.3) 50, Then by Theorem the vector V in Algorithm 4 satisfies
E
| 2 - Es
1s
S d
(1 -1+
113 1oc + (1 cc)
min
8s b
bER4+1
2 b 13 2
(1 + c) |111 + (1)
min beRd
Note that
Z b-
t
+ A b
and similar to the proof of Theorem 12. one can see that
x b tsl + A b
min beRd
Xx X
b
(tg
t
b
2
min BER--1
2
Moreover
2 tg
22 V 53
(Xg
X
v
veti 1 1-tg
+
2
Again similar to the proof of Theorem 12 one can show that
X x v (tg
S
[ X J
1 -g
Combining the above and noting that X | 1 2 0. we have
27
X x²) x3x0 xx)0-(ts. (ts -55)
E
d < _ (1+ 5) 1111// (1 +
min
||
||x3 x x51b }b 1bs
- {5)
l
2
+ A
|b12
a
beRd
$ d
(1+c)]1g ||1 +(1+c) X - x²x bb (tg イ5
2
+ > 6
2
Now since (Xs - x 96 x²}6
(ts - (3) is orthogonal to the vector of all ones, we have
(Xg
x 3 16 (tg { / <
X x²16 (tg - {3)
+
X 3
x3) - (t - is)
-is)
(xg Xz)6- (tz -tz)
2
<
X x)b - (t - t)
2
30

---

Table 4: Variance of different methods for estimating ATE with few samples. = y + y(0)
t y y the maximum & norm
is an imput parameter. 0 E is an enror parameter. Ç is
of X. b' arg gid
Xb pll2 + 5- [bb| and
(x - x-b (t- t) ||: + 6 2 where À G-log(d/o) 5
= arg ben
Method
Variance
Uniform Sampling
m St (1
HT Estimator
mn
GS Walk (trailed
min
Xbb 1al2 +
1 [bll + St
(1
uniform sampling)
212 beR
RAHT estimator
8²
||b*|| + 100d log (d/6)
[Xb
-ull2 m² (1 Φ)
|
Regression on
mn Φ
n²²2
random vector
+(1+
|
(X
X)b (t.
6
GS walk design
Therefore
X3 X v x)
-(t t 5
(21)
12
d (1+ ||/ + (1)
||(x (X- 1x--X -X6 (t (t-
-01 + + b
a
the results follows by picking d/n and combining
18
20
and
21
We compared the variance bounds obtained by our approach for ATE estimation with partial observa- tion with other methods in Table 4
31

---

E Detailed Description of Datasets and Extra Experiments
The costliest operation in our algorithms is compute (ridge) leverage scores. This can be performed input-sparsity time (i.e., O(nnz + d²), where nnz denotes the number of nonzero entries in the covariate matrix X [TO). Therefore, our algorithms are much more efficient than to the GSW design, which has running time of ( [88]. This significant difference the running time is also illustrated in our experiments, as shown in Table
Some details of the datasets used in experiments are presented in Table Moreover, the comparison of our ATE approach with two other methods (Lin's 24 and difference in means) is presented in Table and Figure
Table 5: Detailed description of datasets.
Dataset
Boston
IHDP
Twins
Synthetic ITE Syntethic ATE
samples
506
747
32120
10000
10000
features
13
50
50
ATE
0.0
-4.016
0.0064
6.762
56.911 I
100
17.776
2.0
25.972
150.081 I
|4
48.668
9.045
0.303
12.740
71.713
y (1)
24.334
6.464
0.179
10.186
69.579
v/n
(1)
24.334
2.749
0.160
2.883
2.891
v/1
miin
Xb 14/12
9.831
0.846
0.276
3.865
0.200
Vn
min
Xb y
4.915
0.585
0.167
1.940
1.923
Vn.
min
Xb (0)
4.915
0.407
0.149
1.930
1.912
Vn
min
Xb
(y(1)
0.0
0.548
0.155
0.199
3.830
Table Average time (in milliseconds) for estimating ATE over 1000 trials for different methods.
Dataset Uniform HT
GSW
Classic Reg Adj
Lev Score
4 vecs
Ours
Boston
0.1
9.5
0.4
1.1
2.1
1.4
IHDP 0.1
19.8
0.6
1.6
2.5
2.0
Twins
2.8
27719.0
13.7
43.9
54.8
53.2
Table 7: Results of ATE estimation. For each result, the first number is the average of |T 7 over 1000 trials and the second number is the standard deviation of this quantity.
Dasavet
Unifeem HT
GSW
Classic Reg Adi
Lin's AA
diffenence
means
Lev Score
vecs
Ours
1.736
0.663
0.353
0.655
0.658
1.677
0.628
+1.339
#0.510
#0.265
+0.496
+1.256
20.459
IIDP
4.0.208
0.272
0.0031
0.042
+0.009
0.032
$0.010
0.013
1:0036
0.049
4:0.00
0.5d6
4-0.003
U.264
$.0.33
0.040
Twins
+ +1.025e-3
1.351e
+0.937e-3
200 39.8
1.201e-3
41091706.9
56839.0
10.981e-3
1.301e-3
=0.911e-
11.01 5e-
1.3090 3
10.936e-3
32

---

HT uniform
GSW
0.08
classic regression adjustment
leverage based
vectors
0.06
Lin's
ATE eroor
0.04
0.02
0.00
500
750
1000
1250
1500
1750
2000
2250
2500
Population Size
Figure 2: Results for synthetic ATE dataset. For each number of samples, estimation is performed for 1000 trials. Then the average [ r|/|T| over these trials shown with solid lines and the shades around these lines denote the standard deviation.
Gram-Schmidt Walk Design
Harshaw et al 18 introduced Algorithm for balanced experimental design.
Algorithm 5: Non-Uniform Gram Gram-Schmidt Walk
Input: Matrix X E Raad < < 1, vector p E R"
Initialize the vector of fractional assignments (1)) <- (2p1 1,..., 2pn 1).
Initialize an index j f-
Select an initial pivot unit uniformly at random from [n]].
Sct Ç < max (x] |12
Vo 5 Φ X
Let B (n+d) be matrix such that B =
where e; is the i'th
basis vector of dimension n. while z(5) & +11, do
Create the set S (ie [ If & S, select a new pivot k from uniformly at random.
< 1}.
Compute step direction as (f)) arg min, Bull2 u, for all i $ S, us 1}. Let A {6 TR ( 5 n E [_1,
10
Set 8+ < max A and 5 < min A/.
13
Pick random step size ōy that is equal to 8 with probability 8 /(6t + +8), and is equal - with probability 8+ 5 + 8-)
Update fractional assignmenyt: (j+1) (() + أzu(), Increment index j + j + 1.
14
15
16 Return the assigment vector 2 {f- +11.
33


Reviews:
[{"review_id": "vSGk81BGCN", "reviewer": "Reviewer_zpw6", "summary": "The paper addresses the problem of ATE and ITE estimation in the presence of covariates. In particular, the authors provide finite-sample variance bounds for regression adjustment method-based estimators and novel variants thereof. The core of the methodology is using leverage scores, a randomized numerical linear algebra technique. This approach has been previously employed in selecting a sample of the population on which to perform an experiment.\n", "questions": "* How is the random vector in Key Challenge 2 built given that we don't observe both $\\textbf{y}^{(0)}$ and $\\textbf{y}^{(1)}$?\n* It looks like both potential outcomes are observed in the experiments, is that correct? Again, the experimental section is difficult to read. \n* It appears that, in several experiments, the performance of the proposed method is comparable to that of the leverage scores method. In light of this observation, what is the significance or unique contribution of this approach?", "limitations": "I was unable to find a discussion by the authors about the limitations of this work.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "To the best of my knowledge, using leverage scores for regression adjustment and offering finite-sample variance guarantees are novel contributions to the ATE and ITE estimation literature. The paper is strong in its technical aspects, as many of the theoretical contributions are robustly presented and accompanied by thorough proofs.", "weaknesses": "The idea of the paper is overall good, but the execution is poor. The main issue is that the paper very technically dense and is largely unreadable without the appendix. It reads like a stream of (highly technical) consciousness rather than a scientific report. It also looks like the authors sacrificed intuition in favor of technical details of questionable relevance. \n\nOther notes:\n\n* The algorithm links are broken (as they are in the appendix, but that is not mentioned anywhere)\n* The algorithms referenced from the appendix contain quantities not yet defined which makes the paper difficult to read\n* Many details in experimental section are relegated to the appendix so it's difficult to assess its soundness. \n\nIn my opinion, while the paper has potential to be a valuable addition to the conference, it would require substantial rewriting in order to meet the required standards and expectations.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for positive assessment of our technical contributions and for recognizing the novelty of applying leverage scores to regression adjustment. We agree with the reviewer that the paper can be made more readable and appreciate the reviewer’s suggestions towards doing so, which we will address in the revised version.\n\n\n\n> How is the random vector in Key Challenge 2 built given that we don't observe both $y^{(0)}$ and $y^{(1)}$?\n\nTo construct the random vector v, we only need to observe one of the outcomes $y^{(0)}$ or $y^{(1)}$ for each unit. In particular, for each unit $i$, we toss a fair coin. If it lands on heads, we observe $y^{(0)}_i$ and set $v_i=-2y^{(0)}_i$. If it lands on tails, we observe $y^{(1)}_i$ and set $v_i=2y^{(1)}_i$.\n\n> It looks like both potential outcomes are observed in the experiments, is that correct?\n\nNo, for each unit, we observe at most one of the potential outcomes. The only method that uses both observations is the baseline for ITE which is used for comparison.\n\n> It appears that, in several experiments, the performance of the proposed method is comparable to that of the leverage scores method. In light of this observation, what is the significance or unique contribution of this approach?\n\n\nOur goal was to develop a treatment effect estimation approach that is fast and unbiased, works on finite-population, gives non-asymptotic guarantees, allows only the use of a subset of the population (instead of the whole population), and can be used in different settings (such as online experimental design).\n\nAs we mention in the introduction, each previous method has some shortcomings. So the goal of our method is to address all of these shortcomings simultaneously. For example, the GSW-design is very slow, as our experiments (presented in the appendix) suggest — on the twins dataset with about 32 thousand samples, the GSW design is about 500 times slower than our algorithm. Another example is that classic regression adjustment and using leverage scores to learn two different vectors are biased approaches. \n"}, {"review_id": "IDm9iUFrXz", "reviewer": "Reviewer_HGb9", "summary": "This paper explores the design and analysis of randomized experiments for treatment effect estimation, which is an important problem in causal inference. The goal of treatment effect estimation is to estimate the effect of a specific treatment on individual subjects or the average effect in the population, using only one of the two potential outcomes (either treatment or control) per subject. The objective is to obtain a precise estimator of the treatment effect, which is unbiased and has a smaller variance. However, it is often also desirable to minimize the number of subjects exposed to the experiment due to practical or ethical considerations. To address this, recent research has focused on leveraging linear algebraic techniques to design optimal experiments and estimate treatment effects.\n\nThe authors identify three main limitations in previous approaches: (1) inefficiency resulting from using the entire population in the experiment, (2) sub-optimal variance bound of the resulting estimator, and (3) inflexibility in experimental design, particularly in online settings. To overcome these challenges, the authors propose combining subsampling techniques from numerical linear algebra with classical statistical methods of regression adjustment. This integrated approach aims to obtain an unbiased estimator with a variance comparable to the optimal value, while also being suitable for online experiment settings.\n\nThe paper is organized as follows. In Section 2, the authors demonstrate the effectiveness of their proposed approach by considering four types of problems: (1) mean estimation of a vector, (2) individual treatment effect (ITE) estimation, (3) average treatment effect estimation (ATE) with full observation of both potential outcomes, and (4) ATE with partial observation, where only one potential outcome per subject is available. The primary focus is on problems (2) and (4), while problems (1) and (3) serve as precursors to facilitate the exposition of ideas. In Section 3, the authors describe the key techniques employed in their approach, providing a brief outline of the proofs for some theorems. Subsequently, they present experimental results from synthetic datasets (Section 4.1) and real-world datasets (Section 4.2) to validate the effectiveness of their proposed method.", "questions": "1.  I would like to request the authors to address the the concerns raised in the “Weaknesses” section. \n\n2.  It appears that the error bound in Theorem 2 may be trivial and not particularly informative because $ d \\| y^{(1)} + y^{(0)} \\|_{\\infty}^{2} \\geq \\| y^{(1)} + y^{(0)}\\|_{2}^{2} \\geq \\| y^{(1)} - y^{(0)}\\|_{2}^{2} = \\| t \\|_{2}^{2} $ assuming $\\|  \\| y^{(1)} - y^{(0)}\\|_{2} \\| \\leq \\|  \\| y^{(1)} + y^{(0)}\\|_{2} \\| $. I wonder if the authors clarify the value of Theorem 2? t would be beneficial for the authors to clarify the significance and value of Theorem 2, particularly in relation to its practical implications. \n\n3.  In Section 4.1, it would be helpful to provide information about the population size (n) for each dataset used in the experiments. Additionally, it would be valuable for the authors to evaluate the methods by quantifying the required sample size (i.e., the size of the sub-population) to achieve comparable performance to the method based on the entire population. This analysis would provide insights into the practical efficiency of the proposed method and its advantages in terms of minimizing the required sample size. \n\n4.  In Section 4.2, it appears that the classic regression adjustment method outperforms all the other methods, including the authors' proposed method. In light of these results, it would be important for the authors to provide a clear justification for the value and effectiveness of their proposed method, especially in comparison to the superior performance of the classic regression adjustment method.  \n\n5.  In line 267, the authors mention that the classic regression adjustment technique is biased. While this is true, it is worth noting that there are at least two standard unbiased techniques available, namely the difference-in-means method and Lin's interacted regression adjustment [20]. It would be valuable for the authors to include these unbiased techniques in their experiments to provide a comprehensive comparison and analysis of the proposed method against these established approaches.\n\nAddressing these concerns would greatly enhance the understanding, validity, and practical applicability of the research presented in the paper.", "limitations": "It would be valuable for the authors to provide a short discussion on the technical limitations of their methods and offer insights into the potential adverse consequences that may arise when applying these approaches in real-world settings, specifically in experiment design and treatment effect estimation. Nonetheless, it is important to note that given the paper's primary focus on theoretical aspects, a comprehensive and extensive examination of these limitations may not be deemed critical.", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "The paper demonstrates several strengths that contribute to its overall quality and significance. Firstly, the problem of treatment effect estimation addressed in the paper holds substantial importance within the field of causal inference. The research questions explored are highly relevant and carry practical implications, underscoring the significance of the study.\n\nSecondly, the proposed approach of integrating subsampling techniques rooted in randomized numerical linear algebra with regression adjustment is a well-founded and suitable strategy for tackling the problem at hand. Additionally, the authors' systematic approach in Section 2 effectively illuminates the challenges associated with addressing the problem and provides a rationale for the proposed approach.\n\nFurthermore, the technical tools employed in the paper are well-established and widely utilized in the field. The authors' adept development of methods and analytical arguments reflects a strong grasp of these tools, bolstering the credibility of their findings.\n\nIn summary, the strengths of the paper lie in the significance of the addressed problem, the logical and appropriate nature of the proposed approach, and the authors' proficient utilization of standard technical tools. These strengths collectively enhance the quality and validity of the research presented in the paper.", "weaknesses": "While  the paper demonstrates strengths, there are several areas where improvements can be made to enhance clarity and strengthen the overall presentation. The weaknesses can be broadly categorized into two main areas: clarity and organization, as well as the soundness of the arguments. \n\n\nI.  Concerns regarding clarity and organization\n\nFirstly, the paper lacks clarity, potentially resulting from a lack of organization. To improve clarity, the authors should clearly identify the specific problems they are focusing on and indicate which parts of the paper address these issues. This could be achieved by adding a dedicated subsection or table in the Introduction that summarizes the problems addressed, the solutions proposed, and references to specific sections and theorems.\n\nAdditionally, restructuring sections would contribute to better organization and understanding. Merging Sections 2.1 and 2.2 into a single subsection titled \"ITE Estimation\" and merging Sections 2.3 and 2.4 under the title \"ATE Estimation\" would streamline the presentation and emphasize the two treatment effect estimation problems. Furthermore, Section 3 would benefit from reorganization and clearer subsection headings to provide explicit explanations of the specific techniques being described. For example, Section 3.1, currently titled \"Random vector regression,\" contains a combination of an auxiliary theorem (Theorem 8) used to prove Theorem 2 (specifically, the display equation in line 113) related to ITE estimation, as well as a proof sketch of Theorem 3 in Section 2.3 related to ATE estimation. The rationale for grouping these two contents together is unclear, and it would be beneficial to clarify this and consider separating them into distinct subsections.\n\nMoreover, the paper contains numerous theorem statements that, although likely mathematically correct, pose challenges in interpretation for several reasons. Many of these theorems reference algorithms that are either presented in later parts of the paper or not included in the main text at all. This creates logical inconsistencies and forces the reader to navigate back and forth, hindering comprehension. Furthermore, the authors do not provide explicit remarks on the optimal variance or the origin of the terms used in the error bounds, particularly, neither in close proximity to the theorem statements nor in a designated discussion section. Consequently, it is difficult to grasp the significance of the error bounds in terms of their semantic interpretation and their implications, including how tight or loose they are compared to optimal errors.\n\n\nII. Concerns regarding soundness\n\nMoreover, there are concerns regarding the soundness of the work, which may be further exacerbated by the clarity issues. The authors claim that their proposed method successfully addresses the aforementioned three challenges of previous approaches, but this claim lacks adequate support. To strengthen their arguments, the authors should provide explicit discussions or remarks, either following relevant technical propositions or in a dedicated discussion section. Specifically, they should elaborate on the following aspects in clear and comprehensible language:\n\n1.  Clarify why previous approaches, such as [16], require the entire population in the experiment, while the proposed method does not.\n\n2.  Discuss the optimal values or lower bounds of variance and how the authors' error bounds can be considered \"comparable to the optimal.\"\n\n3.  Explain how the proposed method allows for adaptation to the online experiment design setting.\n\nAdditionally, it would be beneficial to appropriately compare the proposed method to baseline approaches, such as classical regression adjustment estimators, Lin's interacted regression adjustment [20], or the crude difference-in-means estimator. Conducting such comparisons, particularly when the sample size is equal to the entire population, would provide a clearer understanding of the strengths and limitations (if any) of the proposed approach.\n\n\nAddressing these weaknesses would significantly enhance the clarity, soundness, and overall quality of the paper.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their in depth review and valuable suggestions, especially those regarding increasing clarity of the writing and paper organization, which we will incorporate in a revised version.\n\n> To improve clarity, the authors should clearly identify the specific problems they are focusing on and indicate which parts of the paper address these issues.\n\nWe have stated the specific problems in Section 1.2. We consider 4 problems: 1) mean estimation; 2) ITE; 3) ATE with full observation (i.e., for each unit, exactly one of the outcomes are observed); 4) ATE with partial observation (i.e., for each unit, at most one of the outcomes are observed).\n\n> Furthermore, the authors do not provide explicit remarks on the optimal variance or the origin of the terms used in the error bounds…\n\nWe will add a discussion about terms in error bounds. One general way of looking at our bounds is that regression adjustment replaces the population variance with the error of best linear fit.\n\nWe note that the Bernoulli design is min-max optimal (see Lemma 4.1 of [16]).\n\n> Clarify why previous approaches, such as [16], require the entire population in the experiment, while the proposed method does not.\n\nFor [16], this is essentially because of the random walk approach and is inherent to the algorithm. Generally, for previous methods in the literature, including [16], one can uniformly subsample the population and then design an experiment on the subsample and make an estimation, but this leads to worse bounds compared to ours — see the discussion in Section 2.4. Essentially the main reason that our subsampling approach works is the guarantees that leverage score sampling provides for minimizing linear regression error, and the fact that we show the variance of the population can be replaced with the error of the best linear fit.\n\n> Explain how the proposed method allows for adaptation to the online experiment design setting.\n\nNote that in leverage score sampling, we only require upper bounds on leverage scores (not the exact leverage score), as stated in Theorem 7. Then by the constants we have picked in our algorithms/theorems we make sure that the probability of picking a unit for outcome 0 or outcome 1 is less than 0.5. Then essentially, the algorithm only needs to toss a fair coin for each unit arriving in the online setting. The only requirement here is to have the largest singular value of the covariate matrix to be bounded during the online experiment (because of our constant). We note that this is a very mild assumption. We will add a discussion about this. A similar argument applies to our subsampling (partial observation) approach as well.\n\n> It appears that the error bound in Theorem 2 may be trivial and not particularly informative because $d||y^{(1)}+y^{(0)}||\\_{\\infty}^2\\geq||y^{(1)}+y^{(0)}||\\_2^2\\geq||y^{(1)}-y^{(0)}||\\_2^2=||t||\\_2^2$ assuming $||y^{(1)}-y^{(0)}||\\_2^2\\leq||y^{(1)}+y^{(0)}||\\_2^2 $. I wonder if the authors clarify the value of Theorem 2?\n\nNote that $d$ is the number of covariates/features. $y^{(1)}$ and $y^{(0)}$ are n-dimensional vectors, where n is the population size. Therefore the first inequality you mentioned does not hold. Typically, we have $d \\ll n$. Thus, $d||y^{(1)}+y^{(0)}||\\_{\\infty}^2$ is much smaller than $||y^{(1)}+y^{(0)}||\\_2^2$. As an example, pick $\\alpha=\\frac{d}{\\sqrt{n}}$. In this case, the number of observed units will be $O(d\\log(d/\\delta)/\\epsilon^2+d\\sqrt{n})$ and the mentioned term in the variance bound is $\\sqrt{n}||y^{(1)}+y^{(0)}||\\_{\\infty}^2$ which could be much smaller than $||y^{(1)}+y^{(0)}||\\_{2}^2$ when n is large.\n\nIn addition, note that $y^{(1)}+y^{(0)}$ is indeed the vector $\\mu$ which is the vector appearing in Horvitz-Thompson estimator and Gram-Schmidt walk.\n\nIn practice, the size of the population is usually much larger than the number of features/covariates, and we believe in such settings, our bound is very useful.\n\n> In Section 4.2, it appears that the classic regression adjustment method outperforms all the other methods, including the authors' proposed method...\n\nWe would like to note that as opposed to classic regression adjustment, we present non-asymptotic bounds for our approach. To our knowledge, the guarantees on classic regression adjustment are asymptotic. Moreover, as you mentioned, classic regression adjustment is biased. Finally, note that it performs very poorly on our synthetic dataset. In general, we believe that no single approach performs well on all datasets and all settings (similar to the no-free-lunch theorem in machine learning). Our approach is most effective in scenarios where the regression errors are small, i.e., the treatment effect is linearly predicted by covariates well.\n\n> In line 267, the authors mention that the classic regression adjustment technique is biased. While this is true, it is worth noting that there are at least two standard unbiased techniques available, namely the difference-in-means method and Lin's interacted regression adjustment [20]...\n\nThank you for pointing this out. We have added experiments for comparing the performance of the unbiased versions of classic regression adjustment to the pdf accompanying our response. Interestingly, Lin’s approach works very well on all except one of the datasets. However, on the twins dataset, it performs very poorly (but this might be due to numerical stability issues). Note that Lin’s approach is only asymptotically unbiased, while the focus of our work is to provide Non-asymptotic Guarantees and unbiasedness.\n"}, {"review_id": "o5jJtu4n7V", "reviewer": "Reviewer_ZgSr", "summary": "This paper focuses on estimation of individual and average treatment effects (ITE and ATE) with regression adjustment, which is combined with the method of ridge leverage score sampling in order to obtain the desirable variance bounds. The leading case is algorithm 1, which estimates ATE with leverage score sampling and cross adjustment. In this algorithm, for each observational unit, either treated or control unit is observed with equal probability. Then observed outcomes are adjusted via cross regression adjustments for which the regression coefficients are computed using ridge leverage score samples. The paper provides a number of theoretical results regarding the variance bounds and reports experimental results using both synthetic and real-data datasets.   ", "questions": "- Is there any chance to obtain lower bounds for variance? It would complement the existing results on the upper bounds.\n- Would it be possible to estimate the variance bounds?\n- Line 192, is the transpose missing for the last $X$?\n- Steps 5 and 6 in algorithm 1 may contain typos. I think these steps are meant for leverage score sampling but the stated steps are unclear. Please check them.\n- Line 264. It might be helpful to provide a summary of running times in the main text.\n- There is very little discussion on the experimental results in section 4. Probably this is due to the page limit. It would be helpful if there is a short summary of numerical findings in the main text.", "limitations": "One important limitation that is not mentioned in the paper is how to conduct statistical inference on ATE. This is typically done using asymptotic arguments. For example, Harshaw et al [16] provides asymptotic normality and suggests a confidence interval for ATE based on asymptotic normality. Although this paper focuses on finite sample variance bounds, it might be fair to mention the lack of inference methods. ", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The main strength of the paper is that it provides finite sample variance bounds for estimating the sample mean, individual treatment effects, and the average treatment effect with regression adjustment.\n- The research question addressed in the paper is of general nature and is important in many applications.", "weaknesses": "- To my reading of the paper, the most important issue is that it is unclear whether the advances in this paper are substantial relative to Harshaw et al [16]. It seems that theoretically, the variance bounds are comparable between the proposed methods in this paper and the Gram-Schmidt walk (GSW) design method of Harshaw et al [16]. Furthermore, the right panel of Figure 1 shows that GSW outperforms the current paper.\n- The paper criticizes the GSW design, saying that it is not suitable for online experimental design (line 33 and lines 139-142). However, the proposed method is not immediately suitable for online settings. For example, it is unclear how to implement ridge leverage score sampling in a fully online fashion because the leverage score changes as more observations are available.  \n- The one advantage of the proposed method over GSW might be computational speed in view of table 6 in the supplement. However, there is no theoretical analysis in terms of computational aspects. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer very much for their review and their questions/suggestions, which we address below.\n\n> To my reading of the paper, the most important issue is that it is unclear whether the advances in this paper are substantial relative to Harshaw et al [16]. It seems that theoretically, the variance bounds are comparable between the proposed methods in this paper and the Gram-Schmidt walk (GSW) design method of Harshaw et al [16]. Furthermore, the right panel of Figure 1 shows that GSW outperforms the current paper.\n\nNote that although the performance of GSW is similar to our approach, it is much slower than our algorithm. For example, on the twins dataset with about 32 thousand samples, the GSW design is about 500 times slower than our algorithm. \n\n> The paper criticizes the GSW design, saying that it is not suitable for online experimental design (line 33 and lines 139-142). However, the proposed method is not immediately suitable for online settings. For example, it is unclear how to implement ridge leverage score sampling in a fully online fashion because the leverage score changes as more observations are available.\n\nWe would like to emphasize that we did not intend to criticize the GSW design. It certainly is an important work in the intersection of statistics and linear algebra. Our goal was only to address the differences between our approach and GSW design. We will make sure this is reflected in our writing.\n\nRegarding the online setting, we like to note that since we are using ridge leverage score sampling, as we pick the regularization (ridge) parameter, we make sure that the probability of selection of each row for any of the potential outcomes is at most 0.5. More specifically, note that in leverage score sampling, we only require to use upper bounds on leverage scores (not the exact leverage score), as stated in Theorem 7. Then by the constants we have picked in our algorithms/theorems we make sure that the probability of picking a unit for outcome 0 or outcome 1 is less than 0.5. Then essentially, the algorithm only needs to toss a fair coin for each unit arriving in the online setting. The only requirement here is to have the largest singular value of the covariate matrix to be bounded during the online experiment (because of our constant). We note that this is a very mild assumption. We will add a discussion about this. A similar argument applies to our subsampling (partial observation) approach as well. In addition, note that leverage scores can be computed in an online manner. See the following paper.\nCohen, Michael B., Cameron Musco, and Jakub Pachocki. \"Online row sampling.\" Theory of Computing 16, no. 1 (2020): 1-25.\n\n> The one advantage of the proposed method over GSW might be computational speed in view of table 6 in the supplement. However, there is no theoretical analysis in terms of computational aspects.\n\nThanks – we will add a discussion about the running time of our algorithm. Roughly, leverage scores can be computed in input-sparsity time. More specifically in $O(\\text{nnz}(X) + d^{\\omega})$ where $\\text{nnz}(X)$ is the number of nonzero entries of matrix $X$ and $\\omega$ is the matrix multiplication exponent. In the worst case, when $X$ is a dense matrix, this is $O(nd + d^{\\omega})$. Note that this is the bulk of our computation since the running time of other steps of our algorithms is also $O(nd)$. In comparison, the running time of GSW design is $O(n^2 d)$. Note that since typically $n \\gg d$, our running time is much better than GSW.\n\nIn addition, since our operations involve matrix-matrix or matrix-vector multiplications, this is naturally parallelized in computers. In comparison, GSW design requires a random walk with n steps, which makes it less parallelizable. Therefore in practice, one might see even more dramatic differences in terms of running time.\n\n> Is there any chance to obtain lower bounds for variance? It would complement the existing results on the upper bounds.\n\nWe will add a discussion about terms in the error bounds, but one general way of looking at our bounds is that regression adjustment replaces the population variance (which is equivalent to using a vector of all zero for adjustment) with the error of best linear fit.\n\nWe note that the Bernoulli design is min-max optimal (see Lemma 4.1 of [16]), and the following works provide hardness bounds on the discrepancy part of the GSW design. We will add a discussion regarding these to the paper.\nZhang, Peng. \"Hardness Results for Minimizing the Covariance of Randomly Signed Sum of Vectors.\" arXiv preprint arXiv:2211.14658 (2022).\nCharikar, Moses, Alantha Newman, and Aleksandar Nikolov. \"Tight hardness results for minimizing discrepancy.\" In Proceedings of the twenty-second annual ACM-SIAM symposium on Discrete Algorithms, pp. 1607-1614. Society for Industrial and Applied Mathematics, 2011.\n\n>Would it be possible to estimate the variance bounds?\n\nSince the variance bounds depend on the regression error of vectors to which we do not have access (for example, $y^{(1)}+y^{(0)}$ since for each $i$, we only can observe either $y^{(1)}_i$ or $y^{(0)}_i$.), it seems non-obvious how to obtain meaningful variance estimates. However, this is an interesting suggestion to study for future work.\n"}, {"review_id": "QkiNB5Fvc9", "reviewer": "Reviewer_UdXn", "summary": "In this paper, authors present regression adjusted estimators for estimating the average treatment effect under the Bernoulli design.\nIn particular, they show that by using the leverage scores and a ridge regression adjustment, favorable finite sample bounds on the (conditional) variance may be obtained.\nThis contributes nicely to a long-standing debate about the use of regression adjustment in randomized control trials.\nFinally, simulations are performed which coroborate theoretical results.", "questions": "- Can you address my comment on conditional versus unconditional analysis?\n- Can you address my comment on the constant $C$ that is used?\n- How does the Regression Adjusted Horvitz---Thompson estimator relate to double robust estimators?", "limitations": "yes", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The main strength of this paper is a finite sample analysis of regression adjusted estimators under the Bernoulli design.\nThere has been much debate about whether regression adjusted estimators are appropriate for estimates of the average treatment effect, [12, 20].\nHowever, that entire literature has been studied under asymptotic analyses, which provides little insight into finite sample performance.\n\n(Harshaw et al, 2019) have shown that under a sophisticated experimental design, favorable finite sample bounds on the variance may be obtained.\nMoreover, in large samples, (using a well-chosen parameter), the Gram--Schmidt Walk design achieves the same asymptotic variance as Lin's regression adjustment. \nThis puts regression adjustment on the defensive and raises a natural question: can regression adjusted estimators achieve similar finite sample results?\n\nThis paper does a really remarkable job of answering this question in the affirmative, up to some minor details.\nIt will be of interest not only to the NeurIPS community, but also to the statistical causal inference community more broadly.", "weaknesses": "The paper has a few minor issues, but I believe they can all be addressed by the authors in a satisfactory way without too much work.\n\n ## Weakness 1 -- Conditional Analysis\n  \nThe major weakness is a certain ambiguity that arises in a few of the theorems.\nFor example, Theorem 3 does not analyze the estimator in Algorithm 1 unconditionally.\nRather, it analyses it conditioned on a high probability event (i.e. \"with probability $1 - \\delta$, Algorithm 1 computes an unbiased estimator of the average treatment effect with variance...\")\nAlthough this is natural in computer science, this is very unnatural in this area of causal inference.\nThe reason is that we place larger value on the unconditional statistical inference; the validity / width of confidence intervals depends on the bias and variance, not a conditional bias or variance.\nMoreover, Theorem 3 suggests that in order to get small variance, we should be setting $\\delta$ very large; however, this will incur a large bias, which is not discussed in the Theorem nor in the paper.\n\nIn order for the causal inference community to fully appreciate these results (and they will!) the authors should derive bounds on the unconditional bias / variance.\nThis will have the additional benefit of giving experimenters advice on how to choose $\\delta$ to minimize the MSE.\nThis critique applies to Theorem 4 as well.\n\n## Weakness 2 -- Constant $C$\n\nAs it stands, the constant in Theorem 7 is presently unspecified.\nThis means that the Algorithm 1 cannot be run!\nAuthors should give a precise constant here.\n\n## Weakness 3 -- ITE Results are Weak\n\nTheorem 2 shows that the collection of ITEs can be estimated up to a constant error.\nThis is unsurprising because ITE estimation is considered impossible due to the fundamental problem of causal inference.\nI find this estimator and these results not very exciting because it's not clear what a practitioner would do with such a guarantee -- not a single ITE is really known.\nIf authors are pressed on space, I recommend moving this to the appendix; otherwise, it's fine to leave in the paper.\n\n## Weakness 4 -- Comparisons to GSW\n\nI think that the comparison to (Harshaw et al 2019) is great and highlights the relevance of this work.\nHowever, some of the comparisons are not relevant or could be improved.\nHere is a list of improvements:\n\n- (Line 137): Authors write that the variance of GSW is proportional to \"XYZ\". However, that's missing the important aspect of (Harshaw et al 2019), which is that these terms can be traded-off! So, I think you have to address (if only very briefly) that GSW-Design allows for trade-off of these terms, but that for $\\phi = 1/2$ it is proportional to \"XYZ\".\n- (Lines 33, 39, 141, etc): Authors write that GSW Design does not have an online analogue (which is currently true), but I don't think that this is how I would compare your estimator. The practical benefit of your regression adjustment estimator is that the data collection process is just so much simpler -- the Bernoulli design is decentralized, asynchronous, you name it! So although (Harshaw et al 2019) raised an online design as an open problem, it's really only interesting for covariate balancing. Being able to use such a simple design is beneficial much beyond online assignment. I think focusing on the practical simplicity of the Bernoulli design over the practical complexities of the GSW Design is much stronger argument in favor of your estimator, not the fact that Bernoulli design can be implemented online.\n- (Line 264): Authors write that GSW Design \"is a much slower algorithm\". This doesn't seem true. (Harshaw et al 2019) give an implementation of GSW-Design that runs in $O(n^2 d)$, which would seem to also be required by for the linear system solve required to compute the ridge regression and scores.\n\n## Minor Misc Comments\n\nThe following are minor miscellaneous comments\n\n- (Line 126): \"work\" seems like the wrong word, a typo.\n- (Line 242): The notation $\\mathbf{P}_{\\mathcal{P}}$ is used but not defined. I think it's fine to just remove the subscript.\n- (Line 203): I think it is confusing / inappropriate to call Algorithm 1 \"leverage score sampling\" because that seems to suggest that the assignment to treatment and control (i.e. \"sampling\") is based on leverage scores. Can you call this something different? Perhaps \"cross leverage score adjustment\" or something?\n- (Section 4.2): authors don't specify the value of $\\phi$ used for the Gram--Schmidt Walk Design.\n- It's a shame that most of the algorithms are in the appendix. If possible, I would recommend bringing them to the main body for camera ready version, perhaps at the cost of putting ITE estimation results in the appendix.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive assessment of our work and our contributions to finite sample analysis of regression adjustment. We thank them also for their detailed suggestions, especially those on the presentation, which we will address in a revised version of the paper.\n\n> The major weakness is a certain ambiguity that arises in a few of the theorems. For example, Theorem 3 does not analyze the estimator in Algorithm 1 unconditionally...\n\nThank you for pointing this out. We also believe that stating unconditional bounds would improve the paper. We will add the following discussion to the paper.\n\nThe only reason the $1-\\delta$ probability appears is that leverage score sampling only gives a high probability guarantee for the regression error. Moreover our variance only depends on the regression error and the estimator is always unbiased, even if the leverage score sampling does not give the $(1+\\epsilon)$ approximation guarantee on the regression. This “bad event” happens with probability at most $\\delta$. Therefore it is enough to bound the regression error when the bad event happens. \n\nIf A^T S^T S A is not a spectral approximation of A^T A (this can be checked since we have access to matrix A and only happens with probability at most $\\delta$), we just set $b=0$ and in this case, $||Xb-\\mu||\\_2^2 + \\lambda ||b||\\_2^2=||\\mu||\\_2^2$.\nOtherwise let $\\hat{X}$ be the matrix $X$ concatenated with identity and $\\hat{\\mu}$ be the vector $\\mu$ concatenated with a zero vector and $b^*:=\\arg\\min||\\hat{X}b-\\hat{\\mu}||\\_2^2=\\arg\\min ||Xb-\\mu||\\_2^2+\\lambda||b||\\_2^2$. Let $\\tilde{X}$ be the matrix $SX$ concatenated with identity and $\\tilde{\\mu}$ be the vector $S\\mu$ concatenated with a zero vector and $\\tilde{b}:=\\arg\\min || \\tilde{X} b - \\tilde{\\mu} ||\\_2^2$. Then since $\\tilde{X}$ is a spectral approximation of $\\hat{X}$,\n$$ || \\hat{X} \\tilde{b} - \\hat{\\mu} ||\\_2 \\leq || \\hat{X} b^* - \\hat{\\mu} ||\\_2 + || \\hat{X} \\tilde{b} - \\hat{X} b^* ||\\_2 \\leq || \\hat{X} b^* - \\hat{\\mu} ||\\_2 + \\frac{1}{1-\\epsilon}|| \\tilde{X} \\tilde{b} - \\tilde{X} b^* ||\\_2 $$.\nMoreover\n$$||\\tilde{X}\\tilde{b}-\\tilde{X}b^*||\\_2\\leq ||\\tilde{X}\\tilde{b}-\\tilde{\\mu}||\\_2+||\\tilde{X} b^*-\\tilde{\\mu}||\\_2\\leq 2\\cdot||\\tilde{X}b^*-\\tilde{\\mu}||\\_2$$\nNow since leverage score sampling preserves the norms in expectation, we have $\\mathbb{E}\\_S||\\tilde{X}b^*-\\tilde{\\mu}||\\_2=||\\hat{X}b^*-\\hat{\\mu}||\\_2$.\nThen by the law of total expectation,\n$$\\mathbb{E}[(\\widehat{\\tau} - \\tau)^2] \\leq (1-\\delta) \\left( \\frac{32d}{n^2} \\cdot || y^{(1)} - y^{(0)}||\\_{\\infty}^2 + \\frac{8 (1+\\frac{2}{1-\\epsilon})}{n^2} \\min\\_{b} \\left(||X b - \\mu||\\_2^2 + 100 \\log(n/\\delta) \\cdot \\zeta^2  \\cdot ||b||\\_2^2 \\right) \\right) + \\delta \\cdot || \\mu ||\\_2^2.$$\n\n> As it stands, the constant in Theorem 7 is presently unspecified. This means that the Algorithm 1 cannot be run! Authors should give a precise constant here.\n\nThe analysis from other papers (see [9]) shows that it is enough to pick $c$ about 10. We add a note about this to the paper.\n\n> Theorem 2 shows that the collection of ITEs can be estimated up to a constant error. This is unsurprising because ITE estimation is considered impossible due to the fundamental problem of causal inference. I find this estimator and these results not very exciting because it's not clear what a practitioner would do with such a guarantee -- not a single ITE is really known. \n\nAlthough we agree that single ITEs cannot be estimated using our approach (and because of the fundamental problem of causal inference), we note that such bounds can be used to estimate other summary statistics. For example, a variance bound for ATE can be derived from our ITE bounds using the Cauchy-Schwarz inequality. In addition, bounds on median ITE or quantiles can be derived from our approach. Exploring these directions more is an interesting avenue for future work.\n\n> Authors write that the variance of GSW is proportional to \"XYZ\". However, that's missing the important aspect of (Harshaw et al 2019), which is that these terms can be traded-off! So, I think you have to address (if only very briefly) that GSW-Design allows for trade-off of these terms, but that for $\\phi=1/2$ it is proportional to \"XYZ\".\n\nThanks for this suggestion. We will add a note about the trade-off of GSW design to the paper.\n\n> Authors write that GSW Design does not have an online analogue (which is currently true), but I don't think that this is how I would compare your estimator. The practical benefit of your regression adjustment estimator is that the data collection process is just so much simpler -- the Bernoulli design is decentralized, asynchronous, you name it! ...\n\nThank you for recognizing the generality and applicability of our results. We appreciate the suggestions and will add this discussion to the paper.\n\n> Authors write that GSW Design \"is a much slower algorithm\"...\n\nThanks – we will add a discussion about the running time of our algorithm. Roughly, leverage scores can be computed in input-sparsity time. More specifically in $O(\\text{nnz}(X) + d^{\\omega})$ where $\\text{nnz}(X)$ is the number of nonzero entries of matrix $X$ and $\\omega$ is the matrix multiplication exponent. In the worst case, when $X$ is a dense matrix, this is $O(nd + d^{\\omega})$. Note that this is the bulk of our computation since the running time of other steps of our algorithms is also $O(nd)$. In comparison, the running time of GSW design is $O(n^2 d)$. Note that since typically $n \\gg d$, our running time is much better than GSW.\n\nIn addition, since our operations involve matrix-matrix or matrix-vector multiplications, this is naturally parallelized in computers. In comparison, GSW design requires a random walk with n steps, which makes it less parallelizable. Therefore in practice, one might see even more dramatic differences in terms of running time.\n"}]

Author rebuttal:
We would like to thank all reviewers for the positive assessment of our contributions to finite-population treatment effect estimation and its non-asymptotic analysis. We appreciate the comments on the presentation and organization of the paper. We believe these would make the paper stronger, and we will address these comments in the next version of the paper.

Below we address the comments and questions of reviewers individually. The attached PDF contains some extra experiments requested by the reviewers. More specifically, it contains the comparison with the difference-in-means method and Lin's interacted regression adjustment approach.

Meta-review: This paper presents the first finite sample analysis for regression-adjusted estimators under the Bernoulli design in the finite population causal inference setting. The reviewers unanimously agree that these contributions are significant and substantive. The quality of the feedback has been exceptional, I won't repeat them in detail in my meta-review and instead highlight the most important concern raised by the reviewing team: presentation.

Even with the challenges of writing for an interdisciplinary community (causal inference, TCS, ML), the current submission falls short of providing a clear exposition by most measures. The reviewing team represented experienced and mature researchers with diverse backgrounds, and it is concerning that all of them initially found the manuscript inscrutable. They have provided thoughtful feedback on how to improve the presentation, and I strongly urge the authors to carefully follow them to maximize the impact of their work.
Decision: Accept (poster)