Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
98,499 characters · 29 sections · 42 citation commands
Coupling Designs for Randomized Experiments with Complex Treatments
\begingroup \endgroup \setcounter{footnote}{0}
Consider assigning treatments $D_i$ from a distribution $F$ to units $i \in [n]$ in a randomized experiment. It is widely recognized that for simple discrete distributions $F$, experimenters can improve estimation efficiency by stratifying on covariates at design time. For example, if $D_i \in \{0, 1\}$ and $F = \operatorname{Bernoulli}(1/2)$, a matched pairs design would match units with similar baseline covariates $X_i \approx X_j$ into pairs, then set $D_i=1$ and $D_j=0$ or vice-versa, with equal probability. Doing so balances covariates between treatment and control groups, improving precision for generic parametric causal estimators bruhn2009,imai2009,bai2023efficiency. Similarly, if there are $k$ treatments and $F = \operatorname{Unif}([k])$, one can improve efficiency by matched k-tuples randomization: match units into groups of $k$ with similar baseline covariates, then randomly assign one unit to each treatment cochran1957experimental.
Moving beyond these well understood situations, consider a researcher interested in the effect of cash grants $d \in [0, u]$ on future household consumption $Y_i(d)$. They want to estimate an approximation of the average dose-response curve $d \mapsto n^{-1} \sum_i Y_i(d)$ for grant amounts $d \in [0, u]$. Identification of the full curve requires continuous randomization of the treatment, such as $F = \operatorname{Unif}[0, u]$. However, this makes stratification impossible since there are infinitely many treatment levels. One possible solution is to just assign treatments independently between units, forgoing any efficiency improvements.
An alternative solution that potentially improves efficiency is to first discretize the treatment, say into $k=20$ treatment levels, then apply matched $k$-tuples randomization. However, such discretization generally requires changing the causal estimand, since the dose-response is non-identified at intermediate points. Even leaving such identification issues aside, match quality deteriorates rapidly as $k$ increases, reducing the efficiency gains from stratified randomization. Indeed, in moderately sized experiments, it can be challenging to find well-matched pairs of $k=2$ units even for relatively low-dimensional covariates, and will be much more difficult to find well-matched groups of $k=20$ units.
The challenge of improving efficiency in experiments with complex treatments is not limited to continuous treatments. The purpose of this paper is to develop a new family of coupling designs that extends the basic principle of stratification to allow for efficient randomization in experiments with continuous, constrained multivariate, text/image and other irregular treatment spaces. We make three main contributions in this paper:
Experimental design for causal inference problems goes back to at least fisher1926, and has been an active area of research in statistics and econometrics for decades. Recent surveys include athey2017survey and bai2025primer. Several approaches to improving efficiency in randomized experiments have been developed in this literature, including rerandomization morgan2012,li2018rerandomization, pure optimization kasy2016,kallus2017balance, as well as approaches based on discrepancy minimization harshaw2024. However, these approaches are not designed to handle complex treatment spaces. For example, none of these approaches can accommodate continuous treatment distributions.
There is a large literature on improving efficiency by stratified randomization, which is closely related to our work. For example, matched pairs designs assign opposite treatments to pairs of similar units fisher1935,Greevy2004Optimal,bruhn2009. The efficiency properties of matched pairs have been studied by imai2008, fogarty2018, bai2020pairs, bai2021inference, pashley2021, among others. Generalizations to matched $k$-tuples and other stratified designs have been considered by several authors cochran1957experimental,higgins2016blocking,bugni2018inference,bugni2019,cytrynbaum2023,bai2023tuples,bai2023efficiency. Perhaps the closest paper to our work is koo2026incomplete, which provides a modern, potential outcomes based analysis of incomplete block designs yates1936incomplete,kempthorne1956efficiency, which classically were studied using a restrictive, model-based approach. Incomplete block designs refer to stratified randomization where the total number of distinct treatments is larger than the block (stratum) size. Similarly, in our setting it will either be highly inefficient or outright impossible to implement every distinct treatment within each matched group.
The coupling designs we introduce in this paper extend the basic mechanism of stratified randomization to complex treatment spaces by replacing treatment permutations within strata with general negatively dependent couplings that disperse treatment assignments across the treatment space $\mathcal D \subseteq \mathbb{R}^m$. In addition to enabling covariate-balancing randomization in a much broader class of useful experiments, coupling designs also provide new insights into conventional stratified designs, and can even improve on such designs in classical settings by more efficiently trading off between the key forces of dispersion and match quality.
To generate highly dispersed treatments within matched groups, we combine coupling techniques from the Monte Carlo integration literature robert2004monte,Owen2013 with geometry-preserving maps from optimal transport theory brenier1991polar,merigot2011multiscale,carlier2010knothe.
Coupling designs extend the basic mechanism of stratification to allow efficient randomization from any distribution $F$ within tightly matched groups of units. To illustrate this idea, consider extending matched pairs designs to allow randomization of continuous, univariate treatments $D_i \in \mathbb{R}$. The conventional matched pairs design can be understood as drawing $(D_i, D_j) \sim G$ from a coupling $G$ with fixed marginals $G_i = G_j = \operatorname{Bernoulli}(1/2)$ and $D_i = 1-D_j$, which achieves maximal negative correlation: $\operatorname{Corr}_G(D_i, D_j) = -1$. For more general distributions $F$, this perspective suggests first matching pairs of similar units $i$ and $j$, then drawing $(D_i, D_j) \sim G$ from a coupling with fixed marginals $G_i = G_j = F$ and strong pairwise negative correlation: $\operatorname{Corr}_G(D_i, D_j) \ll 0$. One way to construct such a coupling for any $F$ is via a classic idea from Monte Carlo integration theory known as antithetic variates sampling hammersley1956antithetic.
\paragraph{Antithetic Matched Pairs.}
Antithetic variates sampling generates assignments $D_i^* = F^{-1}(U)$ and $D_j^* = F^{-1}(1-U)$ using a common uniform variate $U \sim \operatorname{Unif}[0, 1]$. Since both $U$ and $1-U$ are uniform, the quantile transform produces $D_i^*, D_j^* \sim F$ marginally. By drawing opposing quantiles $U$ and $1-U$, the coupling $(D_i^*, D_j^*) \sim G$ induces strong negative correlation between the treatments. In fact, the results of hoeffding1940masstab show that for any monotone function $y(\cdot)$ of the treatment, antithetic variates achieves minimal correlation:
We can construct an antithetic matched pairs design by assigning treatments $(D_i^*, D_j^*) \sim G$ within tightly matched pairs, inducing strong negative correlation while also implementing the chosen marginal $D_i^* \sim F$ for each unit $i \in [n]$. When $F = \operatorname{Bernoulli}(1/2)$, the design yields $D_i^* = 1 - D_j^*$, recovering the conventional matched pairs design. However, this construction can be used more generally for any univariate distribution $F$, for example $F = \operatorname{Unif}[0, u]$.
Antithetic variates were developed to improve efficiency in Monte Carlo integration problems, such as estimating $\theta_0 = E_F[y(D)]$ using $n^{-1} \sum_i y(D_i)$. In contrast to classic Monte Carlo integration, in causal inference units generally have heterogeneous responses to the treatment with $Y_i(\cdot) \neq Y_j(\cdot)$. We can use matching to enforce approximate homogeneity at the group level. After successful matching, $Y_i(\cdot) \approx Y_j(\cdot)$ within matched pairs, allowing the design to leverage the efficiency improvements from antithetic variates as if the responses were homogeneous.
Simplifying for illustration, consider an estimator $\widehat{\theta} = (1/2)(Y_i(D_i) + Y_j(D_j))$ of the average outcome $\theta_0 = (1/2)(E_F[Y_i(D)] + E_F[Y_j(D)])$. If units are perfectly matched $Y_i(\cdot) = Y_j(\cdot) = y(\cdot)$, the variance relative to independent assignment is
If $y(\cdot)$ is monotone, then by Equation (ref) this antithetic pairs design can significantly improve efficiency. Indeed, it minimizes variance of the estimator $\widehat{\theta}$ among all pair-wise couplings $G$.
The simple example above illustrates how matched pairs randomization can be extended to randomize efficiently from any univariate distribution $F$ when the response function is monotone. In what follows, we generalize this construction to provide coupling-based randomization methods for general treatment spaces $\mathcal D \subseteq \mathbb{R}^m$. We show that these methods improve efficiency under weak smoothness conditions on the potential outcomes $Y_i(\cdot)$.
We construct general coupling designs by first matching the experimental units into homogeneous groups of size $k \ge 2$ using covariates. Next, treatments are drawn within each group from a coupling $(D_i)_{i=1}^k \sim G$ with marginals $G_i = F$ for all $i \in [k]$, for a fixed distribution $F$ over the treatment space $\mathcal D = \operatorname{Supp}(F)$. We can view this as a matched $k$-tuples design with coupling-based randomization within each group. The family of coupling designs can accommodate very general marginal distributions $F$ and spaces $\mathcal D$. When $F$ is discrete, we also recover conventional stratified randomization for appropriate choices of $k$ and $G$. There are many possible couplings that can be used to randomize within groups. We discuss several examples and provide a general construction in Section (ref) below.
Dispersed Treatments. As in Equation (ref), we would like to produce negatively correlated treatments $(D_i)_{i=1}^k$ within matched groups of $k$. For multivariate $D \in \mathbb{R}^m$, this can be achieved by sampling $(D_i)_{i=1}^k$ to be highly dispersed or “spread out” over the treatment space $\mathcal D$. Making treatments dispersed in this way implies that for sufficiently smooth functions $\phi: \mathcal D \to \mathbb{R}$,
When tuple size $k$ is large, it is possible to achieve higher dispersion, since this allows us to coordinate randomization of $(D_i)_{i=1}^k$ to cover more of the treatment space. However, as tuple size $k$ increases, match quality becomes worse and the response functions $Y_i(\cdot)$ and $Y_j(\cdot)$ of matched units $i$ and $j$ are less similar on average. Using a formalization of these concepts introduced in Section (ref), we show that the efficiency gain from a coupling design relative to independent assignment is
Intuitively, by assigning matched $k$-tuples of similar units to highly dissimilar treatments $(D_i)_{i=1}^k$, we prevent spurious in-sample correlations from arising between the treatment assignments and unit-specific heterogeneity. We formalize how this generalizes classical notions of covariate balance to much more complex treatment spaces and distributions.
To show our main result in full generality, Section (ref) defines a coupling operator $U_G$ whose eigenspaces in $L^2(F)$ can be viewed as the principal directions of the coupling $G$ with respect to random sampling. The dispersion $\operatorname{Disp}_G(\phi)$ of any $\phi \in L^2(F)$ decomposes orthogonally over these eigenspaces, so the overall efficiency gain depends on how well the influence functions $s_i(\cdot)$ align with the high-dispersion directions of $G$. This yields a general decomposition of the variance reduction from coupling designs into a weighted sum of eigenspace-specific dispersion $\times$ match quality terms, where the weights reflect the approximation quality of $s_i(\cdot)$ on each eigenspace.
To make the ideas of the paper concrete, we describe an application of an experiment with complex treatments where a coupling design can be used but conventional stratification will either be impossible or have poor performance.
Consider estimating the probability that a product $d$ is purchased, say $Y_i(d) = 1$, when it is shown to unit $i$. For example, the product could be a restaurant on a food delivery platform that provides individual restaurant promotions to its users. While all products are unique, they can typically be compared. For example, we could featurize each restaurant $d_r \in \mathbb{R}^m$ using a vector of attributes like cuisine type, average price, rating, and so on. Then the treatment space $\mathcal D = \{d_1, \dots, d_R\}$ is an irregular, discrete subset of $\mathbb{R}^m$. Suppose the platform uses a discrete choice model to estimate user preferences. They assign $D_i \sim F$ over restaurants $\mathcal D$, observe choices $Y_i(D_i)$, then estimate a Logit discrete choice model with
In Section (ref) below, we show the Logit MLE consistently estimates an approximation to the dose-response $n^{-1} \sum_i Y_i(d)$, the proportion of units who would choose to purchase when shown product $d$. We can use coupling designs to improve precision for estimating such causal parameters. To do so, we construct couplings $G$ that produce high dispersion over the irregular point cloud of possible product types, then sample $(D_i)_{i=1}^k \sim G$ within matched $k$-tuples of similar users. Our construction combines tools from Monte Carlo integration theory and optimal transport (OT), see Section (ref) for details and Figure (ref) for a visualization of the treatments $(D_i)_{i=1}^k$.
For comparison, stratified randomization would randomly assign each of $R$ restaurants without replacement to groups of $k = R$ users, uniformly at random. If the experiment size $n < R$, this is not possible. More generally, if $R$ is large, stratified randomization will behave similarly to iid randomization, due to poor match quality within groups of $k \ge R$ units. By contrast, coupling designs allow experimenters to draw treatments $(D_i)_{i=1}^k$ within matched $k$-tuples for any $k \ge 2$, preserving match quality while also dispersing treatments through the space $\mathcal D$.
Exploiting Smoothness. This example also shows how coupling designs can naturally exploit the smoothness of the outcome functions $Y_i(d)$ over $d \in \mathcal D$. To see this, note that it is natural to assume that units have similar preferences $Y_i(d_r) \approx Y_i(d_l)$ over restaurants with $|d_r - d_{l}| \approx 0$, a form of smoothness. Because of this, assigning similar pairs of units $Y_i(\cdot) \approx Y_j(\cdot)$ to restaurants $d_r$ and $d_{l}$ effectively wastes a sample: we learn the same thing from observing $Y_i(d_r)$ and $Y_j(d_l)$. Coupling designs prevent this by assigning matched groups of similar units $Y_i(\cdot) \approx Y_j(\cdot)$ to highly dissimilar restaurants. Our core theory in Sections (ref) and (ref) explicitly connects the efficiency gain from coupling designs to the smoothness of the responses $Y_i(\cdot)$.
We consider estimators $\widehat{\theta} = n^{-1} \sum_i s_i(D_i)$, where the functions $s_i(\cdot)$ are non-random in a design-based framework. For example, we could have $s_i(d) = s(Y_i(d), d, X_i)$ for a fixed function $s(\cdot)$ of the data. The corresponding finite-population estimand is
We can view $\theta_n$ as a “fully heterogeneous” version of the classic Monte Carlo estimand $\theta_0 = E_F[s(D)]$. In fact, these types of estimators and estimands are ubiquitous in causal inference problems harshaw2025general. Thus, experimental design in causal inference can be viewed as a heterogeneous Monte Carlo integration problem. For brevity, in what follows we will denote $E_n[a_i] \equiv n^{-1} \sum_i a_i$ for any array $(a_i)_{i=1}^n$.
\paragraph{Influence Functions.}
A variety of estimators admit the design-based asymptotic linearization $\widehat \beta - \beta_n = E_n[s_i(D_i)] + o_p(n^{-1/2})$. Because of this, our efficiency analysis of the quantity $E_n[s_i(D_i)]$ under coupling designs also characterizes the first-order efficiency of significantly more general parametric estimators. For example, let $\widehat \beta$ be the coefficient from the OLS regression $Y_i \sim 1 + D_i$. Then, under weak conditions on the design,
We have $s_i(d) = e_i(d)H(d)$ for residual $e_i(d) = Y_i(d) - E_F[\bar Y_n(D)] - \theta_\text{BLP}'(d - E_F[D])$. We provide more details in Appendix (ref). In a slight abuse of terminology, we refer to $s_i(\cdot)$ as unit $i$'s influence function.
Examples (ref) and (ref) can be extended to accommodate regressions $Y_i \sim 1 + t(D_i)$ or logit with a set of basis functions $t(d)$. For example, if $D \in \mathbb{R}^2$, we could choose $t(d) = (d_1, d_2, d_1 d_2)$ for the regression $Y_i \sim 1 + D_{i1} + D_{i2} + D_{i1}D_{i2}$, as in Example (ref) above. Similar to the previous results, this provides the best approximation of $\bar Y_n(d)$ among all linear functions with both main effects and two-way interactions.
A coupling design is implemented in three steps:
The first step of matching is common to all stratified designs, while the second and third steps are unique to the coupling designs we introduce in this paper. Our theory is largely agnostic about the specific details of the matching algorithm. The ultimate aim is to find homogeneous groups of units with similar response functions $Y_i(\cdot)$. One possible proxy for this objective is to minimize a covariate discrepancy:
There are several algorithms for this problem Greevy2004Optimal,bai2021inference,cytrynbaum2023. Due to the curse of dimensionality in matching, this should be done using a small set of covariates expected to be highly predictive for endline outcomes $Y_i$.
A coupling is a joint distribution $G$ over $\mathcal D^k$ with fixed marginals. For tractability and expositional clarity, we consider couplings $G$ that are exchangeable and have identical marginals $G_i = F$ for $i \in [k]$. A coupling is exchangeable if the joint distribution of $(D_i)_{i=1}^k$ is invariant to permutations of the indices $i \in [k]$. We define the set of feasible couplings $\Pi_k(F)$ to be the exchangeable joint distributions $G$ over $\mathcal D^k$ with fixed marginals $G_i = F$.
There are many strategies for sampling highly dispersed uniform random variables $(U_i)_{i=1}^k \sim G_U$. Inspired by the Monte Carlo integration literature, we consider three canonical examples based on the Gaussian copula, stratification, and randomly shifted lattices.
Latin hypercube produces samples $(U_i)_{i=1}^k \sim G_U$ such that the coordinate projections $(U_{ij})_{i=1}^k$ are spread out over the interval $[0, 1]$ for each dimension $j \in [m]$. However, even if the one-dimensional projections $(U_{ij})_{i=1}^k$ are highly dispersed, the samples $(U_i)_{i=1}^k$ may not be jointly well-dispersed through the hypercube $[0, 1]^m$ if $m > 1$. Because of this, Latin hypercube produces strong negative correlations $\operatorname{Corr}_G(\phi(D_i), \phi(D_j)) \ll 0$ only for univariate functions like $\phi(d_1, \dots, d_m) = d_1$, but tends to have weaker effects for jointly varying functions like $\phi(d) = d_1 d_2$.
This problem can be solved with more advanced coupling constructions. Some examples include orthogonal array Latin hypercube sampling owen1992orthogonal, tang1993orthogonal, scrambled digital nets owen1995scrambled, and shifted rank-1 lattice rules cranley1976randomization, sloan1994lattice. For brevity, we only formally describe shifted lattice rules. In the univariate case, such couplings draw $(U_i)_{i=1}^k \sim G_U$ by adding a uniform random shift to a regular grid $(l/k)_{l=0}^{k-1} \subseteq [0, 1]$. For more general $m \ge 1$, they apply a random shift $S \sim \operatorname{Unif}[0, 1]^m$ to a dispersed lattice of points determined by a number-theoretic construction.
For $m=1$ and $z=1$, we have $(U_i)_{i=1}^k = (\pi(i)/k + S)_{i=1}^k \pmod 1$, a random shift and permutation of the regular grid $(l/k)_{l=0}^{k-1}$, which is known as a rotation sampling fishman1983antithetic. For $m \ge 2$, the condition $\gcd(z_j, k) = 1$ implies that the projections $(U_{ij})_{i=1}^k$ onto each coordinate $j \in [m]$ are themselves rotation samples. In addition to being marginally dispersed on each dimension, a well-chosen generating vector $z$ ensures that the $(U_i)_{i=1}^k$ are also jointly well-dispersed through $[0, 1]^m$. For a theoretical analysis of the choice of generating vector $z$, see kuo2003component. We slightly modify the standard construction of sloan1994lattice, adding a random permutation $\pi$ for exchangeability due to heterogeneity of the units in causal inference problems.
In general, we must map the uniform samples $(U_i)_{i=1}^k$ to treatments $D_i = T(U_i)$ that are both dispersed over $\mathcal D$ and have the correct marginal distribution $F$.
For univariate treatments $m=1$, we can use the transport map $T(u) = F^{-1}(u)$, setting $D_i = F^{-1}(U_i)$, so that $D_i \sim F$ by the properties of the quantile function. In the multivariate case, if $F = \otimes_{j=1}^m F_j$ has independent components $D_{ij} \perp \!\!\! \perp D_{il}$ for $j \neq l \in [m]$, then we can similarly enforce $T(U_i) \sim F$ by applying the quantile transform componentwise:
However, this will not work for general multivariate distributions $F$ with dependent components, since $D_i = T(U_i)$ would not have the correct joint distribution. Thus, the componentwise quantile transform can only be used when the treatment space is a product set $\mathcal D = \times_{j=1}^m \mathcal D_j$ and the desired marginal distribution $F$ has independent components. Note that several of the applications described in Section (ref) do not satisfy this structure.
The quantile transform preserves the geometry of the samples $(U_i)_{i=1}^k$ in the sense that if $U_i$ and $U_j$ are far apart in $[0, 1]$, then $F^{-1}(U_i)$ and $F^{-1}(U_j)$ will also be far apart in $\mathcal D$. To accommodate complex treatment spaces and dependent components, we must construct general geometry-preserving maps $T: [0, 1]^m \to \mathcal D$ with $T(U) \sim F$. The geometry-preserving condition naturally leads us to Brenier maps from optimal transport brenier1991polar:
In the univariate case and when $F = \otimes_{j=1}^m F_j$, the Brenier map recovers the quantile transform in Equation (ref) above. However, optimal transport can be used to construct geometry-preserving maps $T^*(U) \sim F$ for much more general spaces $\mathcal D$. For discrete spaces $\mathcal D \subseteq \mathbb{R}^m$, such maps can also be computed efficiently by semi-discrete optimal transport merigot2011multiscale. We discuss the geometric condition motivating this definition and further computational details in Appendix (ref).
This section describes the key mathematical objects for quantifying the relative efficiency of coupling design randomization, defining appropriate measures of match quality and sample dispersion. This allows us to show a simple fundamental relation between these objects and the efficiency gain from a coupling design:
Above, we constructed couplings that “spread out” treatments over the treatment space $\mathcal D$. For intuition about how this improves precision, consider the simple problem of estimating $\theta_0 = E_F[\phi(D)]$ with $\widehat{\theta} = k^{-1} \sum_{i=1}^k \phi(D_i)$. This is a homogeneous version of the more general heterogeneous estimation problems in Section (ref).
If by random chance we sample treatments $D_i \approx D_j$ close together in the space $\mathcal D$, then for smooth enough functions $\phi: \mathcal D \to \mathbb{R}$, the samples $\phi(D_i) \approx \phi(D_j)$ will be quite similar, which effectively wastes an experimental sample. By contrast, if $(D_i)_{i=1}^k$ are dispersed over $\mathcal D$, then we learn more about the function $\phi(\cdot)$ from a given fixed sample size $k$.
To formalize this intuition, let sample variance $\operatorname{Var}_k(a_i) \equiv (k-1)^{-1} \sum_{i=1}^k (a_i - \bar a)^2$ for any $(a_i)_{i=1}^k$ and define the sample dispersion as a normalized measure of how spread out the samples $(\phi(D_i))_{i=1}^k$ are in expectation over $G$.
For the iid design $G_{iid} = \otimes_{i=1}^k F$, we have $E_{G_{iid}} \operatorname{Var}_k(\phi(D_i)) = \operatorname{Var}_F(\phi)$ by unbiasedness of the sample variance, so $\operatorname{Disp}_{G_{iid}}(\phi) = 0$. If $\operatorname{Disp}_G(\phi) > 0$, then the samples $(\phi(D_i))_{i=1}^k$ are more spread out in expectation under $G$ than under iid randomization. For homogeneous estimation problems and smooth $\phi(\cdot)$, this improves efficiency by the mechanism described above. Indeed, the relative efficiency for the homogeneous problem above is
Because $\operatorname{Var}_G(\widehat{\theta}) \geq 0$, we have $\operatorname{Disp}_G(\phi) \leq 1$ for all $\phi(\cdot)$ and $G \in \Pi_k(F)$. We also have the lower bound $\operatorname{Disp}_G(\phi) \ge -(k-1)$, which is attained for any $\phi(\cdot)$ under a clustered coupling with $D_i = D_j$ for all $i, j \in [k]$. Under exchangeable couplings, the dispersion can be interpreted as a normalized measure of negative correlation between the samples $\phi(D_i)$ and $\phi(D_j)$ for $i \neq j$.
Negatively correlated samples tend to “repel” each other, making them more spread out. We work primarily with this formulation in what follows.
\paragraph{Role of Smoothness.}
Recall the discrete choice example in Section (ref), where we argued that coupling designs can exploit smoothness of the purchase decision $y(d)$ in restaurant features $d \in \mathbb{R}^m$. Let $n = k$ and suppose preferences are homogeneous with $Y_i(d) = y(d)$. We can use the Horvitz-Thompson estimator to estimate the best linear approximation of $y(\cdot)$ as in Example (ref). Then $\widehat{\theta} = k^{-1} \sum_{i=1}^k \phi(D_i)$ with influence function $\phi(d) = y(d)H(d)$ for weights $H(d) = \operatorname{Var}_F(D)^{-1}(d - E_F[D])$. When $y(\cdot)$ is smooth, so is $\phi(\cdot)$, and the samples $\phi(D_i) \approx \phi(D_j)$ whenever $D_i \approx D_j$. By spreading $(D_i)_{i=1}^k$ across the treatment space $\mathcal D$, the couplings above produce more dispersed samples $(\phi(D_i))_{i=1}^k$, increasing $\operatorname{Disp}_G(\phi)$ and improving efficiency (Equation (ref)). In Section (ref), we develop the technical machinery to describe exactly how $\operatorname{Disp}_G(\phi)$ is determined by the smoothness and shape of $\phi(\cdot)$.
The direct connection between dispersion and efficiency in Equation (ref) above only holds for homogeneous problems with $\widehat{\theta} = k^{-1} \sum_{i=1}^k \phi(D_i)$. For realistic causal estimation problems, we also need to account for heterogeneity of the functions $s_i(\cdot)$ within matched $k$-tuples of units, due to imperfect matching on only partially predictive covariates. To do so, next we define an appropriate measure of within-group match quality.
The term $v_\Delta(s)$ is a design-based matching discrepancy, with $v_\Delta(s) = 0$ under perfect matching, $s_{ig} = s_{jg}$. The match coefficient $Q_k(s)$ measures how homogeneous the functions $s_{ig}(\cdot)$ are within each matched $k$-tuple, with $Q_k(s) = 1$ under perfect matching. More generally, we have lower and upper bounds $-(k-1)^{-1} \le Q_k(s) \le 1$ for all populations $s = (s_i)_{i=1}^n$ and matching procedures. If units are matched at random, the expected match quality is $E_{\tau}[Q_k(s)] \ge -(n-1)^{-1}$ in the worst case. In theory, it is possible to approach the lower bound by purposefully matching units into $k$-tuples to be as dissimilar as possible, but we do not expect this to arise in practice. See Appendix (ref) for further details.
\paragraph{Covariate Power.}
Recall units are matched into $k$-tuples using observed baseline covariates $(X_i)_{i=1}^n$. The match quality coefficient $Q_k(s)$ will be large if both:
We formalize this observation in Appendix (ref), providing conditions under which $Q_k(s) = R^2_{s|X} \cdot Q_k(\mu) + o_p(1)$ as $n \to \infty$. Here $Q_k(\mu)$ is the matching discrepancy on features of the covariates alone, and $R^2_{s|X}$ measures covariate predictive power for heterogeneity in $s_i(\cdot)$.
Our main theoretical result is that the efficiency gain from coupling designs is proportional to the product of dispersion $\operatorname{Disp}_G(\phi)$ and match quality $Q_k(s)$. To state this result in full generality requires additional technical machinery, which we develop in Section (ref) below. To build intuition, we first state the result in a simple univariate parametric model.
For intuition, note that in the homogeneous special case $s_i(d) = c + a \phi(d)$ for all units $i \in [n]$, relative efficiency is exactly given by $\operatorname{Disp}_G(\phi)$. In general heterogeneous problems, the relative efficiency is dampened by imperfect matching, captured by the match quality coefficient $Q_k(s) < 1$.
Another interesting special case occurs when $\operatorname{Disp}_G(\phi) = 1$, so that relative efficiency is exactly equal to match quality $Q_k(s) = 1 - v_\Delta(s)/v_{iid}(s)$. Rearranging, we find $n\operatorname{Var}_G(\widehat{\theta}) = v_\Delta(s)$, the matching discrepancy from Definition (ref). This shows how $v_\Delta(s)$ can be interpreted as the ideal variance under perfect dispersion, where efficiency is only limited by imperfect matching.
Tuple Size Tradeoff. The quantity $\operatorname{Disp}_G(\phi)$ is generally increasing in tuple size $k$, since for larger $k$ it becomes easier to jointly correlate the treatments $(D_i)_{i=1}^k$ to be spread out over $\mathcal D$. Our analysis in Sections (ref) and (ref) formalizes this effect. By contrast, match quality $Q_k(s)$ is generally decreasing in $k$, since it becomes harder to find many similar units to match together.
Consider the two extremes of this trade-off. At one extreme, we can set $k=2$, sampling $(D_1, D_2) \sim G$ within matched pairs. This makes matches as tight as possible, maximizing $Q_k(s)$, then uses the coupling $G$ to increase dispersion subject to the constraint $k=2$.
At the opposite extreme, we can maximize dispersion by setting $k$ very large and sampling $(D_1, \dots, D_k) \sim G$ to be as dispersed as possible over $\mathcal D$. However, match quality will suffer due to the large tuple size. For general $F$, perfect dispersion is only possible in the limit as $k \to \infty$, but $\operatorname{Disp}_G(\phi) = 1$ is possible for some special distributions $F$ with small support.
\paragraph{Illustration of Tradeoff.}
Consider a researcher estimating the dose-response of welfare outcomes to cash grants, randomizing from $F = \text{Exp}(1)$ to provide many small grants while also trying some larger amounts. Suppose the treatment has no effect on potential outcomes, so $Y_i(d) = \mathds{1}'X_i$ and $s_i(d) = (\mathds{1}'X_i) H(d)$. Theorem (ref) shows that the efficiency under a coupling $G$ is the product of $\operatorname{Disp}_G(H)$ and match quality $Q_k(s)$. The feasible match quality vs.\ dispersion frontier and efficiency gain are shown in Figure (ref) as functions of $k$. Perfect treatment dispersion is impossible here, requiring $k \to \infty$. However, efficiency is actually maximized at moderate $k$, e.g.\ around $k=4$, accepting reduced dispersion in exchange for efficiency gains from higher match quality.
An important motivation for conventional stratified randomization is that it reduces covariate imbalances between the different treatment groups. These imbalances are equivalent to spurious in-sample correlations between the treatment $D_i$ and covariates $X_i$. Here, we show that by assigning similar groups of units to highly dispersed treatments, coupling designs prevent such spurious correlations from arising. Because of this, such designs enable covariate-balancing randomization over complex treatment spaces.
For binary $D_i \in \{0, 1\}$, covariate balance is commonly assessed using the t-statistics from a regression $X_i \sim 1 + D_i$. Up to a normalization, this is equivalent to checking the magnitude of the sample covariance $\operatorname{Cov}_n(D_i, X_i)$. If this covariance is small, then treatments are approximately independent of unit-specific heterogeneity in-sample. Thus, ensuring covariate balance is equivalent to randomizing in a way that enforces $E_G[\operatorname{Cov}_n(D_i, X_i)^2] \approx 0$.
To extend this balance measure to complex treatment spaces $\mathcal D$, let $\phi: \mathcal D \to \mathbb{R}$ and $b: \mathbb{R}^p \to \mathbb{R}$ be basis functions. Write $\operatorname{Var}_n(b_i) = n^{-1} \sum_{i=1}^n (b_i - \bar b)^2$ for $b_i = b(X_i)$ and define the imbalance under a coupling $G$ as the mean-squared sample covariance
We can view $D_i$ as approximately independent of covariates $X_i$ if $\mathcal I_G(\phi, b) \approx 0$ for a rich set of basis functions $\phi(\cdot)$ and $b(\cdot)$. Intuitively, coupling designs prevent such in-sample correlations from arising between $X_i$ and $D_i$ by ensuring that similar units are assigned to highly dissimilar treatments. In particular, the next corollary shows that the imbalance measure $\mathcal I_G(\phi, b)$ is decreasing in the product $\operatorname{Disp}_G(\phi) \times Q_k(b)$ of dispersion and match quality.
This section proves our efficiency result in full generality, without the simplifying parametric assumption we imposed when previewing the results in Section (ref). To do so, we decompose the space of influence functions into a basis of orthogonal subspaces on which $\operatorname{Disp}_G(\cdot)$ is constant, which we view as the principal directions of the coupling $G$.
Define the square-integrable functions $L^2(F) = \{\phi: E_F[\phi(D)^2] < \infty\}$. We assume $s_i(\cdot) \in L^2(F)$ for $i \in [n]$ throughout. Also denote the mean zero subspace $L_0^2(F) \equiv \{ \phi \in L^2(F) : E_F[\phi(D)] = 0 \}$. We define the principal directions of a coupling $G$ to be the eigenspaces of the following linear operator, which captures the pairwise dependence structure of the design.
This operator is well-defined for any choice of $i \ne j$ due to exchangeability of $G$. Since dispersion $\operatorname{Disp}_G(\phi)$ is invariant to constant shifts, it is convenient to work with the mean zero subspace $L_0^2(F)$ in what follows. Note that if $\phi \in L_0^2(F)$, then by tower law $E_F[(U_G \phi)(D)] = 0$, so $U_G$ also maps $L_0^2(F)$ to itself. The coupling operator is self-adjoint and linear on the Hilbert space $L_0^2(F)$, so it has a real spectrum. We additionally require the following condition:
Assumption (ref) holds for any coupling $G \in \Pi_k(F)$ if the marginal $F$ is discrete. For continuous $F$, Lemma (ref) in the appendix shows that it holds for all of the univariate couplings described in Section (ref). More generally, by the spectral theorem Assumption (ref) holds whenever the operator $U_G$ is compact.
A key insight for our analysis is that $\operatorname{Disp}_G(\phi)$ decomposes orthogonally over the eigenspaces of $U_G$ for any $\phi \in L^2(F)$.
We view the eigenspaces of $U_G$ as the principal directions of the coupling $G$ in $L^2(F)$ with respect to sampling, since they control how much dispersion is produced when sampling from $G$. In particular, $\operatorname{Disp}_G(\phi)$ will be large if $\phi(\cdot)$ is well approximated on the high dispersion eigenspaces, e.g.\ with $\operatorname{Disp}_G(E_m) \approx 1$. The following corollary is immediate.
More generally, if $L_0^2(F) = \oplus_{m=1}^M E_m$, then each eigenspace can be obtained by maximizing $\operatorname{Disp}_G(\phi)$ subject to orthogonality to previously found eigenspaces. This is analogous to principal components analysis (PCA), where each principal component can be found by maximizing data variance subject to orthogonality to previous components.
Coupling Analysis. Theorem (ref) also provides a simple recipe for computing $\operatorname{Disp}_G(\phi)$ by analyzing the eigenspaces and eigenvalues of the operator $U_G$. We illustrate this by computing the exact dispersion for the univariate Latin hypercube coupling. Our analysis highlights the role played by smoothness of the function $\phi(\cdot)$ in guaranteeing high dispersion $\operatorname{Disp}_G(\phi)$ under this coupling. We provide a detailed comparison with other couplings in Section (ref) below.
The dispersion basis in Theorem (ref) allows us to generalize the simple relationship between efficiency and the product of dispersion $\times$ match quality from Theorem (ref) to general influence functions $s_i(\cdot) \in L^2(F)$. Impose Assumption (ref) and let $s_i^m(d) = (P_m s_i)(d)$ be the projection of $s_i(\cdot)$ onto $E_m$. Define approximation weights
The weights $w_m(s)$ quantify how well the influence functions can be approximated using functions in eigenspace $E_m$, on average over the experimental units. In what follows, denote $s = (s_i)_{i=1}^n$ and $s^m = (s_i^m)_{i=1}^n$. We also write $\operatorname{Disp}_G(m)$ as the common dispersion on $E_m$.
The weights $w_m(s)$ are non-negative and sum to one. Thus, the theorem shows that the efficiency gain from coupling design randomization is a convex combination of products of dispersion $\operatorname{Disp}_G(m)$ and match quality $Q_k(s^m)$ across eigenspaces $E_m$. The efficiency gain is therefore large when the influence functions $s_i(\cdot)$ align well with eigenspaces with high dispersion $\operatorname{Disp}_G(m) \approx 1$ and good match quality $Q_k(s^m) \approx 1$. We show that many of the couplings introduced above achieve high dispersion generically for smooth functions $s_i(\cdot)$.
\paragraph{Approximate Stratification.}
Recall the definition of the match quality coefficient $Q_k(s) = 1 - v_\Delta(s) / v_{iid}(s)$ from Definition (ref), where $v_{iid}(s)$ is the iid variance and $v_\Delta(s)$ is the average within-group variance of the influence functions. When units are well-matched, we have $Q_k(s) \approx 1$ so that $v_\Delta(s) \ll v_{iid}(s)$.
If $\operatorname{Disp}_G(m) = 1$, we obtain the variance $v_\Delta(s^m)$. For simple discrete treatments, $\operatorname{Disp}_G(\cdot) = 1$ under classic stratified randomization (Example (ref)), so we can regard $v_\Delta(s)$ as the perfectly stratified variance. For general complex treatment spaces $\mathcal D \subseteq \mathbb{R}^m$, perfect stratification is impossible, but the corollary shows a sense in which we can approximate it by using high dispersion couplings.
Next, we illustrate how the theory developed in the previous section can be applied to compare the efficiency and robustness of designs based on LHS with designs based on rotation sampling (RS) and the Gaussian copula. The analysis shows that RS also generically produces high dispersion for large enough $k$, but is less robust to adversarial influence function shapes $s_i(\cdot)$ than LHS. By contrast, for moderate $k$ the Gaussian copula only produces high dispersion for approximately linear functions, a strong parametric restriction.
Recall from Example (ref) that a rotation sample $(U_i)_{i=1}^k \sim G_U$ lies on a randomly shifted equispaced grid, so that $U_i = U_j \oplus l/k$ for some $l \in [k]$ and $i, j \in [k]$, where $a \oplus b \equiv a + b \pmod 1$ for $a, b \in \mathbb{R}$. The canonical marginal for univariate rotation sampling is thus $F = \operatorname{Unif}[0, 1]$. We show that this coupling is efficient if influence functions $s_i(\cdot)$ are smooth, in a sense defined below.
Suppose a function $\phi(\cdot)$ is perfectly $1/k$-cyclic, with $\phi(x) = \phi(x \oplus 1/k)$ for all $x \in [0, 1]$. Then our samples $\phi(U_i) = \phi(U_j)$ for $i, j \in [k]$ are perfectly correlated under rotation sampling, effectively yielding a clustered sample. The low dispersion eigenspace turns out to be exactly this space of cyclic functions:
The acyclic subspace is defined as the orthogonal complement $E_{a} \equiv E_{c}^\perp$ in $L_0^2(F)$. Our analysis shows that $L_0^2(F) = E_c \oplus E_a$, which are eigenspaces of the coupling operator $U_G$ with dispersions $\operatorname{Disp}_G(E_a) = 1$ and $\operatorname{Disp}_G(E_c) = -(k-1)$. Projections on $E_c$ and $E_a$ have closed forms, with $P_c \phi(d) = k^{-1} \sum_{l=1}^k \phi(d \oplus l/k) - E_F[\phi(D)]$ the de-meaned cyclic average and residual $P_a \phi(d) = \phi(d) - k^{-1} \sum_{l=1}^k \phi(d \oplus l/k)$. Then by Theorem (ref), for any $\phi \in L^2(F)$,
Let $w_a(s)$ and $w_c(s)$ denote the approximation weights (Equation (ref)) on $E_a$ and $E_c$ respectively, so that $w_c(s)$ quantifies how cyclic the influence functions $s_i(\cdot)$ are, on average over $i \in [n]$. We apply Theorem (ref) to $G =$ RS, using the facts above.
The low dispersion cyclic space $E_c = E_c(k)$, where the design performs poorly, shrinks as $k$ increases. In particular, we have $E_c(r) \subseteq E_c(k)$ for $k \mid r$, so the weights satisfy $w_c^r(s) \le w_c^k(s)$. This shows a sense in which dispersion is monotonically increasing in tuple size $k$ for rotation sampling.
\paragraph{Robustness Comparison.}
Recall that for LHS, the low dispersion subspace $E_{hist}^\perp$ behaves like an iid design with $\operatorname{Disp}_G(E_{hist}^\perp) = 0$ (Example (ref)). By contrast, the low dispersion space under RS has $\operatorname{Disp}_G(E_c) = -(k-1)$, actually harming relative efficiency through the negative term in Equation (ref). This effect can also be seen in the nominal variance. Let group mean $\bar s_g(\cdot) = k^{-1} \sum_{i \in [k]} s_{ig}(\cdot)$ and define the group variance $v_{g}(s) \equiv (n/k)^{-1} \sum_g \operatorname{Var}_F(\bar s_g)$. By Corollary (ref), we have
This shows RS achieves the ideal perfectly stratified variance $v_\Delta(s^a)$ on $E_a$, but has variance $k \cdot v_{g}(s^c)$ on $E_c$, reducing the effective sample size by a factor of $k$. Thus, the RS coupling is less robust to adversarial influence function shapes than LHS. However, for this worst case to arise in practice, the influence functions $s_i(\cdot)$ must be strongly cyclic with high frequency, which is rare in typical social science applications.
We show that when $k$ is moderate, the Gaussian copula produces high dispersion only for approximately linear influence functions $s_i(\cdot)$, a restrictive condition relative to the general smooth functions that achieve high dispersion under the LHS and RS couplings. This suggests caution when using the Gaussian copula to generate dispersion in experimental design.
For the Gaussian copula, it is convenient to use canonical measure $F = \mathcal{N}(0, 1)$. In this case, the coupling operator $U_G$ coincides with the Mehler kernel operator mehler1866. We use the eigenbasis expansion $L_0^2(F) = \oplus_{m \ge 1} \operatorname{span}(h_m)$ for $U_G$, where $(h_m)_{m \ge 1}$ are the normalized probabilist's Hermite polynomials thangavelu1993. Each $h_m(x)$ is a polynomial of order $m$, for example, $h_1(x) = x$ and $h_2(x) = (x^2 - 1)/\sqrt{2}$. As shown in Theorem (ref), the dispersion of a polynomial $h_m$ can be obtained from its eigenvalue $\lambda_m$ under $U_G$, with $\operatorname{Disp}_G(h_m) = -(k-1)\lambda_m$. The projections onto $E_m = \operatorname{span}(h_m)$ are $s_i^m(\cdot) = \operatorname{Cov}_F(s_i, h_m) \cdot h_m(\cdot)$. An application of Theorem (ref) yields the following result.
Define the space of linear functions $E_L = \{\phi : \phi(d) = a + bd\}$. For any non-constant $\phi \in E_L$, we have $\operatorname{Disp}_G(\phi) = \operatorname{Disp}_G(h_1) = 1$, so this is a high dispersion subspace. By contrast, $|\operatorname{Disp}_G(h_m)| \le (k-1)^{-(m-1)}$ for any Hermite polynomial of order $m \ge 2$, which is rapidly decreasing as tuple size $k$ increases. For larger $k$, the Gaussian copula produces high dispersion only for the linear component of $s_i(\cdot)$, performing no better than iid randomization on $E_L^\perp$.
Corollary (ref) shows that, for moderate $k$, the Gaussian copula is efficient only if the influence functions $s_i(\cdot)$ are approximately linear. Note that $s_i(\cdot)$ may not be linear even if the potential outcomes $Y_i(\cdot)$ are linear in the treatment, since each $s_i(\cdot)$ also depends on the estimand and the design. For example, in the cash transfer experiment of Example (ref), if the responses are linear in the grant amount, $Y_i(d) = a_i + b_i d$, the influence function for the best linear approximation coefficient (Example (ref)) is $s_i(d) = Y_i(d)H(d)$ for $H(d)=(d-E_F[D]) / \operatorname{Var}_F(D)$, which is quadratic in $d$.
Our analysis of the various couplings above suggests that they can be sorted into two categories: non-parametric couplings like LHS and RS, which produce high dispersion under weak smoothness conditions on $s_i(\cdot)$, and parametric couplings like the Gaussian copula, which produce high dispersion only for a restricted class of influence functions with specific shapes. Here we formalize this claim.
For any $\phi: [0,1] \to \mathbb{R}$, define total variation $V_{[0,1]}(\phi) = \sup_\Pi \sum_{j=1}^r |\phi(t_j) - \phi(t_{j-1})|$, where the supremum is over all finite partitions $(t_j)_{j=0}^r$ of $[0,1]$. Define the class of bounded variation functions $\mathcal H(b) \equiv \{\phi: V_{[0,1]}(\phi) \le b\}$. This provides a weak notion of smoothness for functions on $[0,1]$, allowing for discontinuous functions (e.g.\ histograms) provided they don't oscillate too much over the interval $[0, 1]$. This condition is meaningful for both discrete and continuous treatments. For example, if $F$ is discrete, the effective influence function $\tilde s_i(\cdot)$ on $[0, 1]$ defined by $s_i(D) = s_i(F^{-1}(U)) = \tilde s_i(U)$ will be discontinuous, but may have small total variation.
The theorem shows that the Gaussian copula produces high dispersion only for functions $\phi(\cdot)$ that are well-approximated by linear functions, while the non-parametric LHS and RS couplings produce high dispersion for any function $\phi(\cdot)$ with bounded total variation.
To relate this to the variance of the estimator $\widehat{\theta}$, define a smoothness coefficient $\eta_{\operatorname{TV}}(s) \equiv E_n[V_{[0,1]}(s_i)^2] / v_{iid}(s)$, which is a normalized measure of average total variation of the $s_i(\cdot)$.
In particular, the theorem guarantees that LHS and RS improve precision relative to iid randomization when $Q_k(s) \ge \eta_{\operatorname{TV}}(s) / k$, so that the match quality is sufficiently high relative to the roughness of the influence functions $s_i(\cdot)$. This lower bound on efficiency can be improved by using stronger notions of smoothness. For example, by imposing uniform Lipschitz continuity on the $s_i(\cdot)$, the lower bound for LHS can be improved to $Q_k(s) - O(1/k^2)$. However, such smoothness conditions typically rule out discrete treatments, motivating the weaker total variation notion of smoothness that we use here.
We consider an asymptotic regime with a sequence of experimental populations and coupling designs indexed by $n$ with $n \to \infty$. Thus, in this section all variables are implicitly indexed by $n$, denoting their place in the sequence. For example, influence functions $s_i(\cdot) = s_i^n(\cdot)$ for $i \in [n]$, and similarly the design parameters $G = G(n)$ and $k = k(n)$. We often suppress the indexing for brevity.
The experimental design literature has noted an important tradeoff between efficiency and robustness, showing how balancing covariates to improve efficiency can entail a loss of robustness in adverse experimental settings efron1971coin,harshaw2024. Motivated by this, we study a robust notion of consistency, requiring that coupling designs perform well uniformly over a range of possible empirical settings harshaw2025general. Given a family of influence functions $\mathcal{S}$, the uniform mean square error $R_n(G, \mathcal S)$ for a coupling $G$ is
If $R_n(G, \mathcal{S}) = O(1)$, we say that $\widehat{\theta}$ is uniformly $\sqrt{n}$-consistent under the family $\mathcal{S}$ and the coupling sequence $G(n)$.
We can study the worst-case performance of our designs by requiring uniform consistency under weak restrictions on $\mathcal S$. For example, we can set $\mathcal S_{\mathrm{m}}$ to be the set of all $s_i(\cdot)$ with bounded average moments: $\mathcal S_{\mathrm{m}} = \{ s_1, \dotsc, s_n \in L^2(F) : E_n[\operatorname{Var}_F(s_i)] \le 1 \}$. This allows for settings with highly non-smooth influence functions and negative match quality, requiring good performance even if an adversary matched units into groups that are maximally dissimilar. We show that coupling designs are still uniformly $\sqrt{n}$-consistent over $\mathcal S_{\mathrm{m}}$ under weak conditions on the design parameters, providing a strong robustness guarantee.
The influence function configurations that attain the worst case rate in Theorem (ref) are generally pathological, with either highly non-smooth influences $s_i(\cdot)$ or negative match quality, which might not be relevant for empirical practice. For example, for $G = $ LHS, $R_n(G, \mathcal S_{\mathrm{m}})$ is attained by placing perfectly negatively correlated influence functions $s_i(\cdot)$ within each group, which assumes that our matching was not only ineffectual, but actually much worse than random. We can place mild regularity conditions on the set $\mathcal S_{\mathrm{m}}$ to rule out such pathological settings, obtaining bounds that are more informative about the performance of coupling designs in typical applications.
To illustrate this, we consider the nonparametric couplings $G \in \{\text{LHS}, \text{RS}\}$ under a restricted family of influence functions $\mathcal S_{\mathrm{r}}$ that have reasonable match quality and finite total variation. The next result shows that under such conditions, these couplings uniformly dominate the iid design.
We have $R_n(G_{iid}, \mathcal S_{\mathrm{r}}) = 1$, so the LHS and RS coupling designs dominate the iid design in a minimax sense over $\mathcal S_{\mathrm{r}}$ when $\bar \eta / k < q_0$, which holds for large enough $k$. Intuitively, this is because iid randomization has no way of exploiting the structure imposed on the influence functions by $\mathcal S_{\mathrm{r}}$. Moreover, after ruling out arbitrarily non-smooth functions $s_i(\cdot)$, rotation sampling $G=$ RS is uniformly $\sqrt{n}$-consistent over $\mathcal S_{\mathrm{r}}$ even as $k \to \infty$.
We provide conditions under which the estimator $\widehat{\theta}$ is asymptotically normal.
Part (2) of Assumption (ref) is a high level condition ruling out certain degenerate estimation problems. For example, under perfect homogeneity $s_i(d) = s(d)$ with $s(\cdot)$ Lipschitz continuous, one can show that super-efficient estimation is possible using the Latin hypercube coupling with $k=n$. Such perfectly homogeneous settings are not empirically relevant.
Since each of the $n/k$ groups has independent treatments, and each group's contribution to $\widehat{\theta}$ is bounded, the CLT follows from standard Lindeberg-Feller arguments provided the number of groups $n/k$ grows sufficiently fast. Assumption (ref) ensures this by requiring $k = o(n^{1/3})$.
We construct a variance estimator for $\operatorname{Var}_G(\widehat{\theta})$ using a collapsed strata approach hansen1953. Let $\pi: [n / k] \to [n / k]$ be a permutation with no fixed points, $\pi(g) \neq g$ for all $g \in [n / k]$, which associates each group $g$ with a paired group $\pi(g)$. Typically, this is a matching $\pi(\pi(g)) = g$, but we allow for non-matching permutations to accommodate, for example, an odd number of groups.
Let $\widehat{\theta}_{g} = k^{-1} \sum_{i = 1}^{k} s_{ig}(D_i)$ and $\theta_{g} = E_F[\widehat{\theta}_g] = k^{-1} \sum_{i = 1}^{k} \theta_{ig}$ for $\theta_{ig} = E_F[s_{ig}(D)]$. The variance estimator is the scaled squared difference between the mean outcomes in the paired groups:
Let $\Delta_n^2$ be the average squared difference in group effects across the paired groups:
The magnitude of the bias depends on the tuple size $k$ and the heterogeneity in group effects. Under $\sqrt{n}$-consistency, $\sigma_n^2 \asymp n^{-1}$, the normalized bias is $O(k \Delta_n^2)$. We generally expect $\Delta_n^2$ to be asymptotically bounded but not to vanish, so the variance estimator is typically biased upwards also asymptotically. The normalized bias is asymptotically bounded only if $k = O(1)$. However, the unnormalized bias is $O(k \Delta_n^2 / n)$, so it vanishes provided that $k = o(n)$. Because the variance estimator is conservative, confidence intervals constructed using $\widehat{\sigma}_n^2$ are asymptotically valid.
Coupling designs provide a powerful approach for improving the efficiency of experimental designs in complex treatment spaces. The core insight is that the mechanism underlying conventional stratification can be extended to complex treatment spaces by matching units into similar groups, then assigning within-group treatments that are highly dispersed over the space $\mathcal D$. The efficiency gain from coupling designs is proportional to the product of sample dispersion and match quality. The attained dispersion depends on the shape of the influence functions $s_i(\cdot)$, where high dispersion is achieved when the $s_i(\cdot)$ are well-approximated on the high-dispersion eigenspaces of the coupling operator $U_G$ associated with the relevant coupling design.
Several directions for future work remain. First, we focused on exchangeable couplings with fixed marginals for tractability and exposition. Alternative designs that use non-exchangeable couplings or allow for flexible marginal distributions may offer further efficiency improvements, but require new tools to analyze their properties. Second, the key insight of producing high sample dispersion among similar experimental units applies very broadly, and the stratified structure of coupling designs is not essential to leverage this insight. Alternative designs that do not partition units into groups, but still assign highly dispersed treatments to similar units, may be possible and could offer further efficiency improvements. Finally, integrating coupling designs with response-adaptive methods that update the treatment distribution over successive experimental waves is a promising avenue for combining the efficiency gains from both approaches.