Experimental Design

Motivation

Online platforms run randomized experiments on subjects (users, sessions, orders) that arrive one at a time and must be assigned to a control or a treatment group on arrival, before later subjects are seen. The precision of a treatment-effect estimate depends on how similar the two groups are in their observed covariates. When the groups are badly balanced, the estimator's variance grows, and a large literature on covariate balancing (blocking, rerandomization, matched pairs) addresses this for offline experiments, where all covariates are known in advance.

J. Zhao and Z. Zhou (arXiv:2201.12936, Management Science, 2025) cast the online problem as an online matching problem. The cost of an assignment is the discrepancy, the weight of a minimum-cost perfect matching between the two groups, and they compare three designs against an oblivious adversary that fixes the covariate sequence in advance:

• the clairvoyant matched-pair design;
• the completely randomized design;
• their online pigeonhole design.

For one continuous covariate the pigeonhole design attains $O(T^{1/4}(\log T)^{1/2})$ against $\Theta(T^{1/2})$ for complete randomization. This mission concerns the multi-dimensional result, Theorem 9.

Setting

$T$ subjects arrive in order; $T$ is even. Subject $t$ carries a covariate vector $x_t$ with $q$ discrete coordinates, the $i$-th in a finite set $S_i\subseteq[0,1]$ of size $m_i$, and $p\ge2$ continuous coordinates in $[0,1]$. Distances are Euclidean in $\mathbb R^{q+p}$. An assignment puts each subject in control ($A$) or treatment ($B$).

The discrepancy of a split with $|A|=|B|$ is

$$d_{A,B}=\min_{\sigma:A\to B\text{ bijective}}\ \sum_{t\in A}\lVert x_t-x_{\sigma(t)}\rVert_2 ,$$

the value of the assignment problem (1) of the paper.

The completely randomized design picks $A$ uniformly among the $\binom{T}{T/2}$ subsets of size $T/2$; $\mathrm{ALG}^C(x)$ is the mean of $d_{A,B}$.

A pigeonhole design fixes a partition of the covariate space into pigeonholes. When subject $t$ arrives:

1. if one group already has $T/2$ subjects, $t$ joins the other group;
2. otherwise $t$ joins the group that has fewer subjects in its own pigeonhole, and on a tie a fair coin decides.

Theorem 9 uses the pigeonholes

$$\mathcal P(\mathcal S)=S_1\times\cdots\times S_q\times\Big\{[0,c^{1/p}T^{-\varphi}),\ \dots,\ [1-c^{1/p}T^{-\varphi},1]\Big\}^p ,$$

uniform cells of edge $c^{1/p}T^{-\varphi}$ on each continuous axis, for parameters $\varphi\in(0,1)$ and $c>1$. The expected discrepancy is the mean of $d_{A,B}$ over the design's coins.

Formalization targets

Goal: Theorem 9 (p. 26)

There are $C>0$ and $T_0$, depending only on $p,q,S_1,\dots,S_q,\varphi,c$, such that for every even $T\ge T_0$ and every admissible arrival sequence,

$$\mathbb E[d_{A,B}]\le C\Big(T^{1-\varphi}+T^{(p-1)\varphi}\,\lambda_p(T)\Big),\qquad \lambda_2(T)=(\log T)^{3/2},\quad \lambda_p(T)=1\ (p\ge3).$$

Milestones

1. Lemma EC.1: Hoeffding for a binomial, $\Pr(|H(n)-pn|\le\sqrt{n\log n})\ge1-2/n^2$.
2. Lemma 5, $p=2$: if $T=m^2$ and the subjects occupy each cell of the $m\times m$ grid exactly once, then $\mathrm{ALG}^C(x)=O(T^{1/2}\log T)$.
3. Lemma 5, $p\ge3$: on the full grid of $T=m^p$ cells, $\mathrm{ALG}^C(x)=O(T^{(p-1)/p})$.
4. Balance: every run of a pigeonhole design has $|A|=|B|=T/2$.
5. Component 1: the subjects paired within their pigeonhole before the stopping time $\tau$ contribute at most $\tfrac12c^{1/p}\sqrt p\,T^{1-\varphi}$, on every run.
6. Unmatched subjects: at most one unmatched subject per pigeonhole, $H^0_\tau+H^1_\tau$ at most the number of pigeonholes, and $|H^0_\tau-H^1_\tau|=T-\tau$.
7. Component 2: the remaining subjects contribute $O(T^{(p-1)\varphi}\lambda_p(T))$ in expectation.

Corollary 2 (the case $\varphi=1/p$) is an extra theorem after the goal.

Significance

At $\varphi=1/p$ Theorem 9 gives $O(T^{(p-1)/p})$ for $p\ge3$. This is the order $\Theta(T^{(p-1)/p})$ of the matched-pair design (Theorem 7), which sees all covariates in advance. An online design that commits to each assignment on arrival therefore loses at most a constant factor in dimension $p\ge3$, and at most a $(\log T)^{3/2}$ factor at $p=2$. For the completely randomized design no comparable upper bound is known for general sequences (p. 21). Lemma 5 is the paper's bound for complete randomization on the special family of full grids, and is of independent interest as a random-labelling matching estimate.

