Robust Signal Maximization in Spillover Experiments

Robust Signal Maximization in Spillover Experiments

\global\long\def\expec#1#2{\mathbb{E}_{#1}\left[#2\right]} \global\long\def\Pr#1#2{\mathrm{Pr}_{#1}\left[#2\right]} \global\long\def\var#1#2{\mathrm{Var}_{#1}\left[#2\right]} \global\long\def\cov#1#2{\mathrm{Cov}_{#1}\left[#2\right]} \global\long\def\corr#1#2{\mathrm{Corr}_{#1}\left[#2\right]} \global\long\def\V#1#2{\mathcal{V}_{#1}\left[#2\right]} \global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\op#1{\left\Vert #1\right\Vert _{\infty}} \global\long\global\long\global\long\def\Frob#1{\left\Vert #1\right\Vert _{2}} \global\long\global\long\global\long\global\long\global\long \global\long\global\long\def\interior#1{\mathrm{int}\left(#1\right)}

Introduction

Researchers increasingly estimate spillovers in experiments. Sometimes, an intervention is randomized to one set of units while the spillover effects are measured in another set. In such “bipartite” experiments the intervention and outcome units are connected via a network: e.g., of friendship across individuals cai_social_2015 or of supplying relationships across firms best_spillover_2025. In other cases, the intervention and outcome units coincide and direct and spillover effects are jointly estimated in a single experimental sample Miguel2004,chaurey_social_2025. Power is a common issue in such studies: spillover effect estimates can be noisy, while increasing the sample size is costly. It is therefore important to make the best use of the experimental variation, as well as to design the randomization to target the spillover parameter of interest.

This paper develops new tools for experimental design and estimation when spillover effects are of primary interest. Our baseline assumption is that the researcher adopts a linear model of spillovers, in which the key treatment variable is the fraction or number of “friends” (i.e., network connections) who have been randomly selected for the intervention---perhaps augmented with weights reflecting the importance of each connection.\footnote{Specifically, we assume the spillover treatment can be written $x_{i}=\sum_{k}w_{ik}g_{k}$ where $g_{k}$ is an indicator of unit $k$ assigned to the treatment group. The $w_{ik}$ are fully unrestricted, and the intervention units $k=1,\dots,K$ can differ from the outcome units $i=1,\dots,N$. We also discuss how our approach can be used for nonlinear formula treatments $x_{i}=X_{i}(g_{1},\dots,g_{K})$ for known $X_{i}(\cdot)$ that can implicitly depend on some $w_{ik}$ borusyak_design-based_2025.} Such linear specifications are widespread in economics, including in experiments (e.g., Miguel2004,cai_social_2015,atkin_reducing_2024). We assume the researcher has access to the network measure at the experimental design stage. We further assume the researcher wants to identify the spillover effect using the experimental variation; hence, they use some recentered estimator BH1. It is increasingly recognized that recentering is necessary to guarantee a causal interpretation of spillover effect estimates without assuming that the network is exogenous, with appropriate regression or instrumental variable (IV) estimators employed either intuitively Miguel2004 or formally chaurey_social_2025,atkin_reducing_2024.

We then ask: which randomization design and which recentered estimator are likely to yield good power? Our answer reflects two countervailing intuitions. On the one hand, when an outcome unit $i$ is exposed to the shocks of multiple intervention units $k$, it is desirable to positively correlate the treatment assignments of those intervention units. Doing so increases the variance of $i$'s spillover treatment (i.e., the spillover “signal”). On the other hand, when multiple outcome units are exposed to the same intervention unit, their exposures are correlated and this can increase estimation noise if their unobservables are also correlated; correlating assignments of multiple intervention units exacerbates this clustering problem. This makes it more desirable to increase the degree of independence across different exposures, spreading experimental variation out across the sample (“isotropy”). A likely-powerful design is thus one which maximizes spillover signals while ensuring spillover variation is sufficiently diffuse. The estimator choice can also help with isotropy, reweighting the data to spread out the induced spillover variation, at the cost of a lower signal.

We formalize these intuitions in a minimax approach: we look for the experimental design and recentered estimator that minimize the worst-case approximate variance over a large class of error distributions. We parameterize this class parsimoniously, by a scalar that captures the extent to which the researcher is willing to rule out mutual dependencies, dependence on the network, and heteroskedasticity of the errors. This approach has three advantages. First, it reflects the reality that researchers tend to have at best a rough sense of how errors are clustered, especially at the experimental design stage. Second, it allows for general network structures and yields a tractable closed-form characterization of worst-case variance, both as a function of the estimator (holding the design fixed) and as a function of the design (provided the optimal estimator will be used). Although implementing our solution generally involves numerical optimization, it is amenable to practical algorithms. Finally, the minimax approach yields non-degenerate and intuitive solutions in contrast to some corner solutions in the literature (e.g., kiefer_general_1974). For example, independent randomization and cluster-level randomization arise as special cases when, respectively, there are no spillovers and when the exposure is the cluster-level treatment saturation. When exposure is a leave-one-out average of cluster shocks, the solution is to optimally mix between these two designs.

We characterize the minimax-optimal design and estimator in two steps. First, we derive the optimal recentered IV estimator for any experimental design. Its construction entails taking the spillover treatment itself, recentering it by its expectation over the experimental design, and partially “whitening” it to spread out the remaining exogenous variation. While the whitening step weakens the first stage, the induced isotropy improves robustness against non-spherical errors. Second, we characterize the best designs given the optimal recentered instrument will be used. We show this optimal design satisfies natural properties, such as a 50/50 marginal distribution for each shock and separability across independent network blocks. The optimal IV is straightforward to construct for any design, and we propose a computationally efficient approximation to the optimal design; the entire procedure is collected in Algorithm (ref), below. In the parameterization with the highest robustness to adversarial errors, this procedure yields an especially intuitive closed-form solution in which the correlation of shock assignments is determined by the cosine similarity of exposure to them. We show how our estimators are asymptotically normal with a sparse exposure network---even when the errors have strong mutual correlations---and give a simple inference procedure.

We extend these solutions in several directions, including by characterizing minimax-optimal designs and estimators when spillover effects are estimated alongside direct effects or other linear exposure measures. Here the optimal design balances signal and isotropy of the residual variation in the focal spillover treatment after accounting for the other exposures. Other extensions include estimators with predetermined covariates, nonlinear spillover treatments, and optimal design when a researcher faces budget constraints in allocating shocks.

We then illustrate the power gains from using our approach, relative to independent randomization, in semi-synthetic experiments calibrated to the data from the bipartite experiment in cai_social_2015 and the joint estimation of direct and spillover effects in Miguel2004. We find large declines in standard errors in both settings, equivalent to increasing effective sample sizes by around 50--100% or more relative to the independent randomization in the original papers and no whitening. These gains arise from both the optimized experimental design and the use of optimal recentered instruments.

This paper contributes to several related literatures. Most directly, we add to a growing literature on optimal experimental design under interference. Many papers in this literature focus on specific settings, a restricted class of designs, and a pre-specified estimator (with different models of interference and estimands). For example, Baird2018 and Cruces focus on partial interference, where spillovers are confined within clusters, and derive optimal saturation designs for estimating direct and spillover effects, including by regression. pouget-abadie_variance_2019 and harshaw_design_2023 instead focus on bipartite experiments and derive optimal clustered assignments under a linear exposure-response model with heterogeneous treatment effects. thiyageswaran_optimal_2026 study optimal worst-case estimation of the global average treatment effect (GATE) over a general class of designs using the standard Horvitz-Thompson estimator, which does not account for spillovers and is generally biased.\footnote{Other work with general networks includes viviano_experimental_2025 who studies two-wave experiments for estimands including average treatment and spillover effects, using a pilot wave to estimate error variances, and viviano_causal_2025 who choose treatment clusterings to trade off worst-case asymptotic bias and variance when estimating GATEs.}

