The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
78,514 characters
The Limits of Experimental Design: Covariate Balance Beyond Low Dimension
\begingroup
\maketitle
\endgroup
\setcounter{footnote}{0}
\begin{abstract}
We study how fast experimental designs can approach the semiparametric efficiency bound in finite samples, as measured by the excess variance of unadjusted treatment effect estimation.
We prove an impossibility theorem: under weak conditions, no design can approach the variance bound uniformly over smooth outcome models unless covariate dimension $d \ll \log n$.
Even in experiments with thousands of units, this permits only a handful of covariates.
Motivated by this, we propose new designs based on discrepancy minimization that instead attempt to control imbalances over restricted-complexity nonparametric function classes.
Such designs achieve fast rates to their corresponding restricted efficiency targets, permitting $d \ll n$ covariates in an additive nonparametric specification.
They can also be combined with matching to protect against unmodeled outcome variation.
In simulations calibrated to 12 published experiments, our designs reduce variance relative to matched pairs randomization in every empirical setting.
\end{abstract}
\noindent{\emph{Keywords}: Semiparametric Efficiency, Kernel Methods, Discrepancy Minimization.}
\noindent{\emph{JEL Codes}: C14, C21, C90.}
\clearpage
\section{Introduction}\label{section:introduction}
Balancing covariates is a central task of experimental design.
Available methods range from stratified randomization, in use at least since \citet{fisher1926}, to rerandomization and newer discrepancy minimization methods \citep{li2018asymptotic,harshaw2024gsw}.
Effective covariate balance can improve precision, reducing the size and cost of the experiment required to study a given causal effect.
This is well known to practitioners.
For example, a survey of 200 randomized experiments in recent NBER working papers reported in \citet{cytrynbaum2022local} finds that 43\% used some form of stratification.
Recent theoretical work provides a strong justification for finely stratified designs like matched-pairs randomization.
Under appropriate conditions, such designs make simple unadjusted estimators semiparametrically efficient, automatically attaining the \citet{hahn1998} variance bound for estimating treatment effect parameters \citep{bai2021inference,bai2023efficiency}.
However, these guarantees rest on strong assumptions, requiring covariate differences within matched groups to vanish asymptotically.
This may not be possible with the number of covariates encountered in practice.
Note that nearest-neighbor distances for $n$ samples in $d$ dimensions typically decrease at rate $n^{-1/d}$, so halving the average distance between matches requires $2^d$ times as many observations.
To see the practical implications of this curse of dimensionality, consider the OpenResearch Unconditional Income Study, a prominent experiment that randomized 3,000 participants to a universal basic income program, distributing \$40 million in unconditional cash transfers \citep{broockman2024income}.
To randomize treatments, researchers first formed matched triples using several dozen covariates, then assigned one participant in each triple to receive \$1,000 per month for three years.
For the efficiency results above to be relevant here, we would need to obtain near-perfect matches across several dozen covariates, using only $n=3000$ samples.
Motivated by this gap between theory and practice, we develop a new finite-sample efficiency theory for experimental design.
For an unadjusted estimator $\widehat{\theta}$ of the ATE and efficiency bound $V^*$, we study the \emph{excess variance} $\mathcal V_n\equiv n\operatorname{Var}(\widehat\theta)-V^*$.
Let $\psi_i$ denote the covariates of unit $i$ and define the outcome model $m(\psi)\equiv E[(Y(1)+Y(0))/2|\psi]$.
Also let $Z_i\in\{-1,1\}$ be the centered treatment assignment with $E[Z] = 0$.
A calculation shows that excess variance is determined by in-sample imbalances in the outcome model:
\begin{equation}\label{equation:introduction-excess-variance}
\mathcal V_n = 4n \cdot E\bigg[\bigg(\frac{1}{n}\sum_{i=1}^n Z_i m(\psi_i)\bigg)^2\bigg].
\end{equation}
By studying such imbalances, we derive finite sample upper bounds on the excess variance $\mathcal V_n$, showing how rapidly different designs can approach $V^*$ for rough, smooth, and restricted-complexity nonparametric outcome models $m(\psi)$.
In several cases, we also obtain lower bounds establishing that these rates are minimax optimal.
More fundamentally, we prove an impossibility theorem showing that the curse of dimensionality illustrated above affects not only matching but any experimental design seeking to attain the classic variance bound $V^*$.
Under mild conditions, our results imply that no design can guarantee uniform asymptotic efficiency, even over smooth outcome models, beyond the regime $d=o(\log n)$, i.e.\ $d / \log n \to 0$ as $n \to \infty$.
For experiments with thousands of units, this permits only a few covariates.
This suggests an alternative design principle.
Rather than pursuing the classical efficiency bound $V^*$, we propose to balance lower-complexity nonparametric classes $\mathcal G$ that may still approximate the true outcome model well.
We develop new experimental designs based on this principle using nonparametric versions of the Gram-Schmidt walk (GSW) of \citet{bansal2019gram} and \citet{harshaw2024gsw}.
To measure their performance, we define the efficiency target $V^*+\Delta_{\mathcal G}$, where $\Delta_{\mathcal G}$ records outcome variation not captured by a restricted function class $\mathcal G$.
The lower complexity of such classes allows our designs to approach this target at much faster uniform rates and balance many more covariates than stratification can accommodate.
Designs targeting an additive nonparametric specification can attain the restricted efficiency target with $d=o(n)$ covariates, up to log factors.
We show these designs can also be combined with matching to preserve the fast rates above, while providing some protection against unmodeled outcome variation.
In an empirical application based on 12 randomized experiments recently published in economics, our preferred designs reduce estimator variance relative to classical matched pairs in every simulated empirical setting.
We make the following main contributions:
\begin{enumerate}
\item In Section~\ref{section:matching}, we study the finite-sample properties of matched pairs.
We show that no matching scheme can make excess variance $\mathcal V_n$ vanish faster than $n^{-2/d}$, even over linear, hence infinitely smooth, outcome models.
This extremely slow rate shows that matched pairs cannot realistically approach the efficiency bound in finite samples beyond the very low dimensional regime.
Despite this, we show matched pairs is minimax optimal in fixed dimension for rough outcome models, attaining $\mathcal V_n\asymp n^{-2\beta/d}$ over $\beta$-H\"older classes for $0<\beta\leq1$.
Its slow convergence is therefore the price of robustness under very weak assumptions.
\item Section~\ref{section:structure} introduces new randomization methods designed to exploit smoothness.
To do so, we develop a nonparametric Gram-Schmidt Walk that controls the worst-case imbalances over a chosen reproducing kernel Hilbert space.
We exhibit a design that adapts without prior knowledge to every Sobolev smoothness order $s>0$ and attains rate $\mathcal V_n \asymp n^{-(2s/d)\wedge1}$, up to logarithmic factors.
This rate is minimax optimal when $2s\leq d$ and nearly parametric when $2s>d$.
\item Section~\ref{section:structure} also establishes our main impossibility result: under weak conditions, no design can guarantee excess variance uniformly smaller than order $n^{-2s/d}$ over full-dimensional outcome models for any fixed smoothness $s>0$.
This precludes uniform asymptotic efficiency for any experimental design outside the very-low dimensional regime $d=o(\log n)$, even for smooth outcome models.
\item In view of this impossibility result, Section~\ref{section:restricted-efficiency} develops GSW designs using kernels tailored to lower-complexity nonparametric classes, including additive models of the form $g(\psi)=a+\sum_{j=1}^d g_j(\psi_j)$ and richer specifications with nonparametric bivariate interactions.
The additive and bivariate designs approach their restricted efficiency targets at fast rates, allowing $d=o(n)$ and $d=o(\sqrt n)$ covariates, respectively, up to logarithmic factors.
\item Finally, Section~\ref{section:restricted-efficiency} shows how to combine structured GSW with matching, retaining the fast rates above for low-complexity outcome components while also providing protection against unmodeled variation.
We develop a conservative design-based variance estimator for such methods that performs well in our simulations and empirical application.
\end{enumerate}
\subsection{Related Literature}
This paper is motivated by a recent literature showing that finely stratified randomization makes simple estimators attain the \citet{hahn1998} variance bound for the ATE.
Early contributions include \citet{bai2021inference}, \citet{bai2020pairs}, and \citet{cytrynbaum2022local}.
Building on \citet{armstrong2022}, \citet{bai2023efficiency} further show that this classical efficiency bound remains valid for general experimental designs and is attained by fine stratification in a GMM setting.
Related work includes \citet{imai2008}, \citet{fogarty2018}, \citet{pashley2021}, and \citet{kapelner2022}, among others.
Our asymptotic inefficiency result for matched pairs is an experimental analogue of the classical result of \citet{abadie2006large}, who show that slow nearest-neighbor convergence rates prevent matching from effectively debiasing observational estimates.
A related result is \citet{wang2022rerandomization}, who show a $d=o(\log n)$ threshold for balancing linear functions using rerandomization.
For comparison, our designs can balance additive nonparametric function spaces with $d = o(n)$ up to logarithms.
Our finite-sample perspective is related to that in \citet{kallus2018}, who also notes an analogue of Equation~\eqref{equation:introduction-excess-variance}.
His analysis leads to semi-deterministic designs, some of which are NP-hard to compute.
By contrast, we develop a full efficiency theory with explicit convergence rates of finite sample variances to both classical and newly introduced variance bounds.
We show that fast convergence rates can be achieved over restricted-complexity nonparametric spaces by computationally tractable randomized designs.
See the remarks in the text for a detailed comparison.
Our new designs in Section \ref{section:structure} fundamentally rely on the Gram-Schmidt walk vector balancing algorithm of \citet{bansal2019gram}.
\citet{harshaw2024gsw} first studied this algorithm in an experimental setting, proving a design-based oracle inequality for the mean squared error of an unadjusted estimator and using the algorithm to balance linear functions.
\citet{harshaw2021dissertation} first observed that the design can be kernelized, extending the original oracle inequality to this setting.
\citet{chen2026nonlinear} also suggest using GSW to balance nonlinear functions of the covariates.
To the best of our knowledge, we are the first to propose a finite sample efficiency theory for experimental design, using minimax analysis of the excess variance $\mathcal V_n$ as the benchmark.
The behavior of minimax rates over smoothness and restricted-complexity classes is a central organizing theme in nonparametric statistics \citep{stone1982optimal,yang1999information,gine2016mathematical}.
Our use of Sobolev norms to quantify smoothness also follows a long tradition in this literature.
For example, \citet{vdv2011information} show that convergence rates for Gaussian process regression with Mat\'ern priors are driven by a match between prior regularity and the Sobolev smoothness of the truth.
\citet{fischer2020sobolev} derive learning rates for kernel ridge regression in Sobolev norms.
See also \citet{kennedy2024minimax} for minimax rates for heterogeneous treatment effect estimation in observational studies.
The new designs in Section~\ref{section:restricted-efficiency} target balance over restricted function spaces.
Our main examples involve additive and bivariate nonparametric models motivated by functional ANOVA \citep{stone1985additive,hastie1986gam,sobol1993}.
Other prominent approaches to covariate balance include rerandomization \citep{morgan2012,li2018asymptotic} and optimization-based designs \citep{kasy2016,krieger2019}.
A complementary literature studies variance reduction by regression adjustment with moderately high-dimensional covariates \citep{bloniarz2016lasso,wager2016highdimensional,lei2021regression,lu2025debiased}.
Our work can be viewed as a randomization-based analogue of such results.
Finally, our variance estimator combines the design-based quadratic bound approach of \citet{harshaw2026optimized} with the collapsed-strata and pairs-of-pairs methods of \citet{hansen1953}, \citet{abadie2008}, and \cite{bai2021inference}.
\section{Setup}\label{section:setup}
Consider an experiment with $n$ units.
Denote potential outcomes $Y_i(a)$ for $a\in\{0,1\}$ and $i\in[n]$, letting $[r]=\{1,\ldots,r\}$ for any positive integer $r$.
\begin{assumption}\label{assumption:sampling}
Observations $\bigl(\psi_i,Y_i(0),Y_i(1)\bigr)\stackrel{\mathrm{iid}}{\sim} P$ for $i\in[n]$.
The covariate $\psi$ is supported on $[0,1]^d$ and has a density $p$ satisfying $0<\underline p\leq p(\psi)\leq\overline p<\infty$ for every $\psi\in[0,1]^d$.
For $a\in\{0,1\}$, $E_P[Y(a)^2]<\infty$.
\end{assumption}
Unless stated otherwise, every data-generating law $P$ in every model class $\mathcal P$ satisfies Assumption~\ref{assumption:sampling}.
We are interested in estimating the average treatment effect $\theta(P) = E_P[Y(1)-Y(0)]$.
Write $\psi_{1:n}=(\psi_1,\ldots,\psi_n)$ for the realized covariates.
For a realized treatment assignment $A_i\in\{0,1\}$, the observed outcome is $Y_i=Y_i(A_i)$.
For convenience, in what follows, we work with the coding $Z_i=2A_i-1\in\{-1,1\}$ and let $Z=(Z_1,\ldots,Z_n)$ be the full assignment vector.
\begin{defn}[Admissible Designs]\label{definition:admissible-designs}
A conditional assignment law $\sigma = Z | \psi_{1:n}$ is admissible $\sigma \in \mathcal D_n$ if $Z\mathrel{\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}}}(Y_i(0),Y_i(1))_{i=1}^n|\psi_{1:n}$ and $E[Z_i|\psi_{1:n}]=0$ for every $i\in[n]$.
\end{defn}
For example, admissible designs include iid assignment, matched-pairs and other forms of stratified randomization, rerandomization with symmetric acceptance rules, as well as the Gram-Schmidt walk designs introduced in Sections~\ref{section:structure} and~\ref{section:restricted-efficiency}.
We use Horvitz-Thompson estimation
\begin{equation}\label{equation:ht-estimator}
\widehat\theta = \frac{2}{n}\sum_{i=1}^nZ_iY_i.
\end{equation}
When exactly $n/2$ units are treated, as in matched pairs with even $n$, Equation~\eqref{equation:ht-estimator} is the usual difference in means.
Denote $\tau=Y(1)-Y(0)$ for the individual treatment effect and let $v_a^2(\psi)=\operatorname{Var}(Y(a)|\psi)$ for $a\in\{0,1\}$.
Let $\tau(\psi)=E[\tau|\psi]$.
At treatment probability one half, the semiparametric efficiency bound of \citet{hahn1998} for $\theta(P)$ is
\begin{equation}\label{equation:efficiency-bound}
V^*(P) = \operatorname{Var}\bigl(\tau(\psi)\bigr) +2E\bigl[v_1^2(\psi)\bigr] +2E\bigl[v_0^2(\psi)\bigr].
\end{equation}
This bound was originally derived for superpopulation settings with iid data.
Building on \citet{armstrong2022}, \citet{bai2023efficiency} show it holds over a broad class of experimental designs with dependent treatment assignments.
In what follows, we study the gap between the finite sample variance $n\operatorname{Var}_{P,\sigma}(\widehat\theta)$ for $\sigma \in \mathcal D_n$ and the variance bound $V^*(P)$.
To that end, define the \emph{excess variance}
\begin{equation}\label{equation:excess-variance}
\mathcal V_n(\sigma,P)\equiv n\operatorname{Var}_{P,\sigma}(\widehat\theta)-V^*(P).
\end{equation}
A sequence of designs $\sigma_n \in \mathcal D_n$ is pointwise asymptotically efficient at $P$ if the excess variance $\mathcal V_n(\sigma_n,P)\to0$ as $n \to \infty$.
Previous results have established this property for matched pairs designs and more general forms of fine stratification \citep{bai2021inference,bai2023efficiency,cytrynbaum2022local}.
Our results in Section~\ref{section:structure} also add to this list, providing a new family of pointwise efficient designs based on nonparametric GSW.
This classical formulation, however, holds $P$, and hence $d$, fixed as $n$ grows and does not characterize how quickly the variance bound is approached uniformly over a model class.
We show in what follows that this can lead to a very poor approximation to finite sample performance.
Instead, we propose to study finite-sample upper and lower bounds on excess variance $\mathcal V_n$ that make the dependence on sample size, covariate dimension, and outcome-model complexity explicit.
The following simple result identifies excess variance exactly with imbalances in $m(\psi) = E[(Y(1) + Y(0))/2 | \psi]$.
\begin{prop}[Excess Variance]\label{proposition:variance-decomposition}
For every $\sigma\in\mathcal D_n$, the excess variance is
\begin{equation}\label{equation:variance-decomposition}
\mathcal V_n(\sigma,P) = 4n \cdot E \bigg[ \Big( n^{-1} \sum_{i=1}^nZ_im(\psi_i) \Big)^2 \bigg].
\end{equation}
\end{prop}
This result recasts finite-sample efficiency as a function-balancing problem, with performance determined by the mean squared imbalance $n^{-1}\sum_{i=1}^nZ_im(\psi_i)$.
This equivalence was also previously noted in a slightly different form by \citet{kallus2018}.
Since the distribution $P$ is unknown at design time, a design cannot target $m(\psi)$ directly and must instead seek to control imbalance over a class of plausible models $P \in \mathcal P$.
Protecting against a very broad class may cause excess variance $\mathcal V_n(\sigma,P)$ to vanish too slowly for meaningful efficiency gains in finite samples.
Next, we show that matched pairs exhibits just such an extreme case of this tradeoff.
\section{The Limits of Matched Pairs Designs}\label{section:matching}
Matched pairs protects against broad classes of weakly regular outcome models, but this robustness comes at a cost.
In particular, no matched-pairs design can guarantee uniform asymptotic efficiency, even over linear outcome models, beyond the very-low dimensional regime $d=o(\log n)$.
\subsection{A Lower Bound for Matching}\label{subsection:matching-lower-bound}
Recall that a matched-pairs design first partitions the $n$ units into $n/2$ unordered pairs $(i_p,j_p)$, then randomly assigns one unit in each pair to treatment and the other to control with equal probability, independently across pairs.
The pairing itself may be any measurable function of the realized covariates $\psi_{1:n}$ and exogenous randomness $U_M$.
We write $\mathcal M_n$ for the resulting class of designs, with a given instance $M \in \mathcal M_n$.
Matched pairs is clearly admissible, $\mathcal M_n \subseteq \mathcal D_n$.
The definition implies $Z_{i_p}=-Z_{j_p}$ within each pair.
Independence between pairs and Proposition~\ref{proposition:variance-decomposition} then imply
\begin{equation}
\mathcal V_n(M,P)=\frac{4}{n}E \bigg[\sum_{p=1}^{n/2}\bigl(m(\psi_{i_p})-m(\psi_{j_p})\bigr)^2\bigg].
\end{equation}
Under Lipschitz continuity of $m(\psi)$, the within-pair differences $|m(\psi_{i_p})-m(\psi_{j_p})|$ are bounded by a constant multiple of the corresponding covariate distances $\|\psi_{i_p}-\psi_{j_p}\|_2$.
Pointwise asymptotic efficiency results therefore typically show $\mathcal V_n(M_n,P)\to0$ for each fixed law $P$ under a tight-matching condition of the form \citep{bai2021inference}:
\begin{equation}\label{equation:tight-matching}
\frac{1}{n}\sum_{p=1}^{n/2} \|\psi_{i_p}-\psi_{j_p}\|_2^2 \overset{p}{\to}0.
\end{equation}
See also \citet{cytrynbaum2022local} for algorithms guaranteeing such tight-matching conditions for general forms of fine stratification in fixed dimension.
By holding $d$ fixed, these asymptotic guarantees can obscure how quickly matching deteriorates with increasing covariate dimension in finite samples: excess variance may vanish asymptotically while remaining large at empirically relevant sample sizes.
The next result makes this limitation precise, even for linear outcome models.
To state the result, let $\mathcal P_{\mathrm{lin}}$ contain the laws satisfying Assumption~\ref{assumption:sampling} whose outcome model has the form $m(\psi)=a+\gamma'\psi$, with $\|\gamma\|_2\leq1$.
\begin{thm}[Matching Slow Rate]\label{theorem:matching-lower-bound}
For every even $n$ and every $d\geq1$,
\begin{equation}\label{equation:matching-lower-bound}
\inf_{M\in\mathcal M_n} \sup_{P\in\mathcal P_{\mathrm{lin}}}\mathcal V_n(M,P) \geq \frac{1}{2\pi e}n^{-2/d}.
\end{equation}
\end{thm}
Theorem~\ref{theorem:matching-lower-bound} applies to every pairing rule, not only optimal squared Euclidean distance matching as studied in \cite{bai2021inference}.
Changing the matching criterion therefore cannot improve this very slow excess variance rate.
Next, we use this finite-sample result to obtain alternative asymptotics for matched-pairs randomization in which $d=d_n$ may vary with $n$.
This allows us to formalize the sense in which matching is asymptotically inefficient beyond the very low-dimensional regime $d_n\ll\log n$.
\begin{cor}[Asymptotic Inefficiency]\label{corollary:matching-inefficiency}
If $d_n\geq c\log n$ for some fixed $c>0$, then
\begin{equation}\label{equation:matching-inefficiency}
\liminf_{n\to\infty} \inf_{M\in\mathcal M_n} \sup_{P\in\mathcal P_{\mathrm{lin}}}\mathcal V_n(M,P) \geq \frac{1}{2\pi e}e^{-2/c} >0.
\end{equation}
\end{cor}
Equivalently, sample size $n$ may be exponentially larger than covariate dimension $d_n$, yet the worst-case excess variance still remains bounded away from zero.
The same curse of dimensionality drives the non-vanishing bias of matching estimators in observational studies \citep{abadie2006large} and slow convergence rates for nearest-neighbor methods in nonparametric regression \citep{gyorfi2002}.
Figure~\ref{figure:matching-reversal-simulation} illustrates the finite-sample inefficiency suggested by both previous results.
Although recording more covariates steadily lowers the corresponding efficiency bound, the finite sample variance of matched-pairs eventually rises in both panels.
Linear GSW as in \cite{harshaw2024gsw} avoids this reversal under a linear outcome model and rerandomization delays it, but both offer almost no gain under an adversarial nonlinear additive model.
The additive nonparametric GSW design developed in Section~\ref{section:restricted-efficiency} below performs well in both settings.
\begin{figure}[!htbp]
\centering
\includegraphics[width=\textwidth]{figures/section3_matching_reversal.png}
\caption{Variance for $n=240$ units with independent uniform covariates. The left panel uses a normalized linear outcome model. The right panel replaces each linear coordinate with a nonlinear main effect while preserving the same coordinate-level variance profile. In both panels, the coordinate coefficients are proportional to $0.8^{j-1}$ and the dotted line is the efficiency bound using the first $d$ covariates. The secondary online appendix gives the exact models and simulation protocol.}
\label{figure:matching-reversal-simulation}
\end{figure}
\FloatBarrier
\subsection{Matching Targets Rough Outcome Classes}\label{subsection:rough-classes}
For $0<\beta\leq1$, define the H\"older coefficient $[m]_\beta=\sup_{x\neq y}|m(x)-m(y)|/\lVertx-y\rVert_2^\beta$.
In this section, we show that slow excess variance rates of form $n^{-2\beta/d}$ are minimax optimal for fixed $d < \infty$ when the outcome model is $\beta$-H\"older continuous.
We prove that stable matchings generically achieve these rates, even adaptively over $\beta\in(0,1]$, but cannot exploit additional smoothness beyond the Lipschitz endpoint $\beta=1$.
\medskip
\emph{Optimal Matching.} We describe a broad class of pairing rules that attain the excess variance rate $n^{-2\beta/d}$ for fixed $d$.
Optimal Euclidean matching, as in \citet{bai2021inference}, minimizes $\sum_{p=1}^{n/2}\|\psi_{i_p}-\psi_{j_p}\|_2^2$ over all pairings.
More generally, an experimenter may specify any symmetric matching cost $c(x,y)$ and use \citet{derigs1988} algorithm to compute optimal pairs
\begin{equation}
\min_{(i_p,j_p)_{p=1}^{n/2}}\sum_{p=1}^{n/2}c(\psi_{i_p},\psi_{j_p}).
\end{equation}
\emph{Stable Matching.}
In fact, global optimality is stronger than necessary.
A weaker notion is 2-swap stability, which requires that no two pairs can be rematched to reduce their combined cost.
For any two matched pairs $\{i,j\}$ and $\{k,l\}$, require that
\begin{equation}\label{equation:swap-stability}
c(\psi_i,\psi_j)+c(\psi_k,\psi_l) \leq \min\left\{ c(\psi_i,\psi_k)+c(\psi_j,\psi_l), c(\psi_i,\psi_l)+c(\psi_j,\psi_k) \right\}.
\end{equation}
We show that every 2-swap-stable pairing attains the excess variance rate $n^{-2\beta/d}$ for fixed $d$ whenever the matching cost satisfies a weak geometric condition, requiring $c(x,y)$ to be comparable above and below to Euclidean distance.
\begin{assumption}[Matching Cost]\label{assumption:matching-cost}
For some $q>0$ and constants $0<\underline\lambda\leq\overline\lambda<\infty$, the symmetric cost $c$ satisfies $\underline\lambda\|x-y\|_2^q\leq c(x,y)\leq\overline\lambda\|x-y\|_2^q$ for every $x,y\in[0,1]^d$.
\end{assumption}
We define $\Lambda=(\overline\lambda/\underline\lambda)^{1/q}$ as the condition number of the matching cost, which may depend on $d$.
It measures the loss from translating control of the matching cost $c(x,y)$ into control of Euclidean distance.
Thus $\Lambda=1$ for any power of Euclidean distance, while $\Lambda=\sqrt d$ for $c(x,y)=\|x-y\|_1^q$ or $c(x,y)=\|x-y\|_\infty^q$ for any $q>0$, for example.
Let $\mathcal P_\beta$ contain the laws satisfying Assumption~\ref{assumption:sampling} whose outcome model has H\"older coefficient $[m]_\beta\leq1$.
At $\beta=1$, this is the ordinary Lipschitz class, while smaller $\beta$ allow progressively rougher outcome functions.
\begin{thm}[Stable Matching]\label{theorem:stable-matching}
Fix an even $n\geq2$ and $0<\beta\leq1$, and suppose $d\geq3$.
Impose Assumption~\ref{assumption:matching-cost} and let $M_c$ be the matched-pairs design for any pairing that is 2-swap-stable for $c(x,y)$.
There is $C > 0$ depending only on $\beta$ such that
\begin{equation}\label{equation:stable-matching-high-d}
\sup_{P\in\mathcal P_\beta}\mathcal V_n(M_c,P) \leq C(d\Lambda^2)^{2\beta} \cdot n^{-2\beta/d}.
\end{equation}
\end{thm}
This is a finite-sample bound, but we can also interpret it asymptotically.
For example, in the low-dimensional regime with $d<\infty$ fixed while $n\to\infty$, any stable matching as above attains the excess variance rate $n^{-2\beta/d}$.
This holds simultaneously for every $\beta\in(0,1]$, showing that matching adapts to the unknown H\"older smoothness of $m(\psi)$.
Next, we establish a corresponding finite sample lower bound over all admissible designs, showing that no design can improve this rate uniformly over the H\"older class in the regime with low fixed dimension $d < \infty$.
\begin{thm}[H\"older Lower Bound]\label{theorem:design-agnostic-lower-bound}
For every $0<\beta\leq1$, there is a constant $C>0$ depending only on $(\beta, d)$ such that, for every $n\geq2$ and every $d\geq1$,
\begin{equation}\label{equation:design-agnostic-lower-bound}
\inf_{\sigma\in\mathcal D_n} \sup_{P\in\mathcal P_\beta}\mathcal V_n(\sigma,P) \geq C n^{-2\beta/d}.
\end{equation}
\end{thm}
The lower bound in Theorem~\ref{theorem:design-agnostic-lower-bound} is completely design agnostic.
It applies not only to existing procedures, such as stratified randomization, symmetric rerandomization \citep{li2018asymptotic}, and optimization-based balancing \citep{kallus2018}, but to every possible experimental design under the mild admissibility requirements above.
Combining Theorems~\ref{theorem:stable-matching} and~\ref{theorem:design-agnostic-lower-bound} immediately shows stable matched pairs is minimax rate optimal over rough H\"older classes in fixed low dimension $d < \infty$.
\begin{cor}[Minimaxity Over Rough Spaces]\label{corollary:minimax-matching}
Let $0<\beta\leq1$ and fix dimension $3\leq d<\infty$.
Then as $n\to\infty$, stable matching achieves minimax rate
\begin{equation}\label{equation:minimax-matching}
\inf_{\sigma\in\mathcal D_n} \sup_{P\in\mathcal P_\beta}\mathcal V_n(\sigma,P) \asymp n^{-2\beta/d}.
\end{equation}
\end{cor}
Matching therefore provides minimax protection adaptively over the rough H\"older scale $\beta \in (0, 1]$.
Despite this, Theorem~\ref{theorem:matching-lower-bound} also shows that no matched-pairs design can improve upon the rate $n^{-2/d}$ attained at $\beta=1$ even over the infinitely smooth linear class $\mathcal P_{\mathrm{lin}}$.
Thus matching cannot exploit smoothness beyond Lipschitz order.
Motivated by this, the next section develops new GSW designs that can adaptively exploit such higher order smoothness conditions, while also retaining good performance for rough outcome models.
\begin{remark}[Design Agnostic Bounds]
To the best of our knowledge, Corollary~\ref{corollary:minimax-matching} is the first design-agnostic minimax rate for convergence of finite sample variance to the \cite{hahn1998} variance bound under covariate-balancing randomization.
To prove the lower bound, we place a prior on an appropriate hard subfamily of $\mathcal P_\beta$.
We show that, for every treatment allocation $Z$, the squared imbalance $(n^{-1}\sum_{i=1}^nZ_im(\psi_i))^2$, averaged over this family, is bounded below by unavoidable local contributions for a constant fraction of sampled units.
Since this bound holds pointwise in $Z$, it remains valid after averaging over any admissible assignment distribution $\sigma\in\mathcal D_n$, providing a design-agnostic lower bound.
\end{remark}
\begin{remark}
\citet{kallus2018} proves that optimal matched pairs randomization minimizes a related conditional imbalance criterion uniformly over Lipschitz outcome models.
In his framework, the least-favorable outcome model $m$ may depend on both realized covariates $\psi_{1:n}$ and the designer's candidate treatment allocation $Z$.
By contrast, here we study the usual statistical minimax risk over a fixed superpopulation law.
He does not study upper or lower bounds on the convergence rate of the resulting excess variance, as we do here.
\end{remark}
\section{Balancing Smooth Functions}\label{section:structure}
Motivated by the inefficiency results for matched pairs above, we develop nonparametric discrepancy minimization designs that can exploit higher-order smoothness, building on the Gram-Schmidt walk of \cite{bansal2019gram} and \cite{harshaw2024gsw}.
We show that a certain kernelized GSW specification can achieve excess variance rate $n^{-2s/d}$ adaptively over smoothness parameters $0 < s \leq d / 2$ and $n^{-1}$ otherwise, up to logarithmic factors.
Unlike matched pairs, such global balancing schemes continue to improve under additional orders of Sobolev smoothness.
Despite the attractive properties of nonparametric GSW, we also establish our main impossibility result, showing that no admissible design can guarantee uniform convergence to the efficiency bound $V^*(P)$, even over smooth outcome models, outside the very-low dimensional regime $d=o(\log n)$.
\subsection{Sobolev Spaces and RKHS}\label{subsection:sobolev-models}
For nonparametric GSW, it will be convenient to parameterize regularity using Sobolev spaces.
For a multi-index $\alpha=(\alpha_1,\ldots,\alpha_d)$, write $|\alpha|=\sum_{j=1}^d\alpha_j$ and let $D^\alpha F$ denote the weak derivative of a function $F$ on $\mathbb{R}^d$.
Weak derivatives extend ordinary differentiation to functions whose rates of change exist only in an integrated sense, agreeing with ordinary derivatives whenever they exist.
Recall that $L^2(\mathbb{R}^d)=\{F:\mathbb{R}^d\to\mathbb{R}:\int_{\mathbb{R}^d}|F(\psi)|^2d\psi<\infty\}$ is the space of square-integrable functions.
For integers $k\geq1$, define the Sobolev space
\begin{equation}\label{equation:sobolev-whole-space}
H^k(\mathbb{R}^d) \equiv \left\{ F\in L^2(\mathbb{R}^d): D^\alpha F\in L^2(\mathbb{R}^d) \text{ for every }|\alpha|\leq k \right\}.
\end{equation}
The norm on $H^k(\mathbb{R}^d)$ can be expressed in terms of integrated derivatives
\begin{equation}\label{equation:sobolev-whole-space-norm}
\lVertF\rVert_{H^k(\mathbb{R}^d)}^2 \asymp \sum_{|\alpha|\leq k} \int_{\mathbb{R}^d}|D^\alpha F(\psi)|^2d\psi.
\end{equation}
In contrast to the H\"older classes in Section~\ref{section:matching}, which bound local changes uniformly across the covariate space, the Sobolev norm controls the aggregate size of these derivatives.
It therefore imposes an integrated global budget on local fluctuations.
Increasing the Sobolev order $k$ requires more weak derivatives to exist and to be square integrable, producing a progressively more regular class.
The definition extends to smoothness of every real order $s>0$.
Recall the Fourier transform $\widehat F(\omega)=(2\pi)^{-d/2}\int_{\mathbb{R}^d}e^{-i\omega'x}F(x)dx$.
Then $H^s(\mathbb{R}^d)$ consists of functions $F\in L^2(\mathbb{R}^d)$ for which
\begin{equation}\label{equation:sobolev-spectral-norm}
\lVertF\rVert_{H^s(\mathbb{R}^d)}^2 \equiv \int_{\mathbb{R}^d}|\widehat F(\omega)|^2(1+\lVert\omega\rVert_2^2)^sd\omega<\infty.
\end{equation}
At integer orders $s=k$, this yields the same space as the derivative-based definition above.
Increasing the smoothness parameter $s$ requires the high-frequency components of $F$ to decay faster.
The use of Sobolev spaces is convenient rather than essential: the supplementary online appendix shows how the corresponding rates for nonparametric GSW over Sobolev spaces extend to H\"older classes $0<\beta\leq1$ studied above.
\medskip
\emph{RKHSs.} We use reproducing kernel Hilbert space (RKHS) methods to construct a feasible balancing criterion that automatically adapts to the unknown Sobolev order of the outcome model $m(\psi)$.
Let $W$ be a continuous positive semidefinite kernel, so that for every $r\geq1$ and any $(x_i)_{i=1}^r \subseteq \mathbb{R}^d$, the Gram matrix $[W(x_i,x_j)]_{i,j=1}^r$ is positive semidefinite.
Its RKHS $\mathcal H_W$ is a Hilbert space of functions with the reproducing property $h(\psi)=\langle h,W(\psi,\cdot)\rangle_{\mathcal H_W}$ for any $h\in\mathcal H_W$, where $\langle\cdot,\cdot\rangle_{\mathcal H_W}$ denotes its inner product.
By Cauchy-Schwarz, the reproducing property implies $|h(\psi)|\leq W(\psi,\psi)^{1/2}\lVerth\rVert_{\mathcal H_W}$, so point evaluation is a continuous linear functional on $\mathcal H_W$.
The Mat\'ern family $W_d^{\mathrm{Mat},\nu}$ of kernels is indexed by the smoothness parameter $\nu>0$.
Its simplest member, obtained at $\nu=1/2$, is the exponential kernel $W_d^{\mathrm{Mat},1/2}(\psi,\psi')=\exp(-\lVert\psi-\psi'\rVert_2)$.
A formula for general $\nu$ is given in the secondary online appendix.
Write $\mathcal H_d^{\mathrm{Mat},\nu}\equiv\mathcal H_{W_d^{\mathrm{Mat},\nu}}$ for the associated RKHS.
The following standard correspondence links these RKHSs to the Sobolev spaces defined above.
\begin{prop}[Mat\'ern-Sobolev Correspondence]\label{proposition:matern-sobolev-equivalence}
For every $d\geq1$ and $\nu>0$,
\begin{equation}\label{equation:matern-sobolev-equivalence}
\mathcal H_d^{\mathrm{Mat},\nu} = H^{d/2+\nu}(\mathbb{R}^d), \qquad \lVertf\rVert_{\mathcal H_d^{\mathrm{Mat},\nu}} \asymp_{d,\nu} \lVertf\rVert_{H^{d/2+\nu}(\mathbb{R}^d)}.
\end{equation}
\end{prop}
In particular, every Sobolev space $H^s(\mathbb{R}^d)$ with $s>d/2$ is a Mat\'ern RKHS, obtained by taking $\nu=s-d/2$.
This correspondence is standard \citep{wendland2005scattered,kanagawa2018gaussian}.
We rederive it in the online appendix and record a uniform bound on the norm-equivalence as $\nu$ approaches zero, which is used in our later results.
\subsection{Controlling Imbalances for Smooth Functions}\label{subsection:unit-gsw}
Proposition~\ref{proposition:variance-decomposition} shows that excess variance is determined by the expected square of the outcome-model imbalance $n^{-1}\sum_{i=1}^nZ_im(\psi_i)$.
Since $m(\psi)$ is unknown at design time, the design cannot balance it directly.
We can form a feasible design criterion by instead controlling the worst-case imbalance over the unit ball of a chosen RKHS $\mathcal H_W$.
Consider an allocation $Z\in\{-1,1\}^n$ and let $W_n=[W(\psi_i,\psi_j)]_{i,j=1}^n$ be the corresponding kernel matrix.
By the reproducing property, one calculates imbalance
\begin{equation}\label{equation:rkhs-uniform-imbalance}
\mathcal I_W(Z) \equiv \sup_{\lVertg\rVert_{\mathcal H_W}\leq1} \left(\frac1n\sum_{i=1}^nZ_ig(\psi_i)\right)^2 = \frac1{n^2}Z'W_nZ.
\end{equation}
For $W=W_d^{\mathrm{Mat},\nu}$, Proposition~\ref{proposition:matern-sobolev-equivalence} identifies this criterion with the worst-case squared imbalance over a norm ball in $H^{d/2+\nu}(\mathbb{R}^d)$.
We therefore seek an admissible assignment law $\sigma \in \mathcal D_n$ under which the quadratic form $\mathcal I_W(Z)$ is small.
One can obtain such a design using the Gram-Schmidt walk (GSW) of \citet{bansal2019gram}, as adapted to experimental design by \citet{harshaw2024gsw}. This provides a computationally tractable assignment law that keeps $\mathcal I_W(Z)$ small while retaining the randomization needed for robustness and inference.
\emph{GSW Algorithm.} The Gram-Schmidt walk of \citet{bansal2019gram} takes as input an $n\times n$ positive semidefinite matrix $\Gamma_n$ with diagonal entries at most one.
It begins with the infeasible fractional treatment allocation $z^{(0)}=0\in[-1,1]^n$.
At each step $t\geq1$, it selects a pivot index $p_t$, which it retains until $z_{p_t}$ reaches $-1$ or $1$, and chooses an update direction $u$ minimizing $u'\Gamma_nu$, subject to $u_{p_t}=1$ and $u_j=0$ for coordinates that are already integral.
It randomizes between the largest feasible positive and negative steps in direction $u$, with probabilities chosen so the update has conditional mean zero.
At least one fractional coordinate reaches the boundary after each update, so the procedure terminates in finitely many steps in a feasible treatment allocation $Z \in \{\pm 1\}^n$.
See the secondary online appendix for a more detailed discussion.
\emph{Kernelization.}
\citet{harshaw2021dissertation} originally observed that GSW can be kernelized by constructing its input matrix $\Gamma_n$ from a kernel matrix.
We use the \emph{design kernel} $K=1+W$, where $W$ is a continuous positive semidefinite kernel satisfying $\sup_\psi W(\psi,\psi)\leq1$.
The added constant ensures that the associated RKHS can represent arbitrary constant levels of the outcome model $m$.
By the sums-of-kernels theorem \citep{paulsen2016rkhs}, $\mathcal H_K$ consists of functions $g=a+w$, where $a\in\mathbb{R}$ and $w\in\mathcal H_W$, with norm $\lVertg\rVert_{\mathcal H_K}^2=\min_{g=a+w}(a^2+\lVertw\rVert_{\mathcal H_W}^2)$.
Given covariates $\psi_{1:n}$ and a robustness parameter $\varphi\in(0,1)$, define the positive definite matrix
\begin{equation}\label{equation:kernel-gsw-gram}
\Gamma_n = \varphi I_n + \frac{1-\varphi}{2} \bigl[K(\psi_i,\psi_j)\bigr]_{i,j=1}^n.
\end{equation}
Because $K(\psi,\psi)\leq2$, the matrix $\Gamma_n$ is positive definite with diagonal entries at most one and therefore satisfies the GSW input conditions above.
The second term in Equation~\eqref{equation:kernel-gsw-gram} rewards balance over $\mathcal H_K$, while the identity term supplies robustness in outcome directions that are not represented by this RKHS.
Write $\mathrm{GSW}(K,\varphi)$ for the resulting assignment law, suppressing its dependence on $n$.
The mean-zero updates imply $E[Z_i|\psi_{1:n}]=0$ for every unit \citep{harshaw2024gsw}.
The update directions and step sizes depend only on the covariates and auxiliary randomness independent of the potential outcomes, so $Z \mathrel{\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}}} (Y_i(0),Y_i(1))_{i=1}^n | \psi_{1:n}$.
Thus, $\mathrm{GSW}(K,\varphi)$ is admissible.
\begin{remark}[Kernel Allocation]
The RKHS imbalance objective in Equation~\eqref{equation:rkhs-uniform-imbalance} also underlies the kernel allocation design of \citet{kallus2018}, which randomizes between $z^*$ and $-z^*$ for the optimal vector $z^*\in\arg\min z'W_nz$ with $z\in\{-1,1\}^n$.
Computing $z^*$ requires solving an NP-hard mixed integer program, and the resulting design has support of size two.
This limited randomization can have poor robustness properties and preclude standard inference \citep{kallus2021}.
Kernelized GSW avoids these issues through computationally tractable updates and broad randomization support.
\end{remark}
\subsection{Oracle Inequality and Smoothness Gains}\label{subsection:kernel-oracle-rates}
\citet{harshaw2024gsw} show that under GSW, the MSE of the Horvitz-Thompson estimator is bounded above by the loss of an implicit ridge regression of the outcome levels on the covariates.
\citet{harshaw2021dissertation} extends this inequality to kernelized GSW, yielding an implicit kernel-ridge objective.
Extending these results to our current superpopulation setting and combining with Proposition~\ref{proposition:variance-decomposition} yields the following oracle inequality for excess variance.
\begin{prop}[Oracle Inequality]\label{proposition:kernel-oracle}
Let $\sigma=\mathrm{GSW}(K,\varphi)$ as above.
Then
\begin{equation}\label{equation:kernel-oracle}
\mathcal V_n(\sigma,P) \leq \inf_{g\in\mathcal H_{K}} \left\{ \frac4\varphi E\left[(m(\psi)-g(\psi))^2\right] + \frac8{(1-\varphi)n} \lVertg\rVert_{\mathcal H_{K}}^2 \right\}.
\end{equation}
\end{prop}
The first term controls the residual $m-g$ not captured by the RKHS, while the second reflects the difficulty of balancing $g$, as measured by its RKHS norm.
The guarantee is strongest when $m$ can be approximated accurately by functions of moderate norm, which we show below occurs under sufficient smoothness restrictions.
Recall that a sequence of designs $\sigma_n$ is pointwise asymptotically efficient at $P$ if $\mathcal V_n(\sigma_n,P)\to0$.
Previously, this has only been established for matched pairs and other finely stratified designs \citep{bai2021inference,bai2023efficiency}.
Equation~\eqref{equation:kernel-oracle} shows asymptotic efficiency also holds for any kernelized GSW design whose RKHS is dense in $L^2(P_\psi)$.
\begin{cor}[Pointwise Efficiency]\label{corollary:kernel-pointwise}
Fix a design kernel $K$ and $\varphi\in(0,1)$, and, for each $n$, let $\sigma_n=\mathrm{GSW}(K,\varphi)$.
Suppose $\mathcal H_K$ is dense in $L^2(P_\psi)$.
As $n \to \infty$, excess variance
\[
\mathcal V_n(\sigma_n,P) \longrightarrow 0.
\]
\end{cor}
For every $\nu>0$, the Mat\'ern design RKHS is dense in $L^2(P_\psi)$ under every Borel law on $[0,1]^d$.
The Gaussian radial basis function kernel with length scale $\ell>0$ and $W_{\mathrm{RBF},\ell}(\psi,\psi')=\exp\{-\lVert\psi-\psi'\rVert_2^2/(2\ell^2)\}$ has the same density property \citep{scholkopf2002learning}.
Thus, Corollary~\ref{corollary:kernel-pointwise} applies to either kernel family.
However, as we argued in Section~\ref{section:matching}, such pointwise results can fail to yield meaningful guarantees on finite-sample efficiency.
Instead, we use the oracle inequality to derive uniform rates on the excess variance over suitable regularity classes.
The simplest example of such a bound follows from the inequality above when the outcome model belongs to a bounded ball in the design RKHS.
For a design kernel $K$ and $B>0$, define this class by $\mathcal P_K(B)=\left\{P:m\in\mathcal H_K,\ \lVertm\rVert_{\mathcal H_K}\leq B\right\}$.
\begin{cor}[RKHS Fast Rate]\label{corollary:kernel-ball}
Let $\sigma=\mathrm{GSW}(K,\varphi)$ as above.
Then
\begin{equation}\label{equation:kernel-ball}
\sup_{P\in\mathcal P_K(B)}\mathcal V_n(\sigma,P) \leq \frac{8B^2}{(1-\varphi)n}.
\end{equation}
\end{cor}
This shows kernelized GSW attains the fast parametric rate $n^{-1}$ uniformly over any fixed RKHS ball.
The next result is a direct Sobolev specialization: when smoothness $s>d/2$, the Mat\'ern-Sobolev correspondence places a bounded $H^s$ ball inside a bounded Mat\'ern RKHS ball, so the rate in Corollary~\ref{corollary:kernel-ball} transfers immediately.
To state this result on a scale that remains comparable across covariate dimensions, let $H_{\mathrm{av}}^s(\mathbb{R}^d)$ denote $H^s(\mathbb{R}^d)$ equipped with the dimension-normalized norm
\begin{equation}\label{equation:average-sobolev-norm}
\lVertF\rVert_{H_{\mathrm{av}}^s(\mathbb{R}^d)}^2\equiv\int_{\mathbb{R}^d}|\widehat F(\omega)|^2\left(1+\frac{\lVert\omega\rVert_2^2}{d}\right)^sd\omega.
\end{equation}
For fixed $d$, this is equivalent to the ordinary Sobolev norm and hence defines the same function space.
Define the probability model
\begin{equation}\label{equation:sobolev-outcome-class}
\mathcal P_d^s \equiv \left\{ P: m\in H_{\mathrm{av}}^s(\mathbb{R}^d),\ \lVertm\rVert_{H_{\mathrm{av}}^s(\mathbb{R}^d)} \leq 1 \right\}.
\end{equation}
\begin{cor}[Non-adaptive Fast Rate]\label{corollary:matern-supercritical}
Let $s>d/2$ and $\varphi\in(0,1)$.
Choose $\nu>0$ such that $d/2+\nu\leq s$, and let $K=1+W_d^{\mathrm{Mat},\nu}$.
For design $\sigma=\mathrm{GSW}(K,\varphi)$, there is a constant $C$ depending only on $(d,\nu,\varphi)$ such that
\begin{equation}\label{equation:matern-supercritical}
\sup_{P\in\mathcal P_d^s}\mathcal V_n(\sigma,P) \leq C n^{-1}.
\end{equation}
\end{cor}
A fixed Mat\'ern kernelized GSW design thus attains the parametric excess-variance rate $n^{-1}$ uniformly over Sobolev outcome models with smoothness $s>d/2$.
However, this simple result requires both enough smoothness $s>d/2$ and enough prior knowledge of $s$ to correctly choose the Mat\'ern tuning parameter $\nu$.
Also, for fixed $s$ the condition $s>d/2$ eventually fails along any sequence $d=d_n\to\infty$, so this analysis is unsuitable for studying asymptotics beyond low, fixed dimension.
Next, we remove both restrictions by constructing a single sequence of kernelized GSW designs that adapts to every $s>0$ without prior knowledge.
\subsection{Adaptation to Unknown Smoothness}\label{subsection:sobolev-adaptation}
For any fixed $c>0$, set $\nu_n=c/\log n$, so the Sobolev order $d/2+\nu_n$ of the Mat\'ern design RKHS approaches the critical boundary $d/2$ from above, independent of the unknown smoothness $s$.
This yields a design sequence that adapts to every $s>0$.
\begin{thm}[Adaptive Kernel GSW]\label{theorem:sobolev-gsw}
Fix $c>0$ and $\varphi\in(0,1)$, and set $\nu_n=c/\log n$, $K_n=1+W_d^{\mathrm{Mat},\nu_n}$, and $\sigma_n=\mathrm{GSW}(K_n,\varphi)$.
For every $s>0$, there is a constant $C$ depending only on $(s,c,\varphi,\overline p)$ such that, for every $d\geq1$ and $n\geq2$,
\begin{equation}\label{equation:sobolev-gsw}
\sup_{P\in\mathcal P_d^s}\mathcal V_n(\sigma_n,P) \leq C\left(\frac{\log n}{n}\right)^{(2s/d)\wedge 1}.
\end{equation}
\end{thm}
For fixed $d$, the exponent $2s/d$ increases continuously with smoothness until the rate reaches $\log n/n$ at $s=d/2$.
When $s>d/2$, this is only a logarithmic factor slower than the $n^{-1}$ rate in Corollary~\ref{corollary:matern-supercritical}.
Thus, a single design sequence adapts to every $s>0$ without prior knowledge of $s$, paying only a logarithmic factor for adaptivity.
\begin{remark}[Comparison to Matched Pairs]
In the supplementary online appendix, we show that, for fixed $d\geq3$, the same design attains the minimax rate $n^{-2\beta/d}$, up to logarithmic factors, uniformly over $\beta$-H\"older outcome models for every $0<\beta\leq1$.
Unlike matching, however, its rate continues to improve under Sobolev smoothness above order one and becomes nearly parametric once $s\geq d/2$.
\end{remark}
Figure~\ref{figure:smoothness-simulation} illustrates the finite-sample efficiency gains from smoothness for kernelized GSW vs.\ matched pairs.
The relative efficiency improvement grows as the outcome model becomes smoother.
The figure also shows that a design using the fixed parameter $\nu_0=1$ captures much of the gain achieved by Oracle Mat\'ern GSW, which uses the true smoothness parameter, and the advantage increases with $n$.
For comparison, the figure also includes a nonparametric rerandomization design.
In particular, we implement a kernelized version of the linear best-of-$m$ rerandomization in \cite{wang2025bestchoice}.
It selects, out of $m=10{,}000$ independent allocations $Z$, the one minimizing the worst-case imbalance $\mathcal I_W(Z)$ in Equation~\eqref{equation:rkhs-uniform-imbalance}, for the Mat\'ern kernel $W$ with the oracle value of $\nu$.
In contrast to its strong performance for the linear outcome model in Figure~\ref{figure:matching-reversal-simulation}, here rerandomization offers essentially no improvement over matched pairs, reflecting the difficulty of balancing a nonparametric function space through random search.
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{figures/section4_smoothness_simulation.png}
\caption{Prior-averaged excess variance for outcome models $m(\psi)$ drawn from a Mat\'ern Gaussian process with parameter $\nu$. Covariates are uniform on $[0,1]^d$. The left panel fixes $n=240$ and $d=8$. Fixed Mat\'ern GSW uses $\nu_0=1$ independently of the true $\nu$, while Oracle Mat\'ern GSW uses $\nu_0=\nu$. Kernel rerandomization selects, among 10,000 draws, the allocation minimizing $Z'W_nZ$, where $W_n$ is the Mat\'ern kernel matrix constructed using the true $\nu$. The right panel fixes $d=8$ and varies $n$ for $\nu\in\{0.5,1,4,16\}$. The secondary online appendix gives the exact model and simulation protocol.}
\label{figure:smoothness-simulation}
\end{figure}
\subsection{Impossibility Result}\label{subsection:design-agnostic-lower-bound}
The adaptive rate in Theorem~\ref{theorem:sobolev-gsw} still has the basic form $n^{-2s/d}$, which becomes vacuous when $s$ is fixed and $d=d_n\geq a\log n$ for some $a>0$.
Similar to the H\"older lower bound in Theorem~\ref{theorem:design-agnostic-lower-bound}, the next theorem shows that this is not a limitation of kernelized GSW: every admissible design faces the same dimensionality barrier.
\begin{thm}[Design-Agnostic Lower Bound]\label{theorem:sobolev-lower}
Fix $s>0$.
There is a constant $c_0>0$ depending only on $s$ such that, for every $d\geq1$ and every $n\geq2$,
\begin{equation}\label{equation:sobolev-lower}
\inf_{\sigma\in\mathcal D_n} \sup_{P\in\mathcal P_d^s}\mathcal V_n(\sigma,P) \geq c_0n^{-2s/d}.
\end{equation}
\end{thm}
The central implication is the dimensionality barrier.
For fixed $s$ and any sequence $d_n\geq a\log n$, Theorem~\ref{theorem:sobolev-lower} gives design-agnostic asymptotic inefficiency
\begin{equation}\label{equation:sobolev-dimensionality-barrier}
\liminf_{n\to\infty}\inf_{\sigma\in\mathcal D_n}\sup_{P\in\mathcal P_{d_n}^s}\mathcal V_n(\sigma,P)>0.
\end{equation}
Hence, no admissible design can guarantee convergence to the full efficiency bound beyond the low-dimensional regime $d_n=o(\log n)$.
Motivated by this impossibility result, the next section relaxes the pursuit of full efficiency by balancing lower-complexity nonparametric working models.
When such a working model closely approximates the true outcome model, its variance benchmark remains close to the Hahn bound but can be approached much more rapidly in finite samples.
\begin{remark}[Fixed-Dimensional Minimaxity]
For fixed $d < \infty$ and $s\leq d/2$, combining Theorems~\ref{theorem:sobolev-gsw} and~\ref{theorem:sobolev-lower} gives the following finite-sample minimax comparison.
There are constants $0<c_1<C_1<\infty$, depending on $(d,s,c,\varphi,\overline p)$, such that, for every $n\geq2$,
\begin{equation}\label{equation:sobolev-minimax}
c_1 n^{-2s/d} \leq \inf_{\sigma\in\mathcal D_n} \sup_{P\in\mathcal P_d^s}\mathcal V_n(\sigma,P) \leq C_1\left(\frac n{\log n}\right)^{-2s/d}.
\end{equation}
In low-dimensional asymptotics with $d<\infty$ fixed while $n\to\infty$, kernelized GSW is thus minimax rate optimal for every $s\leq d/2$, up to logarithmic factors.
\end{remark}
\section{Efficient Designs Beyond Low Dimensions}\label{section:restricted-efficiency}
In this section, we construct kernelized GSW designs that balance low-complexity nonparametric working models chosen to approximate the outcome model $m(\psi)$, rather than pursuing uniform convergence to the full efficiency bound.
By imposing additive or low-order interaction structure, these designs can accommodate many more covariates than matching or full-dimensional kernelized GSW, while still providing substantial nonparametric efficiency gains.
\subsection{Structured Nonparametric Balance}\label{subsection:structured-balance}
Kernelized GSW controls the worst-case imbalance $n^{-1}\sum_{i=1}^nZ_ig(\psi_i)$ over the unit ball of the RKHS supplied to the algorithm.
Section~\ref{subsection:unit-gsw} used a full-dimensional Mat\'ern space.
Here we instead use structured sum RKHSs, beginning with two examples.
\begin{ex}[Nonparametric Main Effects]\label{example:structured-priority}
We can control each covariate nonparametrically by targeting additive functions of the form $g(\psi)=a+\sum_{j=1}^dg_j(\psi_j)$.
\end{ex}
\begin{ex}[Bivariate Nonparametric Effects]\label{example:structured-bivariate}
Alternatively, we can control every bivariate interaction nonparametrically by balancing functions
\begin{equation}\label{equation:bivariate-nonparametric-class}
g(\psi)=a+\sum_{j=1}^dg_j(\psi_j)+\sum_{j<\ell}g_{j\ell}(\psi_j,\psi_\ell).
\end{equation}
\end{ex}
Both classes can be implemented directly with kernelized GSW through kernel addition.
Let $W_1,\ldots,W_B$ be continuous positive semidefinite kernels satisfying $W_b(\psi,\psi)\leq1$ and set $K(\psi,\psi')=1+B^{-1}\sum_{b=1}^BW_b(\psi,\psi')$.
Then the RKHS $\mathcal H_K$ consists of functions of the form $g=a+\sum_{b=1}^Bf_b$ with $f_b\in\mathcal H_{W_b}$.
For instance, one kernel corresponding to Example~\ref{example:structured-priority} is $K_1(\psi,\psi')=1+d^{-1}\sum_{j=1}^dW_1^{\mathrm{Mat},\nu}(\psi_j,\psi_j')$.
\@startsection{paragraph}{4}{\z@}
{\medskipamount}
{-\fontdimen2\font}
{\normalfont\normalsize\bfseries}{Variance Gap.}
The sum RKHS provides a general formula for low-complexity nonparametric approximations to $m$.
Let $\mathcal G$ denote the $L^2(P_\psi)$ closure of $\mathcal H_K$ and let $m_{\mathcal G}$ be the $L^2(P_\psi)$ projection of $m$ onto $\mathcal G$. We view this as a working model for $m$, with variance gap
\begin{equation}\label{equation:working-model-approximation-gap}
\Delta_{\mathcal G}(P)\equiv4E\bigl((m(\psi)-m_{\mathcal G}(\psi))^2\bigr).
\end{equation}
The corresponding variance target is $V^*(P)+\Delta_{\mathcal G}(P)$.
When $\mathcal G=L^2(P_\psi)$, we have $m_{\mathcal G}=m$ and $\Delta_{\mathcal G}(P)=0$, recovering the full efficiency bound.
More generally, when $\mathcal G$ closely approximates $m$, this target remains close to $V^*(P)$, but can often be approached much more rapidly in finite samples.
\emph{Regularity Conditions.}
As in Section~\ref{section:structure}, obtaining uniform rates requires regularity of the structured approximation $m_{\mathcal G}$.
The component representations in Examples~\ref{example:structured-priority} and~\ref{example:structured-bivariate} are not unique.
For example, main effects can be absorbed into the bivariate components.
Because of this, we impose regularity conditions on the canonical functional ANOVA decomposition \citep{sobol1993}.
In particular, expand
\begin{equation}\label{equation:canonical-low-order-decomposition}
m_{\mathcal G}(\psi)=a+\sum_{j=1}^dm_j(\psi_j)+\sum_{j<\ell}m_{j\ell}(\psi_j,\psi_\ell).
\end{equation}
We require $\int_0^1m_j(u)du=0$ as well as $\int_0^1m_{j\ell}(u,v)du=0$ for every $v$ and vice-versa.
These centering conditions make the corresponding component subspaces mutually orthogonal in $L^2([0,1]^d)$, so the decomposition is unique.
For smoothness $s>0$, we require the canonical components in Equation~\eqref{equation:canonical-low-order-decomposition} to satisfy the pooled Sobolev budget
\begin{equation}\label{equation:low-order-sobolev-budget}
\sum_{j=1}^d\lVertm_j\rVert_{H_{\mathrm{av}}^s(\mathbb{R})}^2+\sum_{j<\ell}\lVertm_{j\ell}\rVert_{H_{\mathrm{av}}^s(\mathbb{R}^2)}^2\leq1.
\end{equation}
The pooled bound prevents the aggregate magnitude and roughness of the structured projection from growing mechanically with the number of recorded covariates.
For the main-effects working model, the bivariate terms are omitted.
We state the main result for the richer bivariate design in Example~\ref{example:structured-bivariate}, leaving the main-effects case to Theorem~\ref{theorem:proof-general-low-order-gsw} in the appendix.
For this result, take $\mathcal G_2$ to be the $L^2(P_\psi)$ closure of the class in Equation~\eqref{equation:bivariate-nonparametric-class} and let $\mathcal P_{d,2}^{s}$ contain the laws satisfying Assumption~\ref{assumption:sampling}, $E[m(\psi)^2]\leq1$, and Equation~\eqref{equation:low-order-sobolev-budget}.
Fix $c>0$ and set $\nu_n=c/\log n$, so each component RKHS lies $c/\log n$ above its critical Sobolev order.
Writing $\psi_{j\ell}=(\psi_j,\psi_\ell)$ and similarly for $\psi'$, define
\begin{equation*}
K_n(\psi,\psi')=1+\binom d2^{-1}\sum_{j<\ell}W_2^{\mathrm{Mat},\nu_n}(\psi_{j\ell},\psi_{j\ell}').
\end{equation*}
For every $n$, the $L^2(P_\psi)$ closure of $\mathcal H_{K_n}$ is $\mathcal G_2$.
\Needspace{9\baselineskip}
\begin{thm}[Bivariate Effects]\label{theorem:low-order-gsw}
Fix $n\geq2$, $d\geq2$, and smoothness $s>0$.
Let $\sigma_n=\mathrm{GSW}(K_n,\varphi_n)$, where $1-\varphi_n=1/2\wedge(n^{-1}d^2\log n)^{1/2}$.
Then, for a constant $C>0$ independent of $n$ and $d$, uniformly over $P\in\mathcal P_{d,2}^{s}$,
\begin{equation}\label{equation:low-order-gsw-approximation}
\mathcal V_n(\sigma_n,P)\leq\Delta_{\mathcal G_2}(P)+C\left(\frac{d^2\log n}{n}\right)^{(s\wedge1)/2}.
\end{equation}
\end{thm}
Theorem~\ref{theorem:proof-general-low-order-gsw} in the appendix extends this result to allow a fixed block of priority covariates, such as baseline outcomes, to enter jointly and have a different smoothness order.
The extension adds the corresponding fixed-dimensional approximation rate while retaining the bivariate rate above for the remaining covariates.
When the working model is accurate, $\Delta_{\mathcal G_2}(P)$ is small, so the design approaches a variance target close to the Hahn bound $V^*(P)$ much more rapidly than the $n^{-2s/d}$ rate obtained by the full-dimensional Mat\'ern designs in Section~\ref{section:structure}, up to logarithmic factors.
The lower-complexity specification also accommodates many more covariates than those designs can.
Whereas all of the preceding designs, including matching and full-dimensional GSW, require the very low-dimensional regime $d_n=o(\log n)$, Equation~\eqref{equation:low-order-gsw-approximation} allows the number of covariates to grow nearly as fast as $\sqrt n$:
\begin{equation*}
\mathcal V_n(\sigma_n,P)\leq\Delta_{\mathcal G_2}(P)+o(1)\quad\text{if}\quad d_n=o(\sqrt{n/\log n}).
\end{equation*}
Note also that, as in Section~\ref{section:structure}, the design does not require prior knowledge of smoothness since the near-critical kernel sequence adapts automatically.
The calibration above uses $\varphi_n\to1$, but the theorem does not suggest that this is necessary.
In simulations, small fixed values of $\varphi$ continue to perform well, and we conjecture that this technical requirement is an artifact of the current analysis.
\emph{Well Specification.} In Theorem~\ref{theorem:proof-general-low-order-gsw} in the appendix, we provide a sharper rate under correct specification of the working model $m_{\mathcal G}$.
We show that if $m=m_{\mathcal G_2}$, then any fixed $\varphi\in(0,1)$ gives $\sigma_n=\mathrm{GSW}(K_n,\varphi)$ satisfying $\mathcal V_n(\sigma_n,P)\lesssim(d^2\log n/n)^{s\wedge1}$.
\smallskip
\emph{Main effects.}
The additive nonparametric design in Example~\ref{example:structured-priority} permits still faster dimension growth.
Let $\mathcal G_1$ denote the $L^2(P_\psi)$ closure of the main-effects models in Example~\ref{example:structured-priority}.
Writing $\sigma_{n,1}$ for the corresponding main-effects GSW design, Theorem~\ref{theorem:proof-general-low-order-gsw} shows that $\mathcal V_n(\sigma_{n,1},P)\leq\Delta_{\mathcal G_1}(P)+o(1)$ whenever $d_n=o(n/\log n)$.
Thus, replacing bivariate interactions by additive main effects increases the allowable dimension from $d = o(\sqrt{n/\log n})$ to $d = o(n/\log n)$, at the cost of the potentially larger variance gap $\Delta_{\mathcal G_1}(P)$.
\subsection{Robustifying Structured Designs with Matching}\label{subsection:matching-structured-balance}
Section~\ref{subsection:structured-balance} showed that targeting a lower-complexity nonparametric working model can deliver strong variance reductions even when $d$ grows far beyond $o(\log n)$.
These gains come with a residual variance gap $\Delta_{\mathcal G}(P)$ from predictable outcome variation outside the working model.
We now robustify structured GSW by combining it with matching.
Structured GSW retains the fast rates above for balancing the modeled component $m_{\mathcal G}$ globally, while matching provides some protection against unmodeled components $m - m_{\mathcal G}$ using local comparisons.
\@startsection{paragraph}{4}{\z@}
{\medskipamount}
{-\fontdimen2\font}
{\normalfont\normalsize\bfseries}{Matched GSW.} Let $M=\{(i_p,j_p):p\in[n/2]\}$ be a matching as in Section~\ref{section:matching}.
Any assignment treating exactly one unit in each pair can be written as $Z_{i_p}=T_p$ and $Z_{j_p}=-T_p$ for a \emph{pair orientation} $T_p\in\{-1,1\}$.
In classic matched-pairs, the $T_p$ are iid random signs.
Here, we propose to use kernelized GSW to generate them.
The imbalance can be written $n^{-1}\sum_{i=1}^nZ_ig(\psi_i)=n^{-1}\sum_{p=1}^{n/2}T_p(g(\psi_{i_p})-g(\psi_{j_p}))$.
Using this formula and a reproducing-property calculation as in Section~\ref{subsection:unit-gsw}, one can show that in this case the squared worst-case imbalance over the RKHS unit ball is $n^{-2}T'G_n^MT$, where $G_n^M\in\mathbb{R}^{(n/2)\times(n/2)}$ is the pair-difference matrix
\begin{equation}\label{equation:matched-pair-gram}
(G_n^M)_{pq}=K(\psi_{i_p},\psi_{i_q})-K(\psi_{i_p},\psi_{j_q})-K(\psi_{j_p},\psi_{i_q})+K(\psi_{j_p},\psi_{j_q}).
\end{equation}
As before, we seek to randomize in a way that makes $T'G_n^MT$ small.
This suggests normalizing $G_n^M$ and using it as the input to the Gram-Schmidt walk to generate the pair orientations $T$.
Let $\widehat\kappa_n=\max_p(G_n^M)_{pp}=\max_p\lVertK(\psi_{i_p},\cdot)-K(\psi_{j_p},\cdot)\rVert_{\mathcal H_K}^2$ be the largest squared RKHS distance within a matched pair and apply GSW to
\begin{equation}\label{equation:matched-gsw-gram}
\Gamma_n^M=\varphi I_{n/2}+\frac{1-\varphi}{\widehat\kappa_n}G_n^M.
\end{equation}
If $\widehat\kappa_n=0$, set $\Gamma_n^M=\varphi I_{n/2}$.
The matrix $G_n^M$ is positive semidefinite, and the diagonal entries of $\Gamma_n^M$ are at most one, as required by the algorithm.
Write $\mathrm{GSW}_M(K,\varphi)$ for the resulting admissible assignment law, which treats exactly half the units.
Let $\kappa_n=E[\widehat\kappa_n]$ be the expected largest squared RKHS distance within a matched pair.
Impose the Section~\ref{subsection:unit-gsw} normalization $K=1+W$, with $\sup_\psi W(\psi,\psi)\leq1$, so that $0\leq\kappa_n\leq4$.
For a square-integrable function $r(\psi)$, define its expected pair energy by
\begin{equation}\label{equation:matched-pair-energy}
Q_n^M(r)=n^{-1}E\bigg[\sum_{(i,j)\in M}\bigl(r(\psi_i)-r(\psi_j)\bigr)^2\bigg].
\end{equation}
The following paired analogue of Proposition~\ref{proposition:kernel-oracle} formalizes the local-global division of labor: matching controls within-pair variation locally, while kernelized GSW controls structured imbalance globally.
\begin{thm}[Matched Oracle Inequality]\label{thm:matched-gsw}
Fix $\varphi\in(0,1)$ and let $\sigma=\mathrm{GSW}_M(K,\varphi)$.
\begin{equation}\label{equation:matched-gsw-oracle}
\mathcal V_n(\sigma,P)\leq\inf_{f\in\mathcal H_K}\left\{\frac4\varphi Q_n^M(m-f)+\frac{4\kappa_n}{n(1-\varphi)}\lVertf\rVert_{\mathcal H_K}^2\right\}.
\end{equation}
\end{thm}
Relative to the unit-level oracle in Equation~\eqref{equation:kernel-oracle}, matching replaces the global approximation error $E[(m(\psi)-f(\psi))^2]$ by the within-pair energy $Q_n^M(m-f)$ and replaces the fixed RKHS penalty coefficient $8$ by $4\kappa_n$.
This can significantly attenuate both terms when the residual varies little within pairs and matched units are close in the kernel geometry.
\begin{ex}[Smooth Plus Rough Decomposition]
Suppose $m=a+g+r$ for some $g\in\mathcal H_K$ with $\lVertg\rVert_{\mathcal H_K}\leq B$ and a residual $r$.
Taking $f=g$ in Theorem~\ref{thm:matched-gsw} gives
\begin{equation}\label{equation:matched-smooth-rough-specialization}
\mathcal V_n(\sigma,P)\leq\frac4\varphi Q_n^M(r)+\frac{4\kappa_nB^2}{n(1-\varphi)}.
\end{equation}
The first term uses matching to control the residual $r$, while the second asks GSW to balance the structured component $g$.
One can show $\kappa_n=o(1)$ under suitable regularity conditions, providing a further rate enhancement.
\end{ex}
We now instantiate this oracle inequality using the main-effects design in Example~\ref{example:structured-priority}.
Recall that $\mathcal G_1$ is the $L^2(P_\psi)$ closure of the class in Example~\ref{example:structured-priority}.
Let $\mathcal P_{d,1}^{s}$ be defined as $\mathcal P_{d,2}^{s}$ with $\mathcal G_2$ replaced by $\mathcal G_1$ and the bivariate components omitted from Equations~\eqref{equation:canonical-low-order-decomposition} and~\eqref{equation:low-order-sobolev-budget}.
Using the same $\nu_n=c/\log n$ as above, define $K_{n,1}(\psi,\psi')=1+d^{-1}\sum_{j=1}^dW_1^{\mathrm{Mat},\nu_n}(\psi_j,\psi_j')$.
\Needspace{8\baselineskip}
\begin{thm}[Matched Main-Effects]\label{theorem:matched-low-order-gsw}
Fix $n\geq2$, $d$, smoothness $s>0$, and $\varphi\in(0,1)$.
Let $\sigma_n=\mathrm{GSW}_M(K_{n,1},\varphi)$.
For each $P\in\mathcal P_{d,1}^{s}$, define the residual $r=m-m_{\mathcal G_1}$.
Then, uniformly over $P\in\mathcal P_{d,1}^{s}$, for a constant independent of $n$ and $d$,
\begin{equation}\label{equation:matched-low-order-gsw}
\mathcal V_n(\sigma_n,P)\lesssim Q_n^M(r)+\left(\frac{d\log n}{n}\right)^{(2s)\wedge1}.
\end{equation}
\end{thm}
\smallskip
Unlike Theorem~\ref{theorem:low-order-gsw}, this bound contains no fixed variance gap: the residual enters through the within-pair energy $Q_n^M(r)$.
If $d_n=o(n/\log n)$, the structured term vanishes and
\begin{equation*}
\mathcal V_n(\sigma_n,P)\lesssim Q_n^M(r)+o(1).
\end{equation*}
If $Q_n^M(r)\to0$ as well, then $\mathcal V_n(\sigma_n,P)\to0$, and the design approaches the full Hahn bound.
If $d\geq3$, the matching $M$ is 2-swap-stable for a cost satisfying Assumption~\ref{assumption:matching-cost}, and $r$ is $\beta$-H\"older with coefficient $L_r$ for some $0<\beta\leq1$, then also
\begin{equation*}
\mathcal V_n(\sigma_n,P)\lesssim L_r^2(d\Lambda^2)^{2\beta}n^{-2\beta/d}+\left(\frac{d\log n}{n}\right)^{(2s)\wedge1}.
\end{equation*}
Thus only the residual pays the rough-class matching rate from Section~\ref{section:matching}, while GSW attains the faster main-effects rate for the structured component.
Figure~\ref{figure:structured-balance-simulation} illustrates the finite-sample gains from targeting a lower-complexity nonparametric function class when many covariates are recorded.
When the nonparametric main-effects working model is well specified, so that $m=m_{\mathcal G_1}$, both main-effects designs strongly outperform matched pairs and full-dimensional GSW.
When the true outcome model additionally contains interactions, $m\neq m_{\mathcal G_1}$ then matched pairs and full-dimensional GSW initially outperform main-effects GSW but deteriorate as the number of covariates increases and are eventually overtaken by it.
The version of main-effects GSW robustified with matching performs the best overall.
It interpolates between the two regimes, attenuating the omitted interaction terms when the dimension is small but also preserving global balance when the dimension is large.
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.88\textwidth]{figures/section5_structured_balance_simulation.png}
\caption{Excess variance at $n=240$ for the main-effects design in Example~\ref{example:structured-priority}. The panels correspond to well specification and additional interactions between the covariates, respectively. The dotted line is the working-model excess-variance target $\Delta_{\mathcal G_1}(P)$. The main-effects designs use the additive kernel $K_1$ from Example~\ref{example:structured-priority}, an average of univariate Mat\'ern kernels. All GSW designs use $\varphi=0.03$ and $\nu=1$. The secondary online appendix gives the exact model and simulation protocol.}
\label{figure:structured-balance-simulation}
\end{figure}
\FloatBarrier
\section{Empirical Application}\label{section:empirical-application}
In this section, we evaluate our new designs using Monte Carlo simulations calibrated to 12 randomized experiments in recently published or forthcoming economics papers.
We compare classical matched-pairs randomization with the robust GSW designs from Section~\ref{section:restricted-efficiency}, with kernels ranging from the full-dimensional Mat\'ern to the reduced complexity nonparametric main effects specification.
The nonparametric additive and bivariate GSW designs in Section~\ref{subsection:structured-balance} reduce estimator variance relative to classical matched pairs in every experiment, while full-dimensional GSW as in Section~\ref{section:structure} performs similarly to matched pairs on average.
The matched versions of GSW in Section~\ref{subsection:matching-structured-balance} provide further robustness: every design reduces estimator variance relative to classical matched pairs in every experiment.
The main-effects specification in Example~\ref{example:structured-priority} delivers the largest variance reduction overall, which is reflected in shorter confidence intervals using our inference methods.
We selected 12 experiments from public replication data made available through the AEA, J-PAL, the World Bank, the \emph{Journal of Development Economics}, and the \emph{Quarterly Journal of Economics}.
The screening considered data availability, sample size, and the existence of meaningful pretreatment covariates and was fixed before any design comparisons.
Appendix~\ref{appendix:empirical-application} describes the selection procedure in detail.
The resulting experiments have sizes $n\in[96,1{,}200]$ and covariate dimensions $d\in[7,51]$.
\emph{Data Calibration.} For each experiment, we calibrate the Monte Carlo population to its empirical covariate distribution and estimated conditional potential-outcome laws.
For continuous outcomes, we model $Y_i(a)=m_a(\psi_i)+\sigma_a(\psi_i)\epsilon_{ia}$, $a\in\{0,1\}$, where each $\epsilon_{ia}$ has a normal distribution.
For binary outcomes, the conditional law is Bernoulli with success probability $m_a(\psi_i)$.
We use cross-validation to select a regression model from linear regression, elastic net with interactions, random forests, and gradient boosting, while a shallow random forest estimates $\sigma_a^2$ from squared cross-fitted residuals.
Clustered experiments are aggregated to the level of randomization using study-specific rules.
We then draw covariates $\psi_{1:n}$ from a smoothed bootstrap and potential outcomes from the fitted conditional laws.
See Appendix~\ref{appendix:empirical-application} for details.
\smallskip
\textbf{Designs and Estimators.}
We compare classical matched-pairs randomization with both unit-level and matched versions of the full-dimensional, bivariate, and additive GSW designs from Section~\ref{section:restricted-efficiency}.
All GSW designs set $\varphi=0.03$ and use the adaptive near-critical Mat\'ern calibration $\nu_n=1/\log n$.
We form a common set of pairs using the matching algorithm in \citet{cytrynbaum2022local}.
Under classical matched-pairs, the pair orientations $T_p\in\{-1,1\}$ are independent, whereas matched GSW chooses them jointly using the Gram-Schmidt walk.
Thus, in the simulation the matched designs differ only in how they jointly orient a common set of pairs.
All designs use difference in means as the point estimator.
\smallskip
\textbf{Inference.} Matched GSW makes the pair orientations dependent, so the usual variance estimators available for classical matched pairs are not valid.
Because of this, in Appendix \ref{appendix:matched-gsw-inference} we develop a new covariance corrected estimator for these designs.
Our main result in Theorem~\ref{theorem:matched-gsw-inference} establishes design-based finite-sample conservativeness for inference on the sample average treatment effect (SATE), while Proposition~\ref{theorem:matched-gsw-ate-calibration} calibrates it for asymptotically non-conservative inference on the ATE.
For each matched GSW design, we report the stabilized estimator $\widehat{\mathcal V}_{\mathrm{stab}}$, implemented using an allocation bank of size $2J = 1024$.
See the appendix for details.
For each simulated empirical environment, we draw 400 independent covariate populations and, within each population, 120 fresh potential outcome schedules and assignments.
We report estimator variance reduction relative to matched pairs, coverage of nominal 95\% ATE confidence intervals, and reduction in CI width.
\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{figures/empirical_application.png}
\caption{GSW relative to matched pairs across 12 simulated experiments, with sample size and covariate dimension shown at left. Panels A and B report reductions in estimator variance for unit-level and matched GSW, respectively. Panel C reports the reduction in mean 95\% ATE confidence-interval width for matched GSW.}
\label{figure:empirical-application}
\end{figure}
\@startsection{paragraph}{4}{\z@}
{\medskipamount}
{-\fontdimen2\font}
{\normalfont\normalsize\bfseries}{Results.}
Figure~\ref{figure:empirical-application} documents broad efficiency gains: aside from unit-level full-dimensional GSW, every GSW design reduces estimator variance relative to classical matched pairs in every experiment, and all three matched designs do so uniformly.
The ranking among the unit-level designs illustrates the benefits from balancing lower-complexity nonparametric function spaces.
Recall that full-dimensional GSW targets the true semiparametric efficiency bound $V^*(P)$, whereas the bivariate and additive designs target the potentially larger efficiency bounds $V^*(P)+\Delta_{\mathcal G}(P)$.
These lower-complexity spaces are much easier to balance in finite samples.
As suggested by our rate theory, this advantage more than offsets the higher asymptotic variance target.
In particular, the nonparametric main effects design performs best in all settings, despite having the largest asymptotic variance target.
Panel C focuses on inference for matched GSW.
We develop new inference methods for these designs in Appendix~\ref{appendix:matched-gsw-inference}, which exploit the variance reduction to provide shorter confidence intervals.
Average coverage is 94.8\% for each matched design, with experiment-level coverage between 94.2\% and 95.2\%.
Taken together, these simulations show how the theoretical efficiency gains studied in the previous sections persist at the sample sizes and under the covariate distributions and potential-outcome relationships in actual experiments.
They support matched additive GSW as our preferred empirical design, since it combines the strongest efficiency gains with shorter confidence intervals and coverage close to nominal across all 12 experiments.
\section{Discussion and Recommendations for Practice}\label{section:discussion}
Both our theoretical and empirical results support a clear practical recommendation.
GSW designs targeting lower-complexity nonparametric function classes reduce estimator variance relative to classical matched pairs in every setting in our empirical application.
This reflects a basic difference between the designs.
Matching is subject to a severe curse of dimensionality, whereas discrepancy minimization designs like GSW can trade off covariate imbalances globally.
Combining the additive design with matching retains fast rates for nonparametric main effects while adding protection against misspecification, making matched additive GSW our preferred design overall.
Our covariance-corrected inference method also translates these efficiency gains into shorter confidence intervals, with coverage close to nominal.
Several directions merit further study.
The Gram-Schmidt walk is only one algorithmic implementation of global discrepancy minimization, and other algorithms different structured function classes may yield further improvements.
A more detailed investigation of inference for nonparametric GSW is beyond the scope of this paper but would be valuable, paralleling recent detailed investigations of inference for classical matched pairs designs \citep{fogarty2018,bai2026graph}.
Finally, it would be interesting to extend these methods beyond binary treatments to settings with multiple treatments or more complex experiments, as in the coupling-design approach of \citet{cytrynbaum2026coupling}.
\bibliographystyle{apalike}
\bibliography{references}