The results are proved on paper; none of them has a machine-checked proof, and the platform holds no item on them. A formalization would close the gaps of the printed argument:

• Lemma 5's induction for $p=2$ ends at $\sqrt T(\log T)^{3/2}$ rather than the stated $\sqrt T\log T$.
• Component 2 applies Lemma 5 to unmatched subjects, which occupy distinct cells but do not fill the grid.
• The case $q\ge1$ is argued per discrete support, although the balancing rule acts across supports.

Every statement in the mission is believed true.

Difficulty

The matched part (Component 1) is deterministic: paired subjects share a pigeonhole of diameter $O(T^{-\varphi})$. The difficulty lies in the subjects left over: the unmatched ones and the subjects after $\tau$. Their labels come from coins, but which subjects end up unmatched, and when $\tau$ occurs, depend on the same coins. So the leftover set is not a fixed set with independent fair labels, and the conditional binomial law used on p. ec22 does not hold literally.

A second obstacle is the random-matching estimate itself. Bounding the expected optimal matching between randomly labelled points of a grid at the rates $T^{1/2}\log T$ ($p=2$) and $T^{(p-1)/p}$ ($p\ge3$) is a multi-scale problem of Ajtai–Komlós–Tusnády type. A single-scale count of imbalances per cell gives only $T^{1/2}$ times a power of $T$ that is too large.

Formalization scope