Relative to this literature, our approach differs in four ways. First, we provide a unified and tractable characterization of optimal experimental designs for general network structures, nesting partial interference and bipartite graphs as special cases. Second, we do not restrict the class of designs a priori. Third, we jointly optimize the design and the estimator to achieve further power improvement. Like harshaw_design_2023, we employ recentering to get unbiased estimators without assuming network exogeneity. Finally, we target the parameters of a linear spillover model, rather than GATE-style estimands. While more restrictive, such constant-effects models are widely used in practice with evidence for linearity in different contexts (e.g., in Miguel2004 and egger2022general).\footnote{In the presence of heterogeneous treatment effects, our baseline estimators identify their convex averages. This is shown in Appendix (ref), in line with earlier results on design-based estimators borusyak_negative_2024.} When both direct and indirect effects are included, our approach focuses on isolating them rather than obtaining the total in a GATE-style analysis. Our problem can also be reformulated to target the GATE; by incorporating additional linear-exposure structure, it is expected to deliver more power than more non-parametric GATE-targeting designs.

We also add to a literature using minimax arguments to justify different experimental designs. In settings without interference, kallus_optimality_2021 and Bai2021 show that complete or blockwise randomization is minimax-optimal absent strong assumptions on error dependence. We extend this logic to spillover estimation with recentered IV, where we show how the minimax-optimal design depends on the allowed adversarial clustering of errors. As in thiyageswaran_optimal_2026, we impose a Schatten $p$-norm constraint on the second-moment matrix of model residuals to define our minimax problem; choosing a smaller $p$ leads to greater robustness to cross-unit dependence and pushes the optimal design towards more diffuse randomization. The ex-ante choice of $p$ by the researcher is similar to choices made in the broader literature on minimax-optimal estimators.\footnote{In non-parametric regression, for example, researchers choose the smoothness class of the data-generating process---determining the tradeoff between robustness and confidence interval length armstrong2018optimal.}

The minimax approach allows us to avoid degenerate solutions that often arise in other analyses. For example, a classic literature on optimal regression design kiefer_general_1974 chooses where to place observations in covariate space to minimize a functional of the ordinary least squares (OLS) covariance matrix (e.g. D-optimality, which minimizes its determinant; kiefer1959optimum). This literature assumes independent and homoskedastic errors, and optimal designs tend to be concentrated on a few extreme points in the design space.\footnote{Our focus on non-spherical errors echoes bickel_robustness_1979 who recover non-degenerate designs when errors are serially correlated.} In the interference setting, the optimal design in pouget-abadie_variance_2019 can similarly produce very coarse randomization: e.g., grouping all units into two clusters only and assigning treatment at the cluster level. harshaw_design_2023 employ a heuristic correction to the optimization criterion to reduce the issue, while our minimax approach resolves it naturally---yielding a guarantee of asymptotic normality with a sparse exposure network and mild regularity conditions.

Finally, we add to a literature on robust and powerful spillover estimation with recentered IV. BH1 show how recentering via knowledge of the (quasi-)experimental design can identify spillover effects under arbitrary endogeneity of the network, while borusyak_optimal_2026 and borusyak_estimating_2025 show how asymptotically optimal recentered instruments can be constructed given a particular design. Here we optimize both the design and instrument, with our applications showing that both levers can meaningfully improve estimator precision. Moreover, while the implementation of the borusyak_optimal_2026 optimal estimator is hindered by its dependence on the (likely unknown) dependence of error terms, our minimax solution provides a parsimonious solution robust to a wide range of error structures.

The remainder of this paper is organized as follows. In Section (ref) we develop our theoretical framework and results. In Section (ref), we propose tractable methods of computing the optimal design and provide a central limit theorem for the recentered IV estimator under the feasible design. In Section (ref) we illustrate these tools in two semi-synthetic experiments. Section (ref) concludes. The proofs of main results are given in Appendix (ref). Additional results are given in Appendix (ref), with proofs in Appendix (ref).

Theory

Setting

We consider estimation of the spillover effects of a set of treatment shocks $g_{k}\in\{0,1\}$, $k=1,\dots,K$, on a set of outcomes $y_{i}$, $i=1,\dots,N$. This might be in a bipartite setting, in which the $K$ intervention units are distinct from the $N$ outcome units, or in a setting like Miguel2004 in which we estimate spillovers within a single experimental sample (with $i=k$).

We parameterize spillovers with a linear exposure model, which we assume is correctly specified:

where the spillover treatment $x_{i}$ combines the shocks $g=(g_{1},\dots,g_{K})^{\prime}$ with exposure weights $w_{i}\in\mathbb{R}^{K}$, which may or may not add to one, and $\varepsilon_{i}$ is an unobserved error.\footnote{Here correct specification includes an implicit exclusion restriction: that the shocks only affect outcomes through $x_{i}$. While we focus on a linear model with constant effects, Appendix (ref) discusses the interpretation of our proposed estimands under heterogeneous treatment effects.} For example, cai_social_2015 estimate a spillover effect $\beta$ from information sessions randomized to a set of $K$ farmers in China on their peers' decisions to adopt weather insurance $y_{i}$; here $x_{i}=w^{\prime}_{i}g$ gives the share of farmer $i$'s neighbors who were assigned to the information session, with $w_{ik}$ being entries of the row-normalized adjacency matrix of the peer graph. In Section (ref) we extend the model to include direct effects of treatment, for settings like Miguel2004. Below we discuss how our main results extend to spillover treatments that are nonlinear functions of the shocks.

In most of our analysis, we condition on the exposure matrix $W=(w^{\prime}_{1},\dots,w^{\prime}_{N})^{\prime}$; we make this conditioning implicit to simplify notation. We allow the error vector $\varepsilon=(\varepsilon_{1},\dots,\varepsilon_{N})^{\prime}$ to be either fixed (as in a standard design-based setup) or drawn from some distribution that can implicitly depend on $W$.\footnote{All of our results hold in both cases but the restriction on the error distribution introduced below is more natural when $\varepsilon$ are viewed as stochastic. Nevertheless, our approach is still design-based in the sense that the moment conditions we leverage derive their validity from the randomization of $g$ alone.} Without loss of generality, we assume that no column of $W$ is entirely zero.

We consider the problem of choosing an experimental design, i.e., a probability distribution for $g$: $\delta\in\mathcal{D}$ where $\mathcal{D}=\Delta(\left\{ 0,1\right\} ^{K})$ and $\Delta(\cdot)$ denotes the simplex. By virtue of randomization, we have $g\protect\mathpalette{\protect\independenT}{\perp}\varepsilon$. We then estimate $\beta$ by recentered IV BH1: i.e., we choose a set of instrument functions $z_{i}(\cdot):\{0,1\}^{K}\rightarrow\mathbb{R}$ where $\expec{\delta}{z_{i}(g)}=0$ for all $i$, and estimate:

We let $\mathcal{Z}_{\delta}$ be the set of all such $z=(z_{i})^{N}_{i=1}$ for a given $\delta$.

Recentering (i.e., the mean-zero property of $z_{i}(g)$, implicitly conditional on $W$) ensures we identify $\beta$ just using the experimental variation in $g$: the network $W$ can be arbitrarily correlated with the unobserved $\varepsilon$. Indeed, the moment condition $\expec{\delta}{\sum^{N}_{i=1}z_{i}\varepsilon_{i}}=0$ is satisfied for any such $z$. Recentered $z\in\mathcal{Z}_{\delta}$ can be obtained from any fixed function, $\tilde{z}_{i}$, of both $g$ and $W$ by subtracting its (known) expectation under a given design: i.e. $z_{i}=\tilde{z}_{i}-\expec{\delta}{\tilde{z}_{i}}$. The set of recentered IV estimators $\hat{\beta}\left[z\right]$, $z\in\mathcal{Z}_{\delta}$, is large, including estimators that weight by functions of $W$ (by rescaling $z$), control for functions of $W$ (by the Frisch--Waugh--Lovell theorem), as well as non-instrumented estimators (e.g., OLS of $y_{i}$ on $\tilde{z}_{i}=x_{i}$ controlling for $\expec{\delta}{\tilde{z}_{i}}$, again by the Frisch--Waugh--Lovell theorem).\footnote{Appendix (ref) considers estimation without recentering. While non-recentered IVs can in principle improve on mean-squared error grounds by introducing a small amount of bias to estimation, we show the corresponding optimal designs can have undesirable properties like unusually high treatment rates. We also show that when $W\mathbf{1}=1$ and an intercept is included in estimation (or, more generally, when $W\mathbf{1}$ is linearly spanned by the included controls), there is no loss in only considering recentered instruments since any bias in non-recentered instruments is absorbed.}