• Subjects are indexed $0,\dots,T-1$; an assignment is a map to Bool (false = control). The covariate space is EuclideanSpace ℝ (Fin q ⊕ Fin p). Admissible sequences have discrete coordinates in $S_i$ and continuous ones in $[0,1]$, with $S_i\subseteq[0,1]$ (the paper's normalisation).
• Discrepancy: an infimum over bijections between the groups. An unbalanced split has the junk value $0$. The designs, the matched set and its complement are balanced on every run, so the junk is never reached.
• Cells: the printed edge $c^{1/p}T^{-\varphi}$ gives a partition only when $c^{-1/p}T^{\varphi}$ is an integer. The formalization uses $m=\lceil c^{-1/p}T^{\varphi}\rceil$ cells per axis of edge $1/m\le c^{1/p}T^{-\varphi}$, with the last cell closed. A pigeonhole is keyed by the exact discrete values and the cell index on each axis. The partition is the printed one, not a quantified family.
• Randomness: the pigeonhole design uses one coin per subject, consulted only on a tie. Its expected discrepancy is the average over all $2^T$ coin vectors, which is exactly the paper's law. The completely randomized design averages over all $\binom{T}{T/2}$ control sets.
• Asymptotics: $O(f(T))$ means $\exists C>0,\exists T_0,\forall T\ge T_0$ even, $\forall x$ admissible. $C$ and $T_0$ are chosen before $T$ and $x$ and depend only on the fixed parameters. $\log$ is the natural logarithm.
• Lemma 5's family $\mathcal X$: $T=m^p$ subjects in $[0,1]^p$, one per cell, stated as injectivity of the cell map.
• Stopping objects: $\tau$ is the first $n$ at which a group holds $T/2$ subjects. Within a pigeonhole the subjects before $\tau$ are paired in arrival order, and the last one is unmatched when their number is odd.
• Corrections to the page: the second cell printed as $[T^{-\varphi},2c^{1/p}T^{-\varphi})$ on p. ec21 is read as $[c^{1/p}T^{-\varphi},2c^{1/p}T^{-\varphi})$. The family's $x'_t\notin\mathbb X$ is read as $x_{t'}\notin\mathbb X$. The bound $H\le\frac1cT^{p\varphi}$ becomes $H\le(\prod_i m_i)\,m^p$ under the ceiling and the discrete supports.

The following formalizations would trivialize the mission and are ruled out:

• a bound on the discrepancy of one run, or of the best run, instead of the average over all coin vectors;
• constants that depend on $T$ or $x$;
• an unbalanced design that exploits the junk value $0$;
• an arbitrary partition in place of the printed grid.

A complete development needs a binomial Hoeffding bound and an exact analysis of the design's recursion (balance, alternation within a pigeonhole before $\tau$). It also needs a multi-scale bound on random bichromatic matchings of grid points. That bound is reusable for Lemma 5 and for matching estimators generally. Contributions on any milestone are welcome, as are proofs of Lemma 5 for arbitrary subsets of grid points, which is what Component 2 needs.

Selected references

• J. Zhao and Z. Zhou, Pigeonhole Design: Balancing Sequential Experiments from an Online Matching Perspective, Management Science, 2025; arXiv:2201.12936v6, 2024. https://arxiv.org/abs/2201.12936
• M. Ajtai, J. Komlós and G. Tusnády, On optimal matchings, Combinatorica 4 (1984) 259–264. https://doi.org/10.1007/BF02579135
• N. Fournier and A. Guillin, On the rate of convergence in Wasserstein distance of the empirical measure, Probability Theory and Related Fields 162 (2015) 707–738. https://doi.org/10.1007/s00440-014-0583-7
• W. Hoeffding, Probability inequalities for sums of bounded random variables, Journal of the American Statistical Association 58 (1963) 13–30. https://doi.org/10.1080/01621459.1963.10500830

Motivation

Firms and public agencies often roll a new product, feature or policy out gradually: a few markets first, more later, every market eventually. When the rollout is itself an experiment, its schedule determines how precisely the effect of the treatment can be estimated. A staggered rollout is the common design in which treatment, once switched on in a unit, stays on. Examples are a ride-sharing feature introduced city by city, a hospital protocol adopted ward by ward, or a stepped-wedge clinical trial. The effect of such a treatment often builds up over several periods, so the analyst wants the instantaneous effect and the lagged ones. The design question is then: in which period should each unit adopt?

Xiong, Athey, Bayati and Imbens (arXiv:1911.03764, Management Science 2024) answer this question in closed form for the trace criterion. Their non-adaptive result, Theorem 3.1, says that an optimal rollout has a treated fraction that follows a fixed S-shaped curve in time, and that each adoption cohort balances the observed and latent covariates.

Setting

There are $N$ units and $T$ periods. A design assigns to every unit $i$ an adoption time $A_i\in[T]\cup\{\infty\}$, where $\infty$ means never treated. The treatment indicator is coded $z_{it}=+1$ if $A_i\le t$ and $z_{it}=-1$ otherwise. The outcome model, (3.1) of the paper, is

$$Y_{it}=\alpha_i+\beta_t+X_i^\top\theta_t+\tau_0z_{it}+\tau_1z_{i,t-1}+\dots+\tau_\ell z_{i,t-\ell}+u_i^\top v_t+\varepsilon_{it}.$$

It has unit effects $\alpha_i$, time effects $\beta_t$, observed covariates $X_i\in\mathbb R^{d_x}$ with period-specific slopes $\theta_t$, and $\ell\ge0$ known lags. Latent covariates $u_i\in\mathbb R^{d_u}$ load on random factors $v_t$ with covariance $\Sigma_v$, and the noise $\varepsilon_{it}$ has variance $\sigma_\varepsilon^2$. The effects $\tau=(\tau_0,\dots,\tau_\ell)$ are estimated by generalized least squares on periods $\ell+1,\dots,T$. The per-period error covariance is $\Psi=U\Sigma_vU^\top+\sigma_\varepsilon^2I_N$, and the errors are uncorrelated across periods.

By the Gauss–Markov (Aitken) theorem, which is proved on Prove2Me as VectorSpaceOpt.gauss_markov, the covariance of the GLS estimator is $(W^\top\Sigma_e^{-1}W)^{-1}$ for the regressor matrix $W$. The precision $\operatorname{Prec}(\hat\tau;A)$ is therefore a deterministic matrix, (EC.2.1)–(EC.2.2): the Schur complement of the nuisance block. The T-optimal design problem (3.2) is

$$\max_{A\in([T]\cup\{\infty\})^N}\ \operatorname{tr}\operatorname{Prec}(\hat\tau;A).$$

With $h=\lfloor\ell/2\rfloor$, the paper defines a vector $a^{(\ell)}=(M^{(\ell)})^{-1}b^{(\ell)}$ of length $\ell-h$ from explicit $M^{(\ell)}$, $b^{(\ell)}$ ((EC.1.2)–(EC.1.3)). The S-shaped fraction $\omega^*_{\ell,t}$ of (3.4) has five stages: it is $-1$ for $t\le h$, equals $a^{(\ell)}_{t-h}$ for $h<t\le\ell$, is linear, $-1+(2t-(\ell+1))/(T-\ell)$, for $\ell<t\le T-\ell$, is antisymmetric in the fourth stage, and is $1$ at the end.

Formalization targets

Goal: Theorem 3.1

Assume $\sigma_\varepsilon^2>0$, $\Sigma_v=\sigma_\varepsilon^2I$, centred rows $\sum_iX_i=0$ and $\sum_iu_i=0$, $\sum_iX_iu_i^\top=0$, and, when $\ell>0$, $T>(\ell^3+13\ell^2+7\ell+3)/(8\ell)$. If

$$A\in\mathbb A_\ell\cap\mathbb A_X\cap\mathbb A_U,$$

then $\operatorname{tr}\operatorname{Prec}(\hat\tau;A')\le\operatorname{tr}\operatorname{Prec}(\hat\tau;A)$ for every design $A'$. Here $\mathbb A_\ell$ is the set of designs whose treated fraction at every $t$ is $(1+\omega^*_{\ell,t})/2$, and $\mathbb A_X$, $\mathbb A_U$ are the sets of designs whose adoption cohorts in periods $2,\dots,T$ have zero sum of $X_i$ (respectively $u_i$).

Milestones

1. Lemma EC.1.2 (corrected): $\operatorname{tr}\operatorname{Prec}(\hat\tau)=-\frac N{\sigma_\varepsilon^2}\sum_{j}\big(f_{j,1}+f_{j,X}+f_{j,U}\big)$. The first term depends only on the treated fractions, the second only on $X$, the third only on $U$.
2. The vector $a^{(\ell)}$ exists, and the resulting $\omega^*$ is feasible: nondecreasing, with values in $[-1,1]$, and the stage-2 values lie in $[-1,-1+(\ell+1)/(T-\ell)]$.
3. $\omega^*$ minimizes $f_1=\sum_jf_{j,1}$ over all nondecreasing paths in $[-1,1]^T$.
4. The covariate terms are nonnegative, and they vanish when $X^\top z_t$ and $U^\top z_t$ are constant in $t$.

Significance

The theorem turns a combinatorial design problem over $(T+1)^N$ schedules into explicit conditions. The fraction of treated units follows a universal curve that depends only on $(\ell,T)$, and the cohorts must balance the covariates. Two features are useful in practice. Without lags the optimal fraction grows linearly. With lags the optimal rollout starts and ends more slowly, which is what the extra stages capture. The paper uses these conditions to build practical designs (rounding the fractions to integers, stratifying on estimated latent factors) and an adaptive version.

The result is proved in the paper; this mission formalizes it. Formalizing it also repairs the record. Two statements in the printed proof are false as printed, and so is one display. The vector $b_\ell$ of Lemma EC.1.2 has an off-by-one, and its observed-covariate term is correct only when $X^\top X=N\,I$. The linear term of (EC.2.6) is summed over the wrong range. The mission states the corrected versions, each with a counterexample to the printed one. Theorem 3.1 itself is true as printed. Planning checked it by brute force on small instances, and checked the optimality of $\omega^*$ with an exact certificate for $\ell\le12$. As far as a search of Prove2Me shows, nothing about panel experimental design or GLS precision has been formalized there before.

Difficulty

The objective is a trace of a Schur complement, and the design enters through $\pm1$ entries in a large regression. Expanding that objective into the separable form of Lemma EC.1.2 requires exact control of the projection onto the two-way fixed effects and the period-specific covariate slopes. This is where the printed lemma goes wrong. The relaxed problem in $\omega$ is a convex quadratic program with monotonicity and box constraints. Its solution is not linear when $\ell>0$: the stage-2 values solve a linear system whose feasibility (bounds and monotonicity) holds only above a threshold in $T$. That is the role of the cubic threshold. A natural first idea is that "spread adoption evenly" is optimal. It is optimal only for $\ell=0$; for $\ell>0$ the KKT system forces nonlinear stages at both ends. A second step must also be checked: the $\pm1$ designs, whose fractions satisfy $\mathbb A_\ell$ and whose cohorts satisfy the balance conditions, minimize all three terms at once.

Formalization scope

• Units and periods are Fin N and Fin T, with periods 0-based. A design is Fin N → WithTop (Fin T), with ⊤ for never treated, and $z_{it}=\pm1$.
• The precision trace is $\sum_{j=0}^{\ell}\min_\gamma$ of the $\Psi^{-1}$-weighted residual sum of squares of the lag-$j$ column regressed on the nuisance columns. These are all $N$ unit dummies, a time dummy for each observed period, and a slope on $X_i$ for each observed period. This equals the trace of (EC.2.2) wherever that formula is defined. It is defined for every design and every $X$, with no matrix inverse of a possibly singular matrix.
• The random model enters only through $\Sigma_e=I\otimes\Psi$. Assumptions 2.1, 2.2 and 3.1 are encoded by the design type and this covariance; $\sigma_\varepsilon^2>0$ is a hypothesis.
• The threshold is read as "if $\ell>0$ then $T>\dots$", because the expression is undefined at $\ell=0$.
• Lemma EC.1.2 is stated with $b_{\ell,t}=(T_\ell+1-2t)/T_\ell$, and its $X$-term is $\frac1N z^\top(P\otimes H_X)z$, with $H_X$ the projector onto the column space of $X$, written without a pseudo-inverse.

The comparison design $A'$ ranges over all designs, never over $\mathbb A_\ell$ or relaxed fractions. The coding is $\pm1$, not $0/1$. The estimator uses only periods $\ell+1,\dots,T$. Coding with $0/1$ or dropping the period-specific slopes would change $\omega^*$ and the covariate term, and would be a different theorem. The set $\mathbb A_{\mathrm{opt}}$ can be empty for some $N$, in which case the goal holds vacuously, exactly as the paper's "if" allows. A sanity file exhibits nonempty instances, at $\ell=0$ with $d_x=d_u=1$ and at $\ell=1$, and reproduces a value of Example EC.1.4.

A complete development needs least-squares projections, Schur complements, centring matrices, and a KKT certificate for a monotone-constrained quadratic program. The first three are reusable well beyond this mission. Contributions to any milestone are welcome.

Selected references

• R. Xiong, S. Athey, M. Bayati, G. Imbens, Optimal Experimental Design for Staggered Rollouts, Management Science, 2024; arXiv:1911.03764v6. https://arxiv.org/abs/1911.03764
• A. C. Atkinson, V. V. Fedorov, The design of experiments for discriminating between two rival models, Biometrika 62(1), 1975. https://doi.org/10.1093/biomet/62.1.57
• J. Bai, Inferential theory for factor models of large dimensions, Econometrica 71(1), 2003. https://doi.org/10.1111/1468-0262.00392
• M. A. Hussey, J. P. Hughes, Design and analysis of stepped wedge cluster randomized trials, Contemporary Clinical Trials 28(2), 2007. https://doi.org/10.1016/j.cct.2006.05.007

Motivation

Online platforms run randomized experiments on users who interact with each other. When the treatment of one user affects the outcome of another, comparing treated and control users is biased. A standard remedy is cluster-based randomization: users are grouped into clusters with few edges between them, and whole clusters are assigned to treatment or control (Ugander et al., KDD 2013). The price is variance: there are far fewer clusters than users, and clusters differ in size, so the Horvitz–Thompson estimator of the treatment effect can be very noisy.

Candogan, Chen and Niazadeh (SSRN 3852100; Management Science, 2024) observe that the experimenter controls more than the marginal treatment probability of each cluster: it also controls the correlation between the assignments of different clusters. Correlating assignments negatively, for example by treating exactly a fixed number of clusters, can reduce the variance. Since the potential outcomes are unknown, the paper designs the correlation against the worst case. This mission formalizes the paper's main asymptotic result: a simple family of designs, independent block randomization (IBR), is asymptotically as good as the best possible design whenever no single cluster dominates.

Setting

There are $n$ clusters. Cluster $i$ is assigned $Z_i\in\{0,1\}$ ($1$ = treatment), and every cluster is treated with the same known probability $q\in(0,\tfrac12]$. An experiment is a joint law $P$ of $Z=(Z_1,\dots,Z_n)$ with $\mathbb P[Z_i=1]=q$ for every $i$; the set of experiments is $\mathcal P_q$. Its correlation matrix is $\Sigma(P)_{ik}=(\mathbb E_P[Z_iZ_k]-q^2)/(q(1-q))$.

Each cluster has a size $w_i\ge0$. Nature chooses a vector $y$ in the box $\times_{i\in[n]}[0,w_i]$, and the worst-case variance of a matrix $S$ is $\max_{y\in\times_i[0,w_i]}y^\top Sy$. Under the paper's Assumption 2.1, the variance of the Horvitz–Thompson estimator reduces to this quantity (Lemma 2.1). The robust correlation-design problem (1) is

$$V^{\mathrm{OPT}}=\min_{P\in\mathcal P_q}\ \max_{y\in\times_{i\in[n]}[0,w_i]}y^\top\Sigma(P)\,y .$$

Relaxing $\Sigma(P)$ to an arbitrary correlation matrix (positive semidefinite, unit diagonal) gives the lower bound $V^{\mathrm{LB}}\le V^{\mathrm{OPT}}$ of (5).

An IBR experiment sorts the clusters by decreasing size, $w_1\ge w_2\ge\dots\ge w_n$, cuts the sorted list into consecutive blocks, and on each block of $m$ clusters independently treats a uniformly random set of $\lfloor qm\rfloor$ clusters with probability $p=\lceil qm\rceil-qm$, and of $\lceil qm\rceil$ clusters otherwise. The best IBR experiment is found by a dynamic program, the Bellman equation (4). Its worst-case variance is $V^{\mathrm{DP}}$. For an integer $k\ge1$, the $k$-partition IBR experiment first groups the clusters into consecutive groups of $k$ and partitions each group optimally; its worst-case variance is $V^k$. Finally, $f_q(k)$ denotes $V^{\mathrm{OPT}}$ for $k$ clusters of size one.

Formalization targets

Goal: Theorem 4.8

For every $q\in(0,\tfrac12]$ there is a constant $C(q)$ such that for all sorted sizes with $w_1>0$,

$$\frac{V^{\mathrm{DP}}-V^{\mathrm{LB}}}{V^{\mathrm{LB}}}\le C(q)\sqrt{\frac{w_1^2}{\sum_{i\in[n]}w_i^2}},$$

and along every sequence of instances with $w_1^2/\sum_iw_i^2\to0$, both this quotient tends to $0$ and $V^{\mathrm{DP}}/V^{\mathrm{OPT}}\to1$. The constant is not fixed: only its dependence on $q$ alone is asserted, so an improved constant does not invalidate the goal.

Milestones

In attack order:

• additivity of the IBR worst case over blocks (§3.2);
• Proposition 3.1, the optimal experiment for equal sizes, with $V^{\mathrm{OPT}}=\tfrac14\tfrac{n^2}{n-1}$, $\tfrac{n+1}{4}$ or $h^*+h^*(h^*-1)\sigma$, and $4V^{\mathrm{OPT}}/n\to1$;
• the Bellman equation (4);
• Lemma 4.1, $V^{\mathrm{LB}}\le V^{\mathrm{OPT}}$;
• Lemma 4.2, $\max\{w_1^2,\tfrac14\sum_iw_i^2\}\le V^{\mathrm{LB}}\le V^{\mathrm{OPT}}\le\sum_iw_i^2$;
• the bound (13), $V^k\le f_q(k)S_k$;
• Lemma 4.3, $V^k/V^{\mathrm{LB}}\le f_q(k)\big(\tfrac4k+\tfrac{k-1}{k}\tfrac{w_1^2}{V^{\mathrm{LB}}}\big)$;
• Theorem 4.5, the approximation ratio of the optimal IBR experiment;
• Lemma 4.7, $f_q(k)\le\tfrac14(k+C(q))$.

Significance

The result. Problem (1) optimizes over all joint laws of $n$ binary variables, an exponentially large polytope, and its optimal design is unknown in general. Theorem 4.8 says that nothing is lost asymptotically by restricting to IBR designs, which are computable in polynomial time and easy to implement and explain. The rate quantifies the loss through a single statistic of the cluster sizes. Example 4.1 of the paper shows that the no-dominating-cluster condition cannot be dropped. Theorem 4.5 adds a finite-$n$ guarantee: for instance, a $7/3$-approximation at $q=\tfrac12$, against a factor $4$ for independent assignment.

Formalizing it. The result is proved in the paper; none of it is machine-checked. A formal proof requires: the second-moment calculus of exchangeable Bernoulli laws behind Proposition 3.1; the identification of $V^{\mathrm{OPT}}$ for equal sizes; the lower bound of Lemma 4.2, obtained through a randomized outcome vector; and the analysis chain of §4. The development also checks the paper's informal steps, such as "the worst case of a group grows with its sizes" in the proof of (13), and the proof of Lemma 4.7, whose displayed equality in Appendix B.7 holds only as an inequality (the lemma is unaffected).

Difficulty

The obvious argument compares $V^{\mathrm{DP}}$ with $V^{\mathrm{OPT}}$ directly. It fails because $V^{\mathrm{OPT}}$ has no closed form beyond equal sizes, so the comparison has to go through the relaxation $V^{\mathrm{LB}}$, and the lower bound on $V^{\mathrm{LB}}$ must hold for every correlation matrix at once. A second obstacle is that the IBR worst case is a maximum of a quadratic over a box, and maxima do not add in general; additivity holds only because the blocks are independent. Finally, the asymptotic statement needs Proposition 3.1 with all its constants: the linear growth $f_q(k)\approx k/4$ is what makes the bound vanish, and that growth needs the exact optimal value for $k$ equal clusters, not just a bound.

Formalization scope

All objects are finite and real-valued; no measure theory is used. Clusters are Fin n, 0-based ($w_1$ is w ⟨0, _⟩). Assignments are Fin n → Bool. A law is a nonnegative real weight function on assignments with total mass one and marginals $q$. The worst case is a real sSup over the whole box (not only its vertices). $V^{\mathrm{OPT}}$ and $V^{\mathrm{LB}}$ are real sInfs of nonempty sets bounded below by $0$. Partitions into consecutive blocks are Mathlib's Composition n, and $V^{\mathrm{DP}}$ is the minimum over them. $V^k$ is the minimum over the compositions refining the $k$-grouping, which equals the paper's construction because the IBR worst case is additive. $f_q(k)$ is $V^{\mathrm{OPT}}$ of $k$ unit clusters, not a closed formula. Every statement assumes $0<q\le\tfrac12$ and $w\ge0$. The §4 statements also assume sizes sorted decreasingly (Remark 4.1) and $w_1>0$, which excludes only the all-zero instance where every ratio is $0/0$. Proposition 3.1 assumes $n\ge2$, and $k\ge1$ wherever $k$ appears.

The statements rule out the trivializing formalizations. The constant $C$ is chosen after $q$ and before $n$ and $w$, never per instance. $V^{\mathrm{DP}}$ ranges over consecutive blocks of the sorted clusters, not over arbitrary partitions. It is defined from the actual IBR law and its correlation matrix, not from a hand-written block matrix.

All definitions live in one file, CorrClusterExp.IBR.Setting. Reusable beyond this mission: the robust correlation-design problem, the exchangeable block law, and the IBR construction. Proofs of any milestone are welcome, as are Lemma 3.3 (the worst case of a block) and the high-multiplicity and logarithmic-partition results (Theorems 4.9, 4.10), which are not part of this mission.

Selected references

• O. Candogan, C. Chen, R. Niazadeh, Correlated Cluster-Based Randomized Experiments: Robust Variance Minimization, Management Science, 2024; author manuscript SSRN 3852100, version of March 22, 2023. https://ssrn.com/abstract=3852100
• J. Ugander, B. Karrer, L. Backstrom, J. Kleinberg, Graph Cluster Randomization: Network Exposure to Multiple Universes, KDD 2013. https://doi.org/10.1145/2487575.2487695
• D. G. Horvitz, D. J. Thompson, A Generalization of Sampling Without Replacement from a Finite Universe, Journal of the American Statistical Association 47, 1952. https://doi.org/10.1080/01621459.1952.10483446

Motivation

A simulation metamodel is a cheap surrogate for an expensive stochastic simulation: after running the simulation at a few design settings, one predicts the mean response at settings never simulated. Kriging, the interpolation method of geostatistics and of the design and analysis of computer experiments (DACE; Sacks, Welch, Mitchell and Wynn 1989; Santner, Williams and Notz 2003), treats the unknown response surface as a realization of a random field and predicts by the best linear predictor. Kriging was developed for deterministic computer codes, where an observation is the exact response. A stochastic simulation returns the response plus sampling noise, and the noise may be correlated across settings when the experimenter uses common random numbers (CRN).

Ankenman, Nelson and Staum, Stochastic Kriging for Simulation Metamodeling (Proc. 2008 Winter Simulation Conference, pp. 362–370), extend kriging to this setting. The mission formalizes §2 of that paper: the form of the optimal linear predictor, called stochastic kriging, and the formula for its mean squared error, which shows how much the simulation noise costs. The source used throughout is the WSC 2008 proceedings paper, not the later Operations Research 58(2) 2010 article, whose numbering differs.

Setting

Fix a probability space $(\Omega,\mathcal F,P)$ and a set of design settings $\mathbf x$. On replication $j=1,2,\dots$ at setting $\mathbf x$ the simulation outputs (display (3), with constant trend)

$$\mathcal Y_j(\mathbf x)=\beta_0+\mathsf M(\mathbf x)+\varepsilon_j(\mathbf x).$$

The constant $\beta_0$ is the overall mean. $\mathsf M$ is a mean-zero random field, the extrinsic uncertainty imposed by the modeller. $\varepsilon_j(\mathbf x)$ is mean-zero sampling noise, the intrinsic uncertainty of the simulation; its variance may depend on $\mathbf x$, and noises at different settings may be correlated (CRN). All of these have finite second moments, and the field is uncorrelated with the noise.

An experiment design runs $n_i\ge1$ replications at each of $k$ design points $\mathbf x_1,\dots,\mathbf x_k$. The data are the sample means (4)

$$\bar{\mathcal Y}(\mathbf x_i)=\frac1{n_i}\sum_{j=1}^{n_i}\mathcal Y_j(\mathbf x_i),\qquad\bar{\mathcal Y}=(\bar{\mathcal Y}(\mathbf x_1),\dots,\bar{\mathcal Y}(\mathbf x_k))^\top,$$

and the target at a point $\mathbf x_0$ is the noise-free response $\mathsf Y(\mathbf x_0)=\beta_0+\mathsf M(\mathbf x_0)$.

Write $\Sigma_{\mathsf M}(\mathbf x,\mathbf x')=\mathrm{Cov}[\mathsf M(\mathbf x),\mathsf M(\mathbf x')]$. Then $\Sigma_{\mathsf M}$ is the $k\times k$ matrix of these covariances at the design points, $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$ is the vector $(\mathrm{Cov}[\mathsf M(\mathbf x_0),\mathsf M(\mathbf x_i)])_{i}$, and $\Sigma_\varepsilon$ is the $k\times k$ covariance matrix of the averaged noises $\sum_{j=1}^{n_i}\varepsilon_j(\mathbf x_i)/n_i$. A linear predictor (5) is $\lambda_0+\lambda^\top\bar{\mathcal Y}$ with arbitrary weights $(\lambda_0,\lambda)\in\mathbb R\times\mathbb R^k$; its mean squared error is $\mathrm E[(\lambda_0+\lambda^\top\bar{\mathcal Y}-\mathsf Y(\mathbf x_0))^2]$. The optimal MSE $\mathrm{MSE}^\star$ is the infimum of this quantity over all weights.

Formalization targets

Goal: display (7)

There is a rule $\Xi$, assigning a positive definite matrix $\Xi(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ to every pair of positive definite matrices, such that in every model as above with $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$, and at every $\mathbf x_0$,

$$\begin{aligned}\mathrm{MSE}^\star&=\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\\&=\Big[\Sigma_{\mathsf M}(\mathbf x_0,\mathbf x_0)-\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\Sigma_{\mathsf M}^{-1}\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\Big]+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top\,\Xi\,\Sigma_{\mathsf M}(\mathbf x_0,\cdot).\end{aligned}$$

The bracket is the usual kriging MSE. The goal also states that $\mathrm{MSE}^\star$ strictly exceeds it whenever $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)\neq0$. The matrix $\Xi$ is existential and depends only on the two covariance matrices, not on $\mathbf x_0$.

Milestone: display (6)

If $\Sigma_{\mathsf M}+\Sigma_\varepsilon\succ0$, the stochastic kriging predictor

$$\widehat{\mathsf Y}(\mathbf x_0)=\beta_0+\Sigma_{\mathsf M}(\mathbf x_0,\cdot)^\top[\Sigma_{\mathsf M}+\Sigma_\varepsilon]^{-1}(\bar{\mathcal Y}-\beta_0\mathbf 1_k)$$

is of the form (5), attains the smallest MSE among all predictors of that form, and is the only one that does.

Further items (not milestones)

• (9)–(10): the two-point example, with closed forms for the predictor and its MSE in terms of $\tau^2$, $r_{12}$, $r_0$, $\mathsf V$, $\rho$ and $n$.
• (14): the equicorrelated example, $\mathrm{MSE}^\star=\tau^2\big(1-k r_0^2/(1+(k-1)r+\gamma/n)\big)$.
• §3, p. 365: under the Gaussian Assumption 1, $\widehat{\mathsf Y}(\mathbf x_0)=\mathrm E[\mathsf Y(\mathbf x_0)\mid\bar{\mathcal Y}]$ almost surely.

Significance

Display (7) quantifies the price of simulation noise. With exact observations the best achievable error is the kriging MSE; with noisy observations it is larger by a positive definite form in the cross-covariance vector. Every point $\mathbf x_0$ that is correlated with the design therefore loses accuracy. The two examples make this concrete. In (10) the MSE increases with the intrinsic correlation $\rho$, which is why the paper concludes that common random numbers, a standard variance-reduction device for comparing systems, do not help prediction. Display (14) is the baseline against which the paper measures the cost of estimating the noise variance. The conditional-expectation statement shows that under Gaussian assumptions stochastic kriging is optimal among all predictors, not only linear ones.

The results are classical in substance: they are the best-linear-prediction formulas of second-order random-field theory, applied to data whose covariance is $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The paper states them without proof ("we can show"). No machine-checked version exists. Formalizing them yields a reusable, measure-theoretic best-linear-predictor theorem for random vectors defined from a model rather than postulated, together with the Loewner-order fact behind the positive definiteness of $\Xi$.

Difficulty

The algebra is short; the difficulty is in deriving the second-moment structure from the model, not assuming it. The covariance of $\bar{\mathcal Y}$ must be computed from the sample-mean definition, using bilinearity of covariance over finite sums of square-integrable variables, and shown to equal $\Sigma_{\mathsf M}+\Sigma_\varepsilon$. The cross-covariance with $\mathsf Y(\mathbf x_0)$ must be shown to equal $\Sigma_{\mathsf M}(\mathbf x_0,\cdot)$. The MSE of an arbitrary affine predictor then has to be written as a quadratic in $(\lambda_0,\lambda)$, and its infimum identified.

A first idea is to take "a random vector with mean $\beta_0\mathbf 1$ and covariance $\Sigma$" as the hypothesis. That proves a different, more abstract statement and leaves the model's content unproved.

For the second line of (7), the positive definiteness of $\Xi$ does not follow from a scalar argument. It is a matrix inequality, $(A+B)^{-1}\prec A^{-1}$ for $A,B\succ0$, and it fails without $\Sigma_\varepsilon\succ0$: when $\Sigma_\varepsilon=0$ the two lines coincide.

Formalization scope

• The model. Settings form an arbitrary type X; design points are x : Fin k → X, so the paper's index $i$ is i.val + 1. Replication $j\ge1$ is index j - 1 of a sequence ε : ℕ → X → Ω → ℝ.
• Moments and covariances. Field and noise are in $L^2$ and have mean zero. Covariances are Mathlib's covariance. The optimal MSE is an infimum over all of $\mathbb R\times\mathbb R^k$, with no unbiasedness constraint and no sign constraint.
• Hypotheses not printed in §2. Three are added and disclosed:
  1. uncorrelatedness of $\mathsf M$ and $\varepsilon$, the second-order content of Assumption 1's "independent of $\mathsf M$", without which the covariance of $\bar{\mathcal Y}$ is not $\Sigma_{\mathsf M}+\Sigma_\varepsilon$;
  2. square integrability;
  3. $\Sigma_{\mathsf M}\succ0$ and $\Sigma_\varepsilon\succ0$ in the goal. Lean's matrix inverse returns $0$ for a singular matrix, and $\Xi\succ0$ is false when $\Sigma_\varepsilon$ is singular.
• Generality kept. No Gaussianity, no independence across replications and no "no CRN" restriction is imposed on (6) and (7).
• Examples. $|r_{12}|<1$ in (10) and $r\le1$ in (14) are added, since both are correlations.
• Non-trivializing. A positive semidefinite $\Xi$, or a $\Xi$ chosen after $\mathbf x_0$, would make the goal trivial or strictly weaker; the statement fixes $\Xi$ as a function of $(\Sigma_{\mathsf M},\Sigma_\varepsilon)$ before the model and the point.

Infrastructure needed and reusable beyond this mission:

• covariance of averages of $L^2$ random variables;
• the MSE of an affine predictor as a quadratic form;
• minimization of a positive definite quadratic;
• the Loewner antitonicity of the inverse;
• for the conditional-expectation item, Gaussian conditioning via uncorrelated-hence-independent residuals.

Contributions to any of these are welcome.

Selected references

• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Proceedings of the 2008 Winter Simulation Conference, IEEE, pp. 362–370, 2008. http://www.informs-sim.org/wsc08papers/042.pdf
• B. Ankenman, B. L. Nelson, J. Staum, Stochastic Kriging for Simulation Metamodeling, Operations Research 58(2), 371–382, 2010 (later journal version, not used here). https://doi.org/10.1287/opre.1090.0754
• J. Sacks, W. J. Welch, T. J. Mitchell, H. P. Wynn, Design and Analysis of Computer Experiments, Statistical Science 4(4), 409–423, 1989. https://doi.org/10.1214/ss/1177012413
• T. J. Santner, B. J. Williams, W. I. Notz, The Design and Analysis of Computer Experiments, Springer, 2003. https://doi.org/10.1007/978-1-4757-3799-8