With this setup, we ask: what design $\delta$ and instrument $z$ should we choose in order to get the most precise estimate of the spillover effect $\beta$?

Optimal Design and Instrument

We follow borusyak_optimal_2026 in studying the finite-sample approximate variance of the recentered IV estimator:

where here we let $\mathcal{E}$ denote the (unknown) distribution of $\varepsilon$, with $\V{\delta,\mathcal{E}}z=\infty$ whenever the denominator is zero. Proposition 2 of borusyak_optimal_2026 shows that $\V{\delta,\mathcal{E}}z$ gives a good approximation to the appropriately scaled asymptotic variance of $\hat{\beta}[z]$ so long as the estimator is well-behaved in a particular sense, which primarily means it converges to $\beta$ at some rate with a non-vanishing first stage.\footnote{As a just-identified IV estimator, $\hat{\beta}\left[z\right]$ may have no finite-sample moments. More precisely, in our setting with discrete $g$ the problems may arise from $\hat{\beta}\left[z\right]$ not well-defined in a small-probability event that $\sum^{N}_{i=1}z_{i}(g)x_{i}=0$.}

We look for the experimental design and instrument that minimize the worst-case $\V{\delta,\mathcal{E}}z$ under a weak restriction on the error distribution. Specifically, we suppose the researcher believes $\mathcal{E}$ belongs to some class of distributions $\mathcal{F}_{p}(\sigma)$ and solve the minimax problem:

The class of error distributions, parameterized by $p\in[1,\infty]$ and $\sigma>0$, is given by:

where $\left\Vert \cdot\right\Vert _{p}$ is the Schatten-$p$ matrix norm defined for positive semi-definite matrices as\footnote{For general matrices, it is defined as $\left(\operatorname{tr}\left(\sqrt{A'A}\right)^{p}\right)^{1/p}$. Schatten norms should not be confused with other matrix norms that can be denoted in the same way.} \[ \left\Vert A\right\Vert _{p}=\left(\operatorname{tr}\left(A^{p}\right)\right)^{1/p}. \] Three important special cases are the nuclear norm (a.k.a., the trace norm) with $p=1$, the Frobenius norm with $p=2$, and the operator norm with $p=\infty$. Denoting the eigenvalues of $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$ by $\lambda_{1},\dots,\lambda_{N}\ge0$, we can rewrite ((ref)) as an upper bound on their power mean: \[ \mathcal{F}_{p}(\sigma)=\left\{ \mathcal{E}\colon\left(\frac{1}{N}\sum^{N}_{i=1}\lambda^{p}_{i}\right)^{1/p}\le\sigma^{2}\right\} , \] with the limit $\mathcal{F}_{p}(\sigma)=\left\{ \mathcal{E}\colon\max\left\{ \lambda_{1},\dots,\lambda_{N}\right\} \le\sigma^{2}\right\} $ for $p=\infty$.

To interpret $\mathcal{F}_{p}(\sigma)$, note first that spherical (i.e., homoskedastic and mutually uncorrelated) errors with mean zero (implicitly conditional on $W$) and variance $\sigma^{2}$, i.e. $\expec{\mathcal{E}}{\varepsilon\varepsilon'}=\sigma^{2}I_{N}$, are at the boundary of this set for any $p$. Different choices of $p$ then place different weights on the correlations among the errors. When $p=1$, only the average second moment of $\varepsilon_{i}$ is constrained, allowing fully unrestricted correlations across observations. At the other extreme, when $p=\infty$, correlations are strongly penalized such that a worst-case scenario in ((ref)) is always with spherical errors. Intermediate values of $p$ correspond to smaller penalties on correlations. Figure (ref) provides a visualization for $N=2$ by plotting the boundary of $\mathcal{F}_{p}(\sigma)$ in terms of the two eigenvalues of $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$; increasing the correlation across errors makes the eigenvalues more asymmetric without changing their sum; see also Appendix (ref) for analytical results with equicorrelated $\varepsilon_{i}$. Overall, $p$ captures the extent to which the researcher is willing to rule out strong correlations across the errors when designing the experiment. Similarly, deviations from homoskedasticity and from network exogeneity (i.e., the zero conditional mean of the errors) are penalized more strongly when $p$ is larger.

\paragraph{Optimal Instrument.}

Solving the problem ((ref)) backwards, we first characterize the optimal instrument for a given design $\delta$:

The Appendix (ref) proof provides an explicit characterization of the worst-case error distribution. In this proposition and for the rest of the paper, we use the convention that $1/\infty=0$, $p/(p+1)=1$ for $p=\infty$, and $A^{0}=A^{\dag}A$ for a matrix $A$.

Intuitively, $z^{*}_{\delta}$ takes the spillover treatment $x$, recenters it by $\expec{\delta}x$, and then reweights the resulting $\tilde{x}_{\delta}$ inversely by a fractional power of $\var{\delta}x$, $S^{1/(p+1)}_{\delta}$. The goal in reweighting is to spread out the useful variation in $\tilde{x}_{\delta}$ across different directions that could be adversarially targeted by clustering in $\varepsilon$. When $p=1$ and $\var{\delta}x$ is an invertible matrix, this fully “whitens” $\tilde{x}_{\delta}$ (i.e., rotates it to make $\var{\delta}{z^{*}_{\delta}}=I$) in order to guard against the fully-unconstrained clustering of $\varepsilon$,\footnote{When $\var{\delta}x$ is not invertible, e.g. because $\delta$ perfectly correlates $x_{i}$ for several $i$, $z^{\ast}$ preserves that dependence.} while as $p$ increases we do less whitening as errors with strong mutual correlations are excluded. As a result, at $p=\infty$, the optimal IV is the unwhitened $\tilde{x}_{\delta}$.\footnote{Theorem (ref) technically only shows this for invertible $\var{\delta}x$; the general optimality of $\tilde{x}_{\delta}$ for $p=\infty$, i.e. when $\mathcal{F}_{p}(\sigma)$ bounds the operator norm of the error second-moment matrix, was previously shown in borusyak_optimal_2026.}

\paragraph{Optimal Design.}

We next characterize the optimal design, assuming the optimal instrument will be used for estimation:

This result follows directly from Theorem (ref). For $p=1$, the optimal design maximizes $\operatorname{tr}\left(\var{\delta}x^{1/2}\right)$. At $p=\infty$, the sum of variances, $\operatorname{tr}\left(\var{\delta}x\right)$, is maximized instead.

In what follows we avoid $p=\infty$ as the optimal design problem has a degenerate rank-1 solution. Specifically, Appendix (ref) shows that there exists $\delta^{*}$ that puts all the mass on just one pair of complementary assignments $g^{*}\in\left\{ 0,1\right\} ^{K}$ and $\mathbf{1}-g^{*}$. Moreover, if $W$ is non-negative, one solution corresponds to $g^{\ast}=\mathbf{1}$, i.e. treating all intervention units with probability 0.5 and treating none otherwise. While this is not problematic under spherical errors, it is highly undesirable when errors can be correlated since the estimator may not be consistent. Choosing $p<\infty$ guards against this possibility and helps with consistency, as we show in Section (ref).\footnote{Another problem with degenerate designs is that the asymptotic approximation that makes the approximate variance a useful concept may fail. Moreover, the estimator itself may not be well-defined with a probability that is not asymptotically small, as $\sum_{i}z_{i}(g)x_{i}=0$ whenever $g=\boldsymbol{0}$; see footnote (ref).}

A further result helps characterize optimal designs:

Intuitively, any design with $\expec{\delta^{\ast}}{g_{k}}\ne0.5$ can be symmetrized by averaging it with its complement and this can never hurt the criterion. Similarly, any across-block dependence in $g_{k}$ can get removed without hurting the criterion, since $S_{\delta^{*}}=W\var{\delta^{\ast}}gW^{\prime}$.\footnote{While $\mathcal{F}_{p}(\sigma)$ allows for error dependence even across blocks, Proposition (ref) implies that the same design would be optimal even if one were to a priori rule out cross-block dependence.} This result is generalized in Appendix (ref): whenever $W$ is invariant under a group of relabelings over rows and columns, there exists a $\delta^{*}$ which is invariant under the same group. In particular, if $W$ is invariant to permutations of rows and columns in a cluster, one may restrict attention to exchangeable designs within that cluster.

Special Cases

We now illustrate the results of Theorem (ref) and Proposition (ref) in several insightful special cases.

\paragraph{No Spillovers.}

A useful benchmark is the standard experimental case with $N=K$ and $W=I$. Since $W$ is diagonal, Proposition (ref) immediately implies that independent random assignment, $g_{i}\stackrel{iid}{\sim}Bernoulli(0.5)$, is optimal in this case, and the optimal instrument is just $z^{*}_{i}(g)\propto g_{i}-0.5$.\footnote{Note that complete randomization, with exactly half of the units treated, is not optimal for $p<\infty$ as it induces a slight negative correlation across assignments which can be exploited by the adversary choosing the errors.}

\paragraph{Group-Specific Shocks.}

This and remaining cases are visualized in Figure (ref). Consider a simple bipartite graph in which observations are partitioned into $J$ groups, $C_{1},\dots,C_{J}$, with shocks assigned at the group level: $K=J$ and $x_{i}=g_{j(i)}$ corresponding to $w_{i}=e_{j(i)}$, where $e_{j}$ is the $j$th standard basis in $\mathbb{R}^{K}$ and $j(i)\in\left\{ 1,\dots,K\right\} $ gives the group of observation $i$. Again Proposition (ref) implies independent $g_{j}\stackrel{iid}{\sim}Bernoulli(0.5)$ is optimal. The optimal instrument is now proportional to $N^{-1/(1+p)}_{j(i)}\cdot(g_{j(i)}-0.5)$ where $N_{j}=\left|C_{j}\right|$. Thus when $p=1$ the estimator down-weights large groups by a factor of $1/\sqrt{N_{j(i)}}$ to ensure the random variation is equally spread out; not doing so would make the estimator noisy if the errors are very correlated within clusters. With higher $p$ less of this reweighting is needed, with no reweighting in the limit $p=\infty$.

\paragraph{Group-Average Shocks.}

Now consider a case in which observations are partitioned into $J$ groups $C_{1},\dots,C_{J}$, but the shocks $g_{i}$ are assigned to individuals ($K=N$) and fully propagate through the cluster-average: $x_{i}=\bar{g}_{j(i)}$ for $\bar{g}_{j}=\frac{1}{N_{j}}\sum_{i\in C_{j}}g_{i}$, which corresponds to $W=\operatorname{diag}(P_{N_{1}},\dots,P_{N_{J}})$ for $P_{n}=\mathbf{1}_{n}\mathbf{1}^{\prime}_{n}/n$. This is isomorphic to the previous case, with shocks that are perfectly correlated within groups and independent across groups being optimal: $g_{i}=\check{g}_{j(i)}$ with $\check{g}_{j}\stackrel{iid}{\sim}Bernoulli(0.5)$. The optimal instrument is again proportional to $N^{-1/(1+p)}_{j(i)}(\check{g}_{j(i)}-0.5)$. Intuitively, there is no need to vary the shocks within groups since treatment is at the group level; perfect within-group shock correlation ensures the maximal variability of the spillover treatment while we still reweight by $N^{-1/(1+p)}_{j(i)}$ to more evenly spread out this variation.

\paragraph{Leave-Out Averages.}

Suppose in the previous case we instead have $N_{j}>1$ for all $j$ and $x_{i}=\bar{g}_{-i}=\frac{1}{N_{j(i)}-1}\sum_{k\in C_{j(i)},k\neq i}g_{k}$. Appendix (ref) shows that an optimal design mixes between group-specific iid $Bernoulli(0.5)$ shocks with probability $\psi_{j}=(1-r_{j})/(1+(N_{j}-1)r_{j})$ for $r_{j}=(N_{j}-1)^{-2p}$, and independent individual-level $Bernoulli(0.5)$ shocks with probability $1-\psi_{j}$, independently across groups.\footnote{Note that this is just one way to induce the optimal $\var{\delta^{*}}x$. Other designs may deliver the same variance matrix and are therefore also optimal.} This design is close to a more conventional saturation design (studied, e.g., in Baird2018) with saturation rates of $0$, $0.5$, and $1$, except with independent (rather than complete) randomization within groups. The optimal instrument is now proportional to $\left(1+(N_{j(i)}-1)\psi_{j(i)}\right){}^{-1/(1+p)}\cdot\left((\bar{g}_{j(i)}-\frac{1}{2})+(N_{j(i)}-1)(\bar{g}_{j(i)}-g_{i})\right)$. It is easy to verify that it is optimal to correlate the shocks more in larger groups: $\psi_{j}$ increases in $N_{j}$, from $\psi_{j}=0$ when $N_{j}=2$ to $\psi_{j}\to1$ when $N_{j}\to\infty$. Intuitively, in the smallest group of size $N_{j}=2$, we are effectively back in the no-spillover case, while with $N_{j}\to\infty$ the leave-out average is equivalent to the group average from the previous case. For intermediate sizes, within-group correlation increases in $p$: as before, it is higher when less clustering is allowed in $\varepsilon$.

\paragraph{Spatial Spillovers.}

Finally, suppose the $N=K$ units are equal-spaced points on a circle, with treatment being an average of shocks to units within the distance of $L$: $x_{i}=\frac{1}{2L+1}\sum^{i+L}_{k=i-L}g_{k}$, where index summation is understood modulo $N$. Appendix (ref) shows that, for $p=1$, the optimal design correlates shock assignments at the same spatial scale as the exposure mapping. Specifically, the correlation between assignments at circular distance $d$ is $\max\left\{ 1-d/(2L+1),0\right\} $, linearly decaying in $d$. The design can be implemented by splitting the circle into blocks of $2L+1$ consecutive units, randomly shifting this partition, and assigning shocks at the block level. The optimal instrument is $z^{\ast}=W^{\dag}\tilde{g}$. Under the shifted-block implementation of the design, it takes values $(2L+1)\tilde{g}_{i}-2L\bar{\tilde{g}}$ for units at the center of a block and $\bar{\tilde{g}}$ for all other units, where $\bar{\tilde{g}}$ is the average of all recentered shocks. Intuitively, the IV undoes the smoothing done by spatial averaging to mitigate the correlations that could otherwise be exploited by the adversary. These results can be extended to a two-dimensional grid, with $x_{i}$ averaging the shocks within a certain distance in each coordinate.

Multiple Exposures

We now consider a model in which the shocks $g$ affect the outcomes through multiple linear exposure measures. For brevity, we consider two exposures but all results generalize immediately to three or more. For a known fixed $N\times K$ matrix $U=(u_{i}')^{N}_{i=1}$, we now have:

Our goal is still to minimize the worst-case approximate estimation variance of the effect of $x_{i}$.

While we consider a general formulation, this setting is especially relevant when the intervention and outcome units are the same and $\tau$ corresponds to the direct effect of the intervention:

which corresponds to $U=I_{N}$.\footnote{Another class of applications is to multiplex networks (cf. zenou2025peer).} In this context, there are three reasons why one would focus on the estimation power of $\beta$. First, the researcher may be most interested in the spillover effect, but still wants to acknowledge the presence of the direct effect. Second, even if direct and spillover effects are of equal interest, the spillover effect is usually more challenging to estimate. A design that yields sufficient power for it would usually deliver precise estimates of the direct effect, too. Finally, if the researcher is interested in the Global Average Treatment Effect (GATE), defined as the effect of switching from $g=\boldsymbol{0}_{K}$ to $g=\boldsymbol{1}_{K}$ via both direct and indirect effects, that problem can be reformulated as the one we solve here.\footnote{Specifically, the GATE equals $\gamma=\bar{w}\beta+\bar{u}\tau$ for known constants $\bar{w}=\frac{1}{N}\sum_{i,k}w_{ik}$ and $\bar{u}=\frac{1}{N}\sum_{i,k}u_{ik}$. One can then rewrite ((ref)) as $y_{i}=\gamma\frac{x_{i}}{\bar{w}}+\tau\left(h_{i}-\frac{\bar{u}}{\bar{w}}x_{i}\right)+\varepsilon_{i}$ and apply our results to this reformulation, optimizing estimation power for $\gamma$ while also including $h_{i}-\frac{\bar{u}}{\bar{w}}x_{i}$ in the model.}

To see how Theorem (ref) and Propositions (ref)--(ref) generalize, collect the parameters and explanatory variables in $\theta=(\beta,\tau)^{\prime}$ and $\boldsymbol{x}_{i}=(x_{i},h_{i})^{\prime}$. Absent restrictions on how the error term can correlate with $W$ and $U$, we again use recentered IV. For some recentered instrument vector $\boldsymbol{z}(\cdot)=(\boldsymbol{z}_{i}(\cdot))^{N}_{i=1}$, where $\boldsymbol{z}_{i}(\cdot):\{0,1\}^{K}\rightarrow\mathbb{R}^{2}$ and $\expec{\delta}{\boldsymbol{z}_{i}(g)}=0$, consider the estimator:

Define the approximate variance of $\hat{\beta}$ as:

where $e_{1}=(1,0)^{\prime}$, with $\boldsymbol{\mathcal{V}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right]=\infty$ whenever $\expec{\delta}{\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)\boldsymbol{x}^{\prime}_{i}}$ is singular. Also, as before, define

An initial result shows that the problem of searching over recentered IV vectors $\boldsymbol{z}$ is equivalent to a constrained version of the original problem of searching over scalar $z$ that requires $z$ to be orthogonal to the “nuisance” exposure $h$:

Imposing the orthogonality constraint yields the analog of Theorem (ref):

The key difference from the optimal instrument in the single-exposure case is that the matrix $W$ is replaced with an appropriate residualization $W-c_{\delta}U$ where $c_{\delta}$ (which plays the role of a Lagrange multiplier on the $\expec{\delta}{z^{\prime}h}=0$ constraint) can be solved for by the scalar optimization ((ref)). The additional technical assumption for $p=1$ is imposed because the objective ((ref)) can be non-differentiable, which may lead to a different analytical form of $z^{*}_{\delta}(g)$. To avoid this issue, we require that the residualized exposure matrix retains as much assignment-induced variation as its dimensions allow. Given orthogonality between the optimal IV for $x_{i}$ and the exposure $h_{i}$, the IV for $h_{i}$ is simply $h_{i}-\expec{\delta}{h_{i}}$.

Analogs of Propositions (ref)--(ref) also follow:

\setcounter{prop}{2} \setcounter{thm}{1}

We note that, unlike with Proposition (ref), the optimal design cannot generally be solved block-by-block in the block-diagonal case because $c_{\delta}$ is common across all blocks.

We illustrate these results in a simple but nontrivial special case: when $x_{i}$ is the leave-out average of shocks in the cluster to which $i$ belongs and $h_{i}=g_{i}$ captures the direct effect. Appendix (ref) shows that, at least when $N_{j}$ is even, it is optimal to mix between cluster-specific $Bernoulli(0.5)$ shocks and complete randomization of the shocks within each cluster, independently across clusters. In the special case where all clusters are of the same size, $N_{j}=n$, one can obtain a solution more similar to the one from Section (ref): mix between cluster-specific shocks with probability $\psi^{\perp}=(1-r^{\perp})/(1+(n-1)r^{\perp})$ for $r^{\perp}=(n-1)^{-2p/(2p-1)}$ and independent $Bernoulli(0.5)$ randomization with probability $1-\psi^{\perp}$. The optimal spillover instrument is proportional to $(\bar{g}_{j(i)}-0.5)+(n-1)^{1/(2p-1)}(\bar{g}_{j(i)}-g_{i})$. Relative to the case without direct effects, there is generally more randomization within clusters ($\psi^{\perp}\le\psi$, with strict inequality for $n>2$ and $p>1$), which helps isolate spillover effects from the direct effects. The exception is when $p=1$: with both homogeneous and heterogeneous $N_{j}$, the design and instrument are then the same as in the original problem, because the objective already creates enough variation to identify direct effects.

Other Extensions

We derive three further extensions to the baseline results. First, Appendix (ref) shows how the results extend when the researcher includes an intercept and possibly other predetermined covariates $r_{i}$ in estimation. We formalize this idea by enriching the class of error distributions to $\varepsilon_{i}=r^{\prime}_{i}\gamma+\tilde{\varepsilon}_{i}$ where $\tilde{\varepsilon}=\left(\tilde{\varepsilon}_{i}\right)$ satisfies $\expec{\mathcal{E}}{\tilde{\varepsilon}\tilde{\varepsilon}'}\in\mathcal{F}_{p}(\sigma)$ while $\gamma$ is unrestricted. The results in Theorem (ref) and Propositions (ref)--(ref) continue to hold, with the exposure matrix $W$ residualized columnwise on the covariates.

Second, Appendix (ref) discusses how Theorem (ref) and Proposition (ref) extend to the general formula treatment setting where the $x_{i}$ are nonlinear (but still known) $i$-specific formulas combining $W$ and $g$. This allows, for example, $x_{i}$ to be an indicator for $i$ having at least one treated friend in their social network or the more elaborate formulas considered in BH1. The proofs to Theorem (ref) and Proposition (ref) turn out to extend verbatim after defining the general $\tilde{x}_{\delta}=x-\expec{\delta}x$ and $S_{\delta}=\var{\delta}x$, making them a robust characterization of optimal design with formula treatments. Proposition (ref), however, does not extend cleanly: it is no longer without loss of generality to consider designs with $0.5$ marginals. The following discussion of feasible design computation and inference is also reliant on linearity of the spillover treatment formula.

Finally, Appendix (ref) considers settings where the researcher faces a constraint of a marginal treatment probability $q\neq0.5$ (common across units), e.g. due to a fixed budget or government target. Theorem (ref) applies to any design and is unaffected by the constraints on feasible designs. We therefore show that Proposition (ref) extends naturally with such constraints, as well as the second claim of Proposition (ref) and the results on computation, below.

Feasible Designs and Inference

We now introduce a relaxed version of the optimal design problem. This relaxation admits a closed-form solution in special cases and a computationally efficient approximation in general. Moreover, this approximation can yield $\sqrt{K}$-consistent recentered IV estimators with $p<\infty$, as well as a simple asymptotic inference procedure.

The Relaxed Problem

Solving the Proposition (ref) problem is computationally intractable for even moderate $K$. While its objective depends on the experimental design $\delta$ only through the $K\times K$ covariance matrix $\var{\delta}g$, the class of all $K\times K$ covariance matrices achievable with binary treatments, $\mathcal{C}_{K}=\left\{ \var{\delta}g\colon\delta\in\mathcal{D}\right\} $, is high-dimensional. The computation time required to find an optimal $\var{\delta}g$, by searching over convex combinations of the $2^{K}$ possible treatment vectors in $\left\{ 0,1\right\} ^{K}$, is generally exponential in $K$.

To make progress, we follow goemans_improved_1995 and thiyageswaran_optimal_2026 in considering a relaxed problem with a polynomial-time solution in $K$. Specifically, we replace the optimization over $\mathcal{C}_{K}$ with optimization over $\mathcal{Q}_{K}=\left\{ \Sigma\in\mathbb{S}^{K}:\Sigma\succeq0,\Sigma_{kk}=1/4\text{ for all }k\right\} $, where $\mathbb{S}^{K}$ is the set of $K\times K$ real and symmetric matrices. Thus $\mathcal{Q}_{K}$ is the set of $K\times K$ covariance matrices with $1/4$ on the diagonal, containing $\mathcal{C}_{K}$ when we restrict attention to designs with $\expec{\delta}{g_{k}}=0.5$ (without loss of generality, by Proposition (ref)). The relaxed Proposition (ref) problem is then:

Searching over $\mathcal{Q}_{K}$ is feasible for moderate $K$ using an interior point solver for conic optimization problems. Two additional results further improve computational efficiency without any additional approximation cost. First, if $W$ can be split into independent blocks $W_{(1)},\dots,W_{(B)}$ with $K(b)$ shocks in each, we can use Proposition (ref) to solve the optimal design problem separately by block. Second, as the following result shows, we can characterize the solution explicitly in terms of a lower-dimensional convex optimization problem over the $K-1$-dimensional simplex $\Delta_{K}=\left\{ \ell\in\mathbb{R}^{K}\colon\ell_{k}\ge0,\sum^{K}_{k=1}\ell_{k}=1\right\} $ rather than $\mathcal{Q}_{K}$; this problem is amenable to gradient-based optimization approaches that scale more efficiently with $K$.

Corollary (ref) captures a simple intuition: shocks to two intervention units should be correlated to the extent that they are connected, in the sense that the same outcome units are exposed to them. This turns out to be the optimal solution to the relaxed problem when $p=1$. Corollary (ref) extends this result to arbitrary $p$, as long as each block of $W$ is sufficiently symmetric---a restrictive condition that nevertheless includes all of the special cases in Section (ref). The corresponding solution is particularly intuitive when $W$ is nonnegative and $p$ is an integer. In that case, the matrix $(W'W)^{p}$ captures $p$th-degree connections between intervention units $k$. When $p=2$, for example, the solution suggests correlating two shocks when they are either connected directly or connected to the same third shock. The worst-case loss from such longer-range shock correlations is lower with higher $p$, as error correlations are more restricted. As $p\to\infty$, this solution correlates all shocks in the same connected component of the bipartite network defined by $W$.

Corollary (ref) further characterizes the IV that is optimal if the relaxed problem's solution can be implemented (an issue we discuss below). This optimal IV satisfies $W'z^{\ast}\propto\tilde{g}$, meaning that partial whitening entails the inverse operation to the exposure mapping, $\tilde{x}=W\tilde{g}$. For instance, when $x$ involves some averaging of the (correlated) shocks, $z^{\ast}$ performs their deconvolution.\footnote{Interestingly, the optimal IV is invariant to $p$, up to the proportionality constant $a$. With larger $p$, the shocks are more strongly correlated; thus, the same IV that fully whitens the shocks when $p=1$ does only partial whitening for larger $p$. We also note that when $W$ consists of independent blocks the constant $a$ varies across blocks; thus, $p$ affects the reweighting of blocks but not the optimal IV within blocks.}

All of these results generalize immediately to settings with an intercept and possibly other predetermined covariates in estimation, by first residualizing $W$ on those covariates column by column. For instance, for $p=1$ and when estimation includes an intercept as the only covariate, Corollary (ref) applied to the demeaned $w_{\cdot k}$ implies that shock correlations under the relaxed problem's solution are equal to the correlations of the $N\times1$ exposures to those shocks (rather than cosine similarities without the intercept). This solution positively correlates shocks with similar exposure while introducing slight negative correlations between shocks with non-overlapping exposure.\footnote{That is, if $w_{\cdot k},w_{\cdot l}\ge0$ and $w_{\cdot k}'w_{\cdot l}=0$ (no exposure overlap), $\Sigma^{\ast}_{kl}$ is proportional to the sample correlation of $w_{\cdot k}$ and $w_{\cdot l}$, which is negative.}

We can also generalize the relaxed problem to multiple exposures, where it remains computationally efficient. By the same arguments as in the Proposition (ref) proof, the relaxed problem is:

Here the inner problem has the same structure as before with a single exposure, with $W$ replaced by $W-cU$; the computationally efficient form of Proposition (ref) and its closed-form special cases thus apply.\footnote{The Proposition (ref) problem becomes $\min_{c\in\mathbb{R}}\min_{\ell\in\interior{\Delta_{K}}}\operatorname{tr}\left\{ \left(D^{-1/2}_{\ell}\left(W-cU\right)'\left(W-cU\right)D^{-1/2}_{\ell}\right)^{p}\right\} $.} Additionally, when $W$ and $U$ are jointly block-diagonal, the inner problem can be split block-by-block for any $c$. The outer problem is then a scalar optimization, amenable to a golden section search over a compact subset of $\mathbb{R}$.\footnote{The minimizer is guaranteed to exist in a compact subset of $\mathbb{R}$ as long as $\max\limits_{\Sigma\in\mathcal{Q_{K}}}\operatorname{tr}\Big(\left((W-cU)\Sigma(W-cU)')\right)^{\frac{p}{p+1}}\Big)\ge\operatorname{tr}\Big(\left((W-cU)I/4(W-cU)')\right)^{\frac{p}{p+1}}\Big)\to\infty$ when $\left|c\right|\to\infty$, which requires only that $U\neq0$.}

Feasible Implementation

Two challenges remain once a solution $\Sigma^{\ast}$ to the relaxed problems in Equation ((ref)) or ((ref)) is found. First, this $\Sigma^{\ast}$ may not correspond to an element of $\mathcal{C}_{K}$: i.e., it may not be implementable via binary shock vectors. Second, even if it is theoretically implementable, there is no known computationally feasible procedure for sampling binary vectors from a generic covariance matrix. These problems can be overcome in simpler special cases: the designs in Section (ref) are indeed implementations of the relaxed problem solution.

Outside special cases, we address both challenges by using a simple Gaussian rounding procedure to sample binary assignments with the variance matrix approximating $\Sigma^{\ast}$, following goemans_improved_1995. Specifically, we draw $\xi\sim\mathcal{N}(0,\Sigma^{\ast})$ and set $g^{GR}_{k}=\mathbf{1}\left[\xi_{k}\ge0\right]$.\footnote{This procedure simplifies further in the special case of Corollary (ref): draw $\eta\sim\mathcal{N}(0,I_{N})$ and set $g^{GR}_{k}=\mathbf{1}\left[(W'\eta)_{k}\ge0\right]$. Indeed, we can write $g^{GR}_{k}=\mathbf{1}\left[\xi_{k}\ge0\right]$ for $\xi_{k}=\sum_{i}w_{ik}\eta_{i}/2\left\Vert w_{\cdot k}\right\Vert _{2}$, where $(\xi_{1},\dots,\xi_{K})\sim\mathcal{N}(0,\Sigma^{\ast})$ for $\Sigma^{\ast}$ from Corollary (ref).} This procedure distorts correlations in a known way: $\corr{}{g^{GR}_{k},g^{GR}_{l}}=\arcsin\left[\corr{}{\xi_{k},\xi_{l}}\right]/\frac{\pi}{2}\ne\corr{}{\xi_{k},\xi_{l}}$; see goemans_improved_1995 and the illustration in Appendix Figure (ref). Given this Gaussian rounding design, we find the optimal IV based on $S_{GR}=W\var{}{g^{GR}}W^{\prime}$ for

with $\arcsin\left[\cdot\right]$ applied entrywise. Algorithm (ref) summarizes the entire procedure.

Beyond simplicity in implementation, this procedure offers two advantages. First, Appendix (ref) shows the cost of using this relaxed and rounded version of the problem is bounded: the worst-case approximate variance of the resulting estimator is at most $\pi/2$ times that of the original problem.\footnote{We find much smaller approximation errors in applications: for all designs in Section (ref), the worst-case approximate variance exceeds the variance of the original problem by no more than 6%. This is computed by comparing the worst-case variance arising from the Gaussian rounded design to that of the unrounded solution to the relaxed problem ((ref)), which is a lower bound for the variance of the optimal solution to the original problem.} Second, as we next show, the fact that our shocks are a simple transformation of a correlated Gaussian vector is helpful for establishing a central limit theorem for the recentered IV estimator and asymptotically valid inference, even when all errors are mutually correlated.

Asymptotic Normality and Inference

We now consider the asymptotic behavior of the estimator $\hat{\beta}=\hat{\beta}\left[z^{\ast}\right]$ corresponding to the Gaussian rounding implementation of the Section (ref) relaxed-optimal design, with the Theorem (ref) optimal IV $z^{\ast}$, as in Algorithm (ref). We show that choosing a small integer $p$ can yield sparsity of the implied shock covariance matrix, which in turn can yield a central limit theorem for $\hat{\beta}$ and a simple inference procedure.

Let $d(\Sigma)=\max_{k}\sum^{K}_{l=1}\mathbf{1}\left\{ \Sigma_{kl}\neq0\right\} $ measure the sparsity of a covariance matrix $\Sigma$ by the maximum number of non-zero covariances across all rows. We consider a sequence of data-generating processes indexed by $K$. For $K\to\infty$ (which under Assumption (ref) will require $N\to\infty$ as well) we assume:

The first part of Assumption (ref)(a) is the key sparsity condition imposed on the solution to the relaxed problem in Equation ((ref)). While the sparsity of $\Sigma^{*}$ is formally an asymptotic condition, the researcher can heuristically check whether $d(\Sigma^{\ast})$ is small compared to $K$; at the cost of more complex primitive conditions, we can accommodate $d(\Sigma^{\ast})$ growing slowly in $K$.\footnote{Theorem (ref) can be extended to cases where $\Sigma^{*}$ is approximately sparse, in that it contains many small entries but has row and column sums that are bounded or grow very slowly. This extension is particularly useful when predetermined covariates are included, since then the effective $W$ and corresponding $\Sigma^{*}$ will not generally be sparse. In such cases we could allow the operator norm of $\Sigma^{*}$ to grow slowly by a smoothing argument and the second order Poincaré inequality of chatterjee2009fluctuations, as applied in lei2018asymptotics. We thank Lihua Lei for pointing this out.} Since $\frac{1}{N}\left\Vert x\right\Vert ^{2}_{2}=\frac{1}{N}\left\Vert Wg\right\Vert ^{2}_{2}\le\lambda_{\max}(W'W)\cdot\frac{1}{K}\left\Vert g\right\Vert ^{2}\cdot\frac{K}{N}$, the second part precludes $x$ from diverging provided $K\asymp N$. Appendix (ref) provides sufficient conditions for both parts of this assumption. It shows that $d(\Sigma^{\ast})$ is upper-bounded by $\left(d(W'W)\right)^{p}$ under the assumptions of Proposition (ref). With $p=\infty$, no sparsity of $\Sigma^{\ast}$ is guaranteed within any connected component of the graph, and so designs that are optimal for $p=\infty$ may not yield consistent $\hat{\beta}$ when $W$ does not consist of many independent blocks and when errors are also strongly dependent. But for smaller integer $p$, the optimal design only correlates units that are connected by paths in $W'W$ of lengths $p$ or less. Thus, for finite integer $p$, both parts of Assumption (ref)(a) are guaranteed when maximum row and column degree of $W$ as well as the maximum absolute value of its elements are all bounded. This set of conditions in turn requires $K\asymp N$ in non-degenerate cases, as also shown in Appendix (ref).

In addition to sparsity, Assumption (ref) imposes regularity conditions to ensure the joint distribution of spillover treatments $x_{i}$ and errors $\varepsilon_{i}$ is well-behaved. Assumption (ref)(b) precludes near-collinearity of exposures to ensure that the whitening matrix remains well-conditioned. Yet, this condition allows $S_{GR}$ to be degenerate, which happens, e.g., when some observations have the same exposure to all shocks. Noting that $h_{K}=\expec{\delta}{\frac{1}{K}x'z^{\ast}\mid W}$, Assumption (ref)(c) implies that $W$ and the experimental design are such that there is sufficient variation in the spillover treatments $x_{i}$ so that the first stage is not degenerate asymptotically. Noting further that $v_{K}=\frac{1}{K}\var{\delta}{\varepsilon'z^{\ast}\mid W,\varepsilon}$, Assumption (ref)(d) ensures that the finite-sample variance of $\varepsilon'z^{\ast}$ conditional on $W$ and $\varepsilon$ (scaled appropriately) converges to a deterministic constant. This is a law of large numbers for the quadratic form $v_{K}$ and it imposes regularity on the joint distribution of $W$ and $\varepsilon$. For example, conditional on a sequence of $W$, if $\varepsilon$ is drawn from a mixture of two deterministic sequences that lead to different limits for $v_{K}$, then this assumption is violated; when dependence in $\varepsilon$ is sufficiently weak or local, then this assumption is satisfied. Assumption (ref)(e) is an anti-concentration condition for the implicit weights of the numerator $\frac{1}{K}\varepsilon'z^{\ast}$; combined with the sparsity condition, it ensures that the numerator can be asymptotically approximated by sums of independent components.

Under these conditions, $\widehat{\beta}$ is $\sqrt{K}$-consistent:

Asymptotically valid inference is straightforward from this result using the normal approximation and a simple plug-in estimator for the standard error, $\sqrt{\hat{v}/(Kh^{2}_{K})}$. Here $S_{GR}$ and therefore $h_{K}$ are known. Moreover, the plug-in estimator for $v^{\ast}$ is straightforward to compute: $\hat{v}=\frac{1}{K}\hat{\varepsilon}\,'S^{(p-1)/(p+1)}_{GR}\hat{\varepsilon}$ for $\hat{\varepsilon}=y-\hat{\beta}x$. The convergence in probability of $\hat{v}$ to $v^{*}$ follows from consistency of $\widehat{\beta}$ and Assumption (ref) (see Appendix (ref)).

Applications

Bipartite Experiment: Cai et al. (2015)

\paragraph{Setup.}

Our first application is the bipartite experiment of cai_social_2015 who study how social networks affect the adoption of weather insurance by randomly assigning rice farmers in villages in China to simple vs. intensive information sessions over two rounds. They collect information on demographics, friendship, and insurance take-up for each farmer. We calibrate a simulation to a simplified version of the social network effect specification in Table 2, Column 2 of their paper. For the sample of farmers assigned to either the simple or intensive information session in the second round, it examines how insurance take-up varies with the fraction of friends that are treated in the first round. Specifically, our outcome units $i$ are those randomized in the second round ($N=1,274$),\footnote{More precisely, outcome units are those used to estimate spillover effects on insurance take-up (from “Type 1” villages) rather than price effects.} and our intervention units $k$ are those randomized in the first round with at least one friend in the second round ($K=995$). As in cai_social_2015, we define $x_{i}=\sum_{k}w_{ik}g_{k}$ where $w_{ik}$ equals an indicator that $i$ and $k$ are friends divided by the total number of friends $i$ has (including friends outside the first-round sample). Thus, $x_{i}$ measures the percentage of the second round unit's friends who were treated in the first round. Each $i$ and $k$ is assigned to one of 44 administrative villages, which are administrative units that each combine multiple local communities; throughout this analysis we drop 10 of 6,251 friendships that are across administrative villages so $W$ has block-diagonal structure that eases computational burden. The top row of Figure (ref)(a) shows the exposure matrix $W$ for one example village, along with the $W'W$ matrix that captures similarity in exposure weights across intervention units.

We estimate the specification $y_{i}=\alpha+\beta x_{i}+\upsilon_{i}$ using recentered IV without whitening. Our simple specification has a smaller estimated coefficient of $\beta^{\ast}=0.1408$ compared to the specification in cai_social_2015, which uses OLS as an estimator and includes an additional set of controls and village fixed effects, rather than recentering.

For $S=1000$ simulations, we generate the outcome as $y_{i}=\beta^{*}x_{i}+\varepsilon_{i}$, where $g_{k}$ varies by the experimental design, $W$ is fixed, and we draw $\varepsilon=(\varepsilon_{i})$ from six different data-generating processes (DGPs), i.e., joint distributions that vary the structure of the $N\times N$ matrix $\expec{\mathcal{E}}{\varepsilon\varepsilon^{\prime}}$.\footnote{While the original outcomes in cai_social_2015 are binary, our generated outcomes are continuous.} For the first five DGPs, the distribution of the errors is calibrated so that the standard error of the unwhitened recentered IV estimator of $\beta$ under a 50% Bernoulli RCT matches the standard error of the estimated coefficient in the data. The error distributions for the simulation are as follows:

\paragraph{Results.}

Following Algorithm (ref), we solve for the feasible optimal design and recentered IV estimator of the indirect effect for several choices of $p$. The bottom row of Figure (ref)(a) plots a correlation heatmap for these designs. As our theory suggests, $p=\infty$ involves perfectly correlating all treatments within clusters of connected farmers (with independent randomization across villages), while lower $p$ implies more diffuse correlation structures.\footnote{Appendix Figure (ref) shows how the optimal designs in both applications trade off isotropy, as measured by $\operatorname{tr}(S_{\delta^{*}})/\lambda_{\max}(S_{\delta^{*}})$, and signal $\operatorname{tr}(S_{\delta^{*}})$. As expected, the frontier moves down and to the right as $p\rightarrow\infty$. Notably, the baseline independent Bernoulli design (indicated by the red RCT dot) has higher isotropy than any optimal design but also much lower signal.}

Table (ref) reports root mean-squared error (RMSE) for the optimal recentered IV estimator from Theorem (ref); RMSE represents the standard deviation of the recentered IV estimator since the bias is negligible. All estimators include an intercept.\footnote{Correspondingly, Algorithm (ref) is applied with $R=\mathbf{1}$. We use the optional part of step 2 to block-diagonalize the residualized $W$. This allows us to solve for the optimal design village by village, with independent randomization across villages.} The columns correspond to the error structures described above, while the rows correspond to the optimal designs, along with a benchmark of a simple Bernoulli experiment with a 50% treatment probability and a recentered IV estimator without whitening. In brackets, we also include the gain in effective sample size, defined as the squared RMSE ratio relative to the benchmark.

All optimal designs beat a standard RCT, which is not optimized to estimate spillover effects. Under homoskedastic errors, the design with $p=\infty$ (i.e., perfectly correlating treatments within blocks) works best. But under more complex and realistic error structures, the advantage of extreme shock correlations deteriorates: the $p=2$ design is best for the DGP in columns 3--5 with network-dependent or correlated errors. This choice of $p$ also shows robust power gains across all columns, equivalent to increasing sample size of 41--75% in all columns.

To isolate the gain from the optimal design vs. the optimal estimator, Appendix Table (ref) reports RMSE of the recentered IV estimator without whitening for the $p$-optimal designs. On one hand, the $p=\infty$ design involves no whitening yet achieves substantial improvements across most error distributions---suggesting a key role of design. Moreover, for simple (e.g., homoskedastic) errors whitening only reduces precision gains. On the other hand, column 5 shows that with more challenging error distributions extra power gains due to the $p=1$ and $p=2$ designs are only realized if the instrument is also whitened. Appendix Table (ref) also includes an analysis of coverage for confidence intervals built using the normal approximation from Theorem (ref). Since the network of cai_social_2015 is made up of many independent villages, all designs and simulations have close to nominal coverage, even when shocks are perfectly correlated within villages under $p=\infty$.

Direct and Indirect Effects: Miguel and Kremer (2004)

\paragraph{Setup.}

Our second application is in the setting of Miguel2004 who estimate the direct and spillover effects of a randomized school-level deworming treatment on helminth infection rates among children in Kenya. The intervention shocks $g_{i}$ are at the school level; slightly simplifying the original analysis, we define the outcome $y_{i}$ as the average infection rate in school $i$. Thus, the outcome and intervention units are the same set of $N=K=49$ schools (specifically, Groups 1 and 2 of the experiment for which health outcomes are observed). Following Miguel2004, the spillover treatment $x_{i}$ is the number (rather than the share) of students in nearby schools that were treated, rescaled by 1,000. This is formalized by the $N\times N$ adjacency matrix $W$ with entry $w_{ik}=E_{k}/1000$ if school $i$ is within 3km of school $k\ne i$, where $E_{k}$ is the number of pupils eligible for the treatment in school $k$, and 0 otherwise (including $w_{ii}=0$ for all $i$). The top row of Figure (ref)(b) shows the exposure matrix $W$ and the matrix of similarity in exposure weights $W'W$.

We calibrate our simulation to a simplified version of column (1) in Table VII of Miguel2004. We estimate a school-level regression $y_{i}=\alpha+\beta x_{i}+\tau g_{i}+\upsilon_{i}$ using recentered IV without whitening. This yields coefficients $\beta^{\ast}=-0.226$ and $\tau^{\ast}=-0.258$.\footnote{Our specification differs from Miguel2004 in several ways. We use a linear model rather than probit, we recenter rather than including additional controls, and we drop a third treatment variable based on exposure to schools between 3 and 6 km away. Despite these simplifications, our estimates are very similar to the corresponding average marginal effects in the original paper, of -0.26 for the indirect effect and -0.25 for the direct effect.}

For $S=1,000$ simulations, we generate the outcome as

where again the distribution of $g$ depends on the experimental design, $W$ is fixed, and we vary $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$ according to the same six DGPs described above, with a few minor changes. We replace the village-correlated DGP with a “cluster-correlated” one by creating five clusters of schools using spectral clustering on $W$ and use random cluster shocks. For the network-correlated DGP, we set $W_{\text{raw}}=W$. We calibrate the simulations like in the first application: to match the estimated standard error of the unwhitened recentered IV estimator and targeting the same shares of variation explained by different error components.\footnote{The corresponding parameters are{ as follows. Homoskedastic: $\nu=0.17$; heteroskedastic: $\gamma=0.4$, $\nu=0.12$; degree-in-mean: $\gamma=0.056$, $\nu=0.085$, cluster-correlated: $\nu=0.10$; network-correlated: $\gamma=0.12,\nu=0.05$.}}

\paragraph{Results.}

We follow Algorithm (ref) to solve for the feasible design optimized for the indirect effect in a linear model ((ref)) that also includes a direct effect. As in the first application, we include the intercept for all estimators.\footnote{We again use the algorithm with $R=\mathbf{1}$. Here, however, we do not need the optional part of step 2.}

The bottom row of Figure (ref)(b) presents the correlation heatmap for the optimal designs used in the simulations. Now that the direct effect is in the model, we no longer have a perfectly correlated design when $p=\infty$. As $p$ decreases, the optimal design from Proposition (ref) weakens both positive and negative correlations of the assignments to protect against adversarially correlated (again, positively or negatively) error distributions.

The RMSE results for the indirect and direct effects are in Table (ref), including in brackets the effective sample size gain compared to the baseline RCT design with an unwhitened recentered IV estimator. For estimating the spillover effect, Panel (a) shows that all optimal designs beat the simple RCT under all error structures. Selecting $p=\infty$ gives the best performance in all columns. Error correlations are not adversarial enough for $p=2$ to dominate the other designs, as it did in the previous application,\footnote{A possible reason is that here $x_{i}$ is the number (rather than share) of treated neighbors, so there is more scope for increasing signal by correlating the assignments---which is what the $p=\infty$ design targets.} but it remains a good choice, especially with more complex error structures. Overall, the power gains are even more substantial than in the first application: with the $p=2$ choice, they range from +76% in column 2 to +241% in column 3, while with $p=\infty$ they reach +570% in column 3.\footnote{To isolate the performance impact of the optimal design vs. the optimal estimator, Appendix Table (ref) reports RMSE of the recentered IV estimator without whitening for the $p$-optimal designs.} Yet, Appendix Table (ref) provides a reason to prefer a conservative choice of $p$: under $p=\infty$, the resulting covariance matrix of the shocks is too dense to meet the conditions of Theorem (ref). For the simulations with correlated errors, confidence intervals built using the normal approximation under the $p=\infty$ design severely undercover the target parameter, while the $p=1$ design has close to nominal coverage under all DGPs.

Since our proposed designs are optimized for the indirect effect, they can lead to some deterioration of performance of the direct effect estimator, as reported in Panel (b) of Table (ref). For $p=\infty$ this deterioration can be substantial, equivalent to a 59% reduction in sample size in column 6. However, for $p=1$ and $p=2$ the deterioration is mild, of at most 22% and 27%, respectively.

Conclusion

We have shown how researchers can design experiments and choose recentered estimators to more precisely estimate spillover effects, under weak restrictions on the distribution of unobservables. The core insight of this approach is that spillover signal can be increased, by correlating the assignments of connected units, without overly concentrating variation in directions where the unobservables may be adversarially clustered. Balancing this tradeoff between signal and isotropy yields a tractable minimax problem with intuitive special cases and practical approaches to computation. In two semi-synthetic experiments, based on high-profile applications, we find large decreases in spillover effect standard errors relative to conventional randomization and estimation that are possible without large deterioration in the standard errors for the direct effects in the application with such effects.

Several open questions remain. First, while we have focused on linear spillover treatments, and derived theoretical extensions to nonlinear spillover treatment constructions, computation of optimal designs in nonlinear settings remains a challenge. Second, though we have characterized the optimal design and estimator for spillovers in the presence of direct effects, we have not presented the more general problem of minimizing a weighted average of their variances. Third, we have not shown how the optimal design problem changes when a researcher has access to some prior information on the distribution of unobservables, such as from a first-wave pilot (viviano_experimental_2025). Finally, we have not studied how the minimax approach changes when a researcher is interested in estimating a more complicated parametric model of spillovers, such as from a structural economic model (borusyak_estimating_2025). We hope to address some of these extensions in future drafts.