Keywords: degenerate U-statistics; dependency graph; hypergraph theory; kernel tests.

1 Introduction↩︎

U-statistics are a broad class of statistical estimators that extend the concept of the sample mean. Instead of simply averaging the observations, they average a symmetric, measurable function (kernel) of \(k > 1\) arguments over all \(\binom{n}{k}\) possible subsets of a sample of size \(n\). To mitigate the computational burden, Incomplete U-statistics [@blom1976some] consider only a subset, referred to as a design, of these \(\binom{n}{k}\) elements. While in principle the design can be strategically chosen to minimize the variance of the estimator, for ease of implementation, a random selection—either with or without replacement—is commonly used, despite not being the most efficient approach (see [@lee1990u], ch. 4.3).

The theoretical properties of U-statistics have been extensively studied in both Statistics [@lee1990u] and Probability [@korolyuk2013theory]. In particular, non-degenerate U-statistics have a Gaussian limiting distribution [@hoeffding1948u], while degenerate U-statistics converge to non-standard asymptotic distributions ([@lee1990u], ch. 3.2.3). For degenerate second-order U-statistics, i.e., when \(k=2\), [@gregory1977large] showed that the limiting distribution is an infinite weighted sum of centered \(\chi^2\) random variables. [@weber1981incomplete] extends the analysis to second-order degenerate incomplete U-statistics, demonstrating that both normal and weighted chi-square limiting behaviors can occur depending on the choice of the design. The latter non-standard distribution is often considered intractable for practical purposes; this is because the weights depend on the eigenvalues of a kernel-specific integral equation and, except for a few “nice” choices of the kernel function ([@lee1990u], ch. 6.2), this equation cannot be solved analytically.

On the applications side, U-statistics have found widespread use across various fields. These range from classical statistical estimation and inference tasks: e.g., in nonparametric hypothesis testing [@bergsma2014consistent; @yao2018testing], robust statistics [@joly2016robust; @MinskerWei2020], bootstrap and resampling methods [@leucht2013dependent; @bastian2024testing], to modern machine learning applications: e.g., in empirical risk minimization [@clemenccon2008ranking; @chen2023distributed], supervised learning ensembles [@mentch2016quantifying; @peng2022rates] and kernel methods [@gretton2005measuring; @gretton2012kernel; @liu2016kernelized]. Second-order U-statistics are the most common and widely used in practice. This class includes estimators of variance, covariance, Kendall’s Tau, Spearman’s Rho, Wilcoxon statistics, Cramer-Von Mises statistics. Additionally, they encompass novel estimators of the Maximum Mean Discrepancy (MMD) and the Kernel Stein Discrepancy, which are kernel methods designed for two-sample testing [@gretton2012kernel] and goodness-of-fit testing [@liu2016kernelized], respectively (see [@schrab2022efficient] for a discussion). Higher-order U-statistics are also relevant in applications, as they include estimators of dependence measures such as distance covariance [@szekely2007measuring], sign covariance [@bergsma2014consistent], and the Hilbert–Schmidt Independence Criterion (HSIC) [@gretton2005measuring].

Given the extensive applicability of U-statistics, there has been growing interest in studying incomplete U-statistics, both from a statistical perspective as in [@chen2019randomized; @sturma2024testing; @leung2026berry] and also from an applied machine learning perspective as in [@papa2015sgd; @clemenccon2016scaling; @schrab2022efficient], but there also remain significant challenges. On the one hand, selecting a design that minimizes variance is highly desirable, but it is a notoriously difficult combinatorial problem ([@lee1990u], ch. 4.3.2). Due to this challenge, researchers generally resort to random designs [@papa2015sgd; @clemenccon2016scaling; @chen2019randomized; @leung2026berry], which are easier to implement and analyze, but come at the cost of a suboptimal variance. On the other hand, the problem of determining the limiting distribution remains open. For instance, in the case of second-order degenerate incomplete U-statistics, no study has yet characterized how the limiting distribution transitions between normal and weighted chi-square laws, depending on the growth rate of the size of a minimum variance design. Addressing this gap is crucial for gaining a deeper understanding of the statistical properties of incomplete U-statistics and ensuring the validity of inferential procedures. Moreover, these theoretical insights would be essential for guiding practical applications, providing clear criteria for when a Gaussian approximation can be used or when resampling methods are required to approximate the limiting distribution.

In this paper, we present an efficient design construction and its asymptotic framework, to address the computational and statistical challenges posed by incomplete U-statistics. Our main contributions can be structured into four distinct points:

1. We show that the dependency graph [@baldi1989normal] of an incomplete U-statistic is the line graph of the design hypergraph.

2. We derive a new Berry–Esseen bound for incomplete U-statistics of deterministic designs, valid under the weakest assumptions known in the literature, using Stein’s method [@chen2004normal]. Moreover, we show that when both the maximum degree of the design hypergraph and the order \(k\) of the U-statistic are \(O(\log^q(n)\) for some \(q>0\), the bound converges to zero, ensuring asymptotic normality, even in the degenerate case and when the order diverges. Interestingly, these results are obtained without relying on the Hoeffding decomposition [@Hoeffding1961].

3. We develop efficient algorithms for constructing incomplete U-statistics based on equireplicate designs ([@lee1990u], ch. 4.3.2), a subclass of deterministic designs that tightly controls the maximum degree of the design hypergraph. The algorithms run in linear time in the design size and achieve minimum variance when \(k=2\). Our construction for \(k>2\) may be of independent interest to the combinatorial design and hypergraph communities.

4. We validate our theoretical results through numerical experiments on kernel-based test statistics, specifically the MMD and the HSIC. In a real-data application, we showcase a permutation-free variant of the MMD test that achieves substantial computational gains while matching the power and type I error control of its permutation-based counterpart, built on the same incomplete U-statistic.

Section 2 reviews background material and sets the notation for the paper. Section 3 derives a Berry–Esseen bound for incomplete U-statistics of deterministic designs and establishes conditions for Gaussian limiting distributions, even in the degenerate case and when the order diverges. Section 4 presents efficient algorithms for constructing equireplicate designs. Section 5 validates our theoretical results through extensive simulations on kernel-based test statistics and a real-data application presenting a permutation-free variant of the MMD test. Section 6 concludes with a summary and directions for future research. All proofs and technical details corresponding to each section of the main paper are provided in the supplemental materials at the end, where sections are labeled with an “S”. The R code for this paper can be found at https://github.com/CaesarXVII/IUstat_equireplicate_designs.

Existing works addressing convergence rates to limiting distributions of incomplete U-statistics that are most relevant to our setting and contributions include: [@rinott1997coupling], which examine a Markov-type dependence framework with applications to incomplete U-statistics; [@chen2019randomized], which consider randomized incomplete U-statistics in high-dimensional settings; [@sturma2024testing], which extend [@chen2019randomized] to a mixed degenerate setting with applications to testing a null hypothesis defined by equality and inequality constraints; [@kim2024dimension], which propose a new class of incomplete U-statistics that enable valid inference regardless of how the dimension scales with the sample size; [@shao2025u] that develop higher-order approximations for the sampling distribution of studentized non-degenerate incomplete U-statistics; and [@leung2026berry], which analyzes incomplete U-statistics of random designs—generated via Bernoulli sampling—across different regimes relating the sample size \(n\) and the expected design size. Among these prior works, none has provided a comprehensive framework encompassing: (i) finite-sample results on the distance to normality that account for both degenerate and non-degenerate cases; (ii) asymptotic analyses allowing the order \(k\) to grow with \(n\) and the design size to increase superlinearly in \(n\); and (iii) efficient algorithms for constructing minimum-variance designs when \(k=2\). Additional related work is discussed in 7.

2 Background and Notation↩︎

In this section, we introduce the necessary background and notation used throughout the paper. We strive for clarity and consistency in our notation, aligning it with the standard conventions in each of the domains we explore, namely U-statistics, combinatorial designs, and hypergraph theory. In particular, our notation and terminology for U-statistics closely follow [@lee1990u], for combinatorial designs we refer to [@ColbournDinitz2006] and [@wallis2016introduction], and for hypergraphs we primarily follow [@bretto2013hypergraph].

2.1 Incomplete U-statistics↩︎

Let \(X_1, \ldots, X_n\) be i.i.d. random variables taking values in a measurable space \((T, \mathcal{F})\) with common distribution \(P\). Consider the couple \((V,\mathcal{B}_k)\), where \(V = \{1, \ldots, n \}\) is the index set of the random variables and \(\mathcal{B}_k = \{ S \subseteq V, k \in V \mid \vert S \vert = k \}\) the collection of all subsets of \(V\) of size \(k\). Then, an incomplete U-statistic of order \(k\) can be written as:

\[U_{n,D}^{(k)} = \frac{1}{|D|} \sum_{ \{i_1, \ldots, i_k \} \in D} h( X_{i_1}, \ldots, X_{i_k}), \label{eq:IUk}\tag{1}\] where \(h: T^k \rightarrow \mathbb{R}\) is a fixed measurable kernel function that is symmetric in its arguments, i.e., \(h\left(x_1, \ldots, x_k\right)=h\left(x_{i_1}, \ldots, x_{i_k}\right)\) for every permutation \(i_1, \ldots, i_k\) of \(\{1, \ldots, k\}\) and \(D\subseteq \mathcal{B}_k\) represents the design of the incomplete U-statistics of size \(|D|\). If \(D = \mathcal{B}_k\), we obtain the complete \(k\)th-order U-statistic with \(|D| = \binom{n}{k}\). Likewise, in the combinatorial designs literature [@ColbournDinitz2006], the previously defined couple \((V,\mathcal{B}_k)\) identifies the complete (or trivial) \(k\)-design that considers only blocks of size \(k\).

For \(c\in \{0,1,2,\ldots, k\}\), let \(\mathcal{J}_c = \{ S_1, S_2 \in D \mid |S_1 \cap S_2 | = c \}\) be the set of all pairs of elements in the design which have \(c\) indices in common and denote \(f_c = |\mathcal{J}_c|\) its cardinality. Then, we can write the variance of an incomplete U-statistic of order \(k\) as:

\[\operatorname{Var} U_{n,D}^{(k)} =|D|^{-2} \sum_{c=0}^k f_c \; \sigma_c^2,\label{eq:var95ustat}\tag{2}\] where \(\sigma_c^2=\operatorname{Cov}\left(h\left(S_1\right), h\left(S_2\right)\right)\) with \(S_1, S_2 \in \mathcal{J}_c\) and \(h(S)\) is shorthand for \(h\left(X_{i_1}, \ldots, X_{i_k}\right)\), where \(S=\left\{i_1, \ldots, i_k\right\} \in D\). Moreover, by Theorem 4 in [@lee1990u] (page 15), we have:

\[\frac{\sigma_{b}^2}{b} \leq \frac{\sigma_{c}^2}{c} \label{eq:ineq}\tag{3}\]

for \(0 < b \leq c \leq k\) and \(b \in \{0, 1, \ldots, k\}\). In this work, we assume that \(\sigma_{k}^2 < \infty\). Combined with inequality 3 and the identity \(\sum_{c=0}^k f_c = |D|^{2}\), this assumption guarantees that \(\operatorname{Var} U_{n,D}^{(k)}\) is finite. Moreover, still by inequality 3 , if \(\sigma_c^2 = 0\), then it follows that \(\sigma_1^2 = \cdots = \sigma_{b}^2 = 0\) and the U-statistic is called degenerate of order \(c\). In this work, we assume \(\sigma_k^2 > 0\), which ensures that \(U_{n,D}^{(k)}\) can be at most degenerate of order \(k-1\). More generally, the order of the degeneracy determines the asymptotic distribution of the complete U-statistic. This is because \(\sigma_c^2 = 0\), implies that the first \(c\) terms vanish in the Hoeffding decomposition [@Hoeffding1961], which is a representation of a U-statistic of order \(k\) with a sum of uncorrelated U-statistics of order \(1, \ldots, k\). For example, when \(k=2\) and \(\sigma_1^2 = 0\), we have a first order degeneracy and [@gregory1977large] proved that \(n \;( U_{n,\mathcal{B}_2}^{(2)} - E[U_{n,\mathcal{B}_2}^{(2)}])\) converges in distribution to an infinite weighted sum of centered \(\chi^2\) random variables. However, if \(D \subset \mathcal{B}_k\), the limiting behavior of \(U_{n,D}^{(2)}\) can vary and crucially depends on the choice of the design [@weber1981incomplete].

2.2 Equireplicate designs↩︎

The class of equireplicate designs is the main focus of our paper. Within this class, every index occurs in the same number of elements of the design, usually called blocks. Any \(r\)-equireplicate design \(D\) based on \(V\), where the positive integer \(r\) denotes the replication parameter, satisfies:

\[|D| \; k = n \; r. \label{eq:equirep}\tag{4}\]

When \(k = 2\), incomplete U-statistics based on equireplicate designs achieve minimum variance, as demonstrated by Theorem 1 on page 195 of [@lee1990u]. When \(k > 2\), additional conditions are needed to ensure minimum variance, such as balanced incomplete block designs (BIBDs) and cyclic designs (see [@lee1982incomplete], Examples 2 and 7, respectively). Further details are provided in 8.1.

The existence of equireplicate designs requires specific divisibility conditions. In general, given a nonzero integer \(a\) and an integer \(b\), we write \(a \mid b\) to indicate that \(a\) divides \(b\), and \(a \nmid b\) otherwise. For our design constructions, we use modular arithmetic and some basic group theory (see [@gallian2021contemporary], ch. 2 for an overview). We write \(b \pmod n\) to denote the remainder upon dividing \(b\) by \(n\), but set \(b\pmod n\vcentcolon=n\) if the remainder is zero; we also write \(\text{gcd}(b,n)\) as the greatest common divisor of \(b\) and \(n\). In Section 4, we present efficient algorithmic constructions of equireplicate designs that use these algebraic concepts. Since we relabeled 0 as \(n\) when calculating values \(\pmod n\), we consider \(\mathbb{Z}_n=\{1,2,\ldots, n\}\) as the group of integers modulo \(n\) with additive identity element \(n\). This notation was chosen to ensure that the elements of our designs are labeled from the index set \(V=\{1,\ldots, n\}\).

2.3 Hypergraphs and deterministic designs↩︎

A hypergraph \(H\) is a couple \((V, E)\), where \(V\) is the vertex set, and \(E \subseteq \mathcal{P}(V)\setminus\{\varnothing\}\) is the hyperedge set, with \(\mathcal{P}(V)\) denoting the power set of \(V\). Each hyperedge \(e \in E\) is a non-empty subset of \(V\) and the singleton \(e_{v} = \{v\}\) represents a loop for all \(v \in V\). We call \(H(v)\) the set of hyperedges containing the vertex \(v\) and define the degree as \(d(v)=|H(v)|\) for all \(v \in V\). The average degree of a hypergraph \(H\) is denoted by \(\Bar{d}(H)\) and its maximum degree by \(\Delta(H)\). If each vertex has the same degree, we say that the hypergraph \(H\) is \(r\)-regular, which implies that \(\Bar{d}(H) = \Delta(H) = r\).

When every hyperedge contains exactly \(k\) vertices, i.e., \(|e| = k\) for all \(e \in E\), the hypergraph is called \(k\)-uniform. We can extend the classical handshaking lemma—that relates the number of edges in a graph with the sum of the degrees—to \(k\)-uniform hypergraphs. Indeed, for any \(k\)-uniform hypergraph \(H = (V, E)\), it holds that:

\[|E| \; k = \sum_{v \in V} d(v)\;, \label{eq:handshaking95kuni}\tag{5}\]

because each hyperedge contributes exactly \(k\) incidences, one for each vertex it contains.

To derive the Berry-Esseen bound and the asymptotic properties in section 3, we rely on a particular line graph construction. For a given hypergraph \(H = (V, E)\), its line graph \(L(H)\) is the couple \((E,E^{\prime})\) where \(E^{\prime} \;=\;\bigl\{\{e_i,e_j\} \;\bigm|\; e_i, e_j \in E \;\text{and}\; |e_i \cap e_j| \;\neq\;0 \bigr\}\). Note that our definition automatically includes loops, i.e., self-edges. In contrast, alternative definitions of the line graph explicitly exclude this possibility by imposing the additional condition \(e_i \neq e_j\) for any element of \(E^{\prime}\). In any case, \(L(H)\) is a graph, with vertex set that coincides with the hyperedge set of \(H\) and where an edge connects two vertices if and only if the corresponding hyperedges in \(H\) have at least one vertex in common.

A correspondence exists between \(k\)-uniform hypergraphs and deterministic designs¹:

Because of this equivalence, any structural or theoretical insight gained in the hypergraph setting translates directly into properties of the corresponding designs and viceversa. Thus, for instance, the couple \((V, \mathcal{B}_k)\) can be interpreted either as the complete \(k\)-uniform hypergraph—commonly denoted by \(K_n^{(k)}\)—or as the complete \(k\)-design that underlies the complete \(k\)th-order U-statistic. More generally, given a deterministic design \(D\subseteq \mathcal{B}_k\), the couple \((V,D)\) defines the hypergraph of the design \(D\), which we denote by \(\mathcal{D} = (V,D)\). When \(D\subset \mathcal{B}_k\), \(\mathcal{D}\) is a sub-hypergraph of \(K_{n}^{(k)}\) and its hyperedge set \(D\) is the design of an incomplete U-statistic of order \(k\). By equation 5 , any deterministic design \(D\) defined over \(V\) must satisfy:

\[|D| \; k = n \; \Bar{d}(\mathcal{D})\;, \label{eq:deter}\tag{6}\]

where \(\Bar{d}(\mathcal{D})\) is the average degree of the hypergraph of the design \(D\). We rely on this result to express the size of a deterministic design as \(|D|=n\;\Bar{d}(\mathcal{D})/k\). Clearly, for an \(r\)-equireplicate design, \(\Bar{d}(\mathcal{D}) = r\) and equation (6 ) boils down to equation (4 ).

3 Normal Approximations of Incomplete U-statistics↩︎

An incomplete U-statistic can be seen as a classical mean computed in the space induced by the kernel function \(h\). Therefore, equation 1 can be expressed more concisely as:

\[U_{n,D}^{(k)} = \frac{1}{|D|} \sum^{}_{S \in D} h(S), \label{eq:IUk95red}\tag{7}\] where \(h(S)\) stays for the random variable \(h\left(X_{i_1}, \ldots, X_{i_k}\right)\) with \(S=\left\{i_1, \ldots, i_k\right\} \in D\).

The random variables in the set \(\left\{ h(S), \;S \in D \right\}\) are identically distributed since the data, \(X_1, \ldots, X_n\), are i.i.d., and \(h\) is a fixed, symmetric and measurable kernel function. Consequently, we denote \(E\left[h(S)\right] = \mu_{k}\) for all \(S \in D\) and, following the notation introduced in section 2.1, we can conclude that \(\operatorname{Var}[h(S)] = \sigma^{2}_{k}\) for all \(S \in D\). However, unlike in a standard i.i.d. mean estimation setting, the set may contain random variables which are dependent. The main insight that underpins all the results of this section is that we can tame this dependence by controlling \(\Delta(\mathcal{D})\), the maximum degree of the hypergraph of the design. This quantity has a simple interpretation: for a given design \(D\), \(\Delta(\mathcal{D})\) represents the index that appears most frequently in the design.

More specifically, in this section we characterize the dependence structure of an incomplete U-statistic through its dependency graph, which allows us to establish a Berry-Esseen bound that holds for all incomplete U-statistics of deterministic designs, even in the degenerate case. Leveraging this bound, we study the asymptotic properties of this class of statistics, when key parameters—such as the maximum degree \(\Delta(\mathcal{D})\) and the order \(k\)—are allowed to grow with \(n\). Overall, our results indicate that, among all deterministic designs of a given size, equireplicate ones should be preferred because they provide precise control over \(\Delta(\mathcal{D})\).

3.1 The dependency graph of incomplete U-statistics↩︎

The set \(\left\{ h(S), \;S \in D \right\}\) contains the random variables whose average defines \(U_{n,D}^{(k)}\), as shown previously in equation 7 . To characterize the dependence structure of this set, we define its dependency graph, a concept originally introduced in [@baldi1989normal] to derive normal approximations of distributions:

In a dependency graph, there is an edge between two random variables only if they are dependent (absence of an edge implies independence). In the case of U-statistics, we include an edge between two random variables \(h(S_1)\) and \(h(S_2)\) if and only if \(S_1\cap S_2\neq\emptyset\).

Note that this corresponds with the definition of \(L(\mathcal{D})\), given in section 2.3:

The maximum degree of the dependency graph \(\Delta(L(\mathcal{D}))\), bounds how many random variables each \(h(S)\) can be dependent with. Controlling this quantity is key for establishing a normal limiting distribution, as we demonstrate in the following section.

3.2 Berry-Esseen bounds for incomplete U-statistics↩︎

We begin by deriving an interpretable and tight upper bound for the maximum degree of the dependency graph \(\Delta(L(\mathcal{D}))\), expressed solely in terms of \(k\) and \(\Delta(\mathcal{D})\). This result is then combined with key findings from the literature on dependency graphs to establish a Berry–Esseen bound for incomplete U-statistics based on deterministic designs.

Intuitively, the maximum value of \(\Delta(L(\mathcal{D}))\) is attained when the block \(S\), which identifies the variable \(h(S)\), consists of indices that each appear with the highest possible frequency in the design; i.e., each has degree \(\Delta(\mathcal{D})\) in the hypergraph \(\mathcal{D}\). In this case, each of the \(k\) indices in \(S\) can contribute up to \((\Delta(\mathcal{D}) - 1)\) edges in the line graph (excluding the self-edge), yielding a total of at most \(k \;(\Delta(\mathcal{D}) - 1)\) dependencies for \(h(S)\). In contrast, the minimum value of \(\Delta(L(\mathcal{D}))\) is reached when only one index in the block \(S\) has degree \(\Delta(\mathcal{D})\), while the remaining \(k - 1\) indices do not appear in any other blocks. In this setting, \(h(S)\) is connected to exactly \((\Delta(\mathcal{D}) - 1)\) other variables in \(L(\mathcal{D})\). We formalize the upper and lower bounds—both of which are tight—in the following lemma.

With this foundation, we can leverage powerful existing results from the literature on dependency graphs. In particular, we rely on [@chen2004normal] that couples Stein’s method with a concentration inequality to derive normal approximations under local dependence. Combining their findings with Lemma 1, we obtain a Berry-Esseen bound valid for all incomplete U-statistics based on deterministic designs.

As either \(k\) or \(\Delta(\mathcal{D})\) increase, while all other quantities are held fixed, the distance to normality grows. A similar effect occurs when the gap between the average degree \(\bar{d}(\mathcal{D})\) and the maximum degree \(\Delta(\mathcal{D})\) widens, which corresponds to deterministic designs whose associated hypergraph exhibits an increasingly skewed degree distribution.

Our Berry-Esseen bound distinguishes itself from existing results in the U-statistics literature by being valid in both degenerate and non-degenerate cases. It holds under minimal moment conditions on the kernel function \(h\) and makes explicit the critical role played by the order \(k\) and the maximum degree \(\Delta(\mathcal{D})\), both of which must be controlled to prevent deviations from normality. In contrast, existing works on convergence rates to limiting distributions for incomplete U-statistics either exclude the degenerate case altogether [@shao2025u], or address it under assumptions that are less interpretable than ours or stronger. For instance, [@rinott1997coupling] study the more general setting of weighted U-statistics but address the degenerate case only when \(k=2\), under comparatively less transparent assumptions. Likewise, [@kim2024dimension] focus on the second-order case and require the kernel to have a finite fourth moment to establish a Berry–Esseen bound for their cross U-statistics. [@chen2019randomized] derive Gaussian approximations valid for arbitrary order \(k\) and in high dimensions, but their results still rely on bounded polynomial moments of degree at least four and also [@sturma2024testing], who extend [@chen2019randomized] to a mixed-degenerate setting more in line with our unified framework, assume the kernel to be sub-Weibull. Finally, [@leung2026berry] establishes Berry–Esseen bounds for incomplete U-statistics under Bernoulli sampling—assuming a finite third moment of the kernel—however the analysis is limited to a fixed order \(k\).

In addition, our approach departs from traditional methods in the U-statistics literature by focusing on the dependency graph of the random variables \(\left\{ h(S), \;S \in D \right\}\), thus operating directly in the space induced by the kernel function \(h\). This perspective allows us to entirely bypass the Hoeffding decomposition, which is the standard analytical tool but is sensitive to degeneracy. Notably, when the U-statistic is degenerate of order \(k-1\), the first \(k-1\) terms of the Hoeffding decomposition vanish, leaving only the highest-order component. In contrast, the dependency graph is not affected by such degeneracy, as it encompasses all types of dependencies, including those of order \(k\).

Note that the presence of any type of degeneracy still impacts the variance of an incomplete U-statistic. To establish 8 , we relied on a lower bound for \(\operatorname{Var}U_{n,D}^{(k)}\) that remains valid even under extreme degeneracy, specifically when \(\sigma^2_{k-1} = 0\). A comprehensive discussion on the variance of incomplete U-statistics of deterministic designs—including tight upper and lower bounds that involve \(\Bar{d}(\mathcal{D})\) and \(\Delta(\mathcal{D})\)—can be found in 9.5.

Among deterministic designs, equireplicate ones offer tight control over \(\Delta(\mathcal{D})\). In an \(r\)-equireplicate design \(\mathcal{D}^{\dagger}\), we have \(\Delta(\mathcal{D}^{\dagger}) = r\) since every index appears exactly \(r\) times. In the next Corollary of Theorem 1, we establish a Berry-Esseen bound specific to all incomplete U-statistics based on equireplicate designs.

The result follows by substituting \(\Delta(\mathcal{D}) = \bar{d}(\mathcal{D}) = r\) in 8 , implying that there is no gap between the maximum and average degrees. Indeed, the hypergraph of an \(r\)-equireplicate design is \(r\)-regular, thus its degree distribution is a point mass at \(r\).

Moreover, Corollary 1 is valid for both degenerate and non-degenerate cases, and the bound in 9 is tighter for incomplete U-statistics based on \(r\)-equireplicate designs than the bound in 8 for U-statistics based on non-equireplicate deterministic designs of the same size (see Remark 2 for a detailed explanation).

3.3 Asymptotic results in the finite and infinite order regimes↩︎

Theorem 1 provides a Berry-Esseen bound for all incomplete U-statistics based on deterministic designs. To ensure convergence to a normal distribution, both the order \(k\) and the maximum degree \(\Delta(\mathcal{D})\) must grow slowly with the number of observations \(n\), as suggested by 8 . This requirement is formalized in the central limit theorem below, where \(k\) represents a sequence of natural numbers indexed by \(n\).

The proof, provided in 9, is based on the Berry-Esseen bound of Theorem 1. To our knowledge, Theorem 2 is the first asymptotic result in the U-statistics literature that simultaneously covers both the degenerate and non-degenerate cases, and is valid in both the finite-order regime (where \(k\) is fixed) and the infinite-order regime (where \(k\) grows with \(n\)), thus providing a unified theoretical framework with minimal assumptions.

However, a generic sequence of deterministic design \(D^{(k)}_n\) may not respect the condition \(\Delta(\mathcal{D}^{(k)}_n)=O(\log^q(n))\) required for Theorem 2 to hold. Specifically, this happens when a hypergraph \(\mathcal{D}\) has an unbalanced degree distribution, as shown in the following example.

On the other hand, for an equireplicate design, the parameter \(r\) determines the order of growth of the entire degree distribution of its associated hypergraph. This allows for considerable flexibility in the design construction as \(r\) can be chosen to grow at a desired rate with \(n\)—such as the \(O(\log^q(n))\) of Theorem 2—and even be fixed as \(n\) diverges.

In the next paragraphs, we compare Theorem 2 with existing results in the literature, highlighting its key features in both order regimes and further justifying the choice of equireplicate designs among deterministic ones. In our simulation studies of Section 5, we focus on the finite-order regime, which is relevant for kernel-based hypothesis testing.

Finite-order regime. When \(k\) is fixed, we omit the subscript and the superscript \(k\), referring simply to the kernel function \(h\) and the design \(D_n\). Thus, with minimal moment conditions on \(h\) and by requiring that the maximum number of dependencies in the set \(\left\{ h(S), \;S \in D_n \right\}\) grows at most with a logarithmic rate, Theorem 2 ensures asymptotic normality—even in the degenerate case—while preserving the classical computational advantages of incomplete U-statistics. This is because, any deterministic design \(D_n\) satisfies \(|D_n| = n \Bar{d}(\mathcal{D}_n)/k\) by equation 6 . Thus, since we assumed that \(\Delta(\mathcal{D}_n)=O(\log^q(n))\), then \(|D_n| = O(n \log^q(n))\) because \(\Bar{d}(\mathcal{D}_n) \leq \Delta(\mathcal{D}_n)\) by definition of the maximum degree.

Moreover, from the Berry–Esseen bound in 8 with \(p = 3\), it follows that under the assumptions of Theorem 2, the convergence rate to the standard normal distribution is \(n^{-1/2}\), up to logarithmic factors (see the proof of Theorem 2 for details). This rate matches the classical Berry–Esseen bound for the complete U-statistic in the finite-order non-degenerate case (see [@lee1990u], ch. 3.3 and references therein), again up to logarithmic terms. It is also only slightly slower than the \(o(n^{-1/2})\) rate established in [@korolyuk1989convergence] for the complete U-statistic in the finite-order degenerate case, where the limiting distribution is non-Gaussian. In that setting, the variance scales as \(n^{k/2}\) for the complete statistic, while for the incomplete one it scales as \(|D_n|^{1/2}\) (see Corollary 2). This last result confirms—and extends to deterministic designs—the findings of Remark 3.2 in [@chen2019randomized], which were stated for random designs.

Besides, existing results on normal approximations in the finite-order regime—such as [@chen2019randomized]—have to rely on bootstrap methods in practice to estimate the variance, which can increase the computational burden. On the contrary, our central limit theorem in 10 is readily applicable in the most extreme case of degeneracy—namely, of order \(k-1\)—for incomplete U-statistics based on equireplicate designs, as shown in Corollary 2. This setting is particularly relevant for hypothesis testing problems, including kernel-based tests, that we tackle in our simulation studies of Section 5.

The result follows from Theorem 2. This is because \(\Delta(\mathcal{D}_n) = r_n = O(\log^q(n))\) by assumption, and the variance simplifies to \(\operatorname{Var}U_{n,D^{\dagger}_n}^{(k)} = \sigma^2_k/|D^{\dagger}_n|\) by applying the lower bound in Lemma 2. The final expression then follows by substituting \(\operatorname{Var}U_{n,D^{\dagger}_n}^{(k)}\) in 10 , with \(|D^{\dagger}_n|\) determined by equation 4 . If the kernel is further assumed to be bounded, we can allow \(r_n = O(n^{1/4})\) in 11 , as explained in the proof of Theorem 5 of 9.

Proposition 2 shows that \(\sigma^2_k\) can be consistently estimated with the sample variance \(s^2_k\) calculated over the set \(\left\{ h(S), \;S \in D^{\perp}_n \right\}\), which contains i.i.d. random variables that do not have indices in common. More specifically, if \(k \mid n\), \(D^{\perp}_n\) is an equireplicate design with \(r=1\) and thus of size \(n/k\) by equation 4 ². For instance, fixing \(n=9\) and \(k=3\), we can consider the set \(\left\{ h\left(X_1, X_2, X_3\right),h\left(X_4, X_5, X_6\right), h\left(X_7, X_8, X_9\right) \right\}\) to calculate \(s^2_k\). The choice of this estimator for \(\sigma^2_k\) makes the implementation of 11 straightforward in practice, with negligible additional computational cost due to the estimation step.

Infinite-order regime. Recently, infinite-order U-statistics (IOUS) have attracted growing attention in the statistics and machine learning communities, particularly due to their relevance in uncertainty quantification for supervised learning ensembles (see e.g., [@mentch2016quantifying; @peng2022rates]). Theorem 2 contributes to this increasingly active and important area of research by providing the first result on incomplete U-statistics that considers a diverging order \(k\), even in the presence of degeneracy. The only related works on IOUS in the incomplete setting are [@SongChenKato2019], which develops distributional approximations for high-dimensional, non-degenerate IOUS and [@sturma2024testing], that extends [@SongChenKato2019] to a mixed-degenerate setting. However, their analyses rely on the Hoeffding decomposition, which becomes problematic in the IOUS regime. Specifically, when assuming \(\sigma^2_k < \infty\), the first-order term in the decomposition—known as the Hájek projection—vanishes as \(k \to \infty\), significantly complicating the analysis. This issue arises from inequality 3 , which implies \(k\sigma^2_1 < \sigma^2_k\). Hence, if \(\sigma^2_k\) is bounded and \(k\) diverges, we necessarily have \(\sigma^2_1 = O(k^{-1})\). As a result, the variance of the Hájek projection shrinks to zero, and controlling the moments of an increasing number of degenerate terms becomes challenging, as noted by [@SongChenKato2019].

Our approach avoids these complications as it does not rely on the Hoeffding decomposition and instead operates directly in the space induced by the kernel function \(h\). Consequently, it remains valid in both degenerate and non-degenerate cases without imposing assumptions on the order of the \(\sigma^2_c\) terms in equation 2 . The only additional requirement in the infinite-order regime is a logarithmic growth condition on \(\max\{k,\theta_k\}\), making our framework particularly interesting for the IOUS setting. Under the conditions of Theorem 2, the computational efficiency of incomplete U-statistics is preserved even in the infinite-order regime, by the same reasoning outlined for the finite-order case.

4 Efficient Construction of Equireplicate Designs↩︎

In this section, we highlight another major advantage of equireplicate designs: they can be constructed in linear time with respect to the design size.

For second-order incomplete U-statistics, equireplicate designs of any given size can be efficiently constructed and have minimum variance. Therefore, equireplicate designs should be preferred when \(k=2\), the most common and widely used case in practice.

When \(k > 2\), constructing equireplicate designs becomes challenging, as it relates to open problems in discrete mathematics. In Section 4.2 we present a construction which has still a linear computational complexity in the design size—when \(k\) is fixed—and that allows \(|D| = O(n^2)\) (see Remarks 9 and 8 for further details). This construction, based on cyclic permutations (see [@lee1982incomplete], Example 7), may be of independent interest for the combinatorial design and hypergraph community.

4.1 Construction of r-equireplicate designs when k = 2↩︎

For second-order incomplete U-statistics, an \(r\)-equireplicate design is a subset of \(\mathcal{B}_2\) of size \(|D|=nr/2\), such that each index appears in exactly \(r\) pairs. Our approach relies on a partition of \(\mathcal{B}_2\) into disjoint \(1\)-equireplicate designs when \(n\) is even (\(2 \mid n\)), and into disjoint \(2\)-equireplicate designs when \(n\) is odd (\(2 \nmid n\) and thus \(2 \mid r\) by equation 4 if the design exists). In both cases the whole partition can be generated sequentially, allowing us to obtain an \(r\)-equireplicate design by taking the union of a prescribed number of subsets from it. Algorithms 1 and 2 implement this procedure for even and odd \(n\), respectively. Theorem 3 establishes that these algorithms have linear computational complexity in the design size and that the resulting incomplete U-statistics achieve minimum variance. Additional details are provided in 10.1. That section also discusses the connections between our algorithms and factorizations of the complete graph \(K^{(2)}_n\).

4.2 Construction of r-equireplicate designs when k > 2↩︎

For incomplete U-statistics of order \(k\), an \(r\)-equireplicate design is a subset of \(\mathcal{B}_k\) of size \(|D|=nr/k\), such that each index appears in exactly \(r\) blocks. Our approach relies on a partial partition of \(\mathcal{B}_k\) into disjoint \(k\)-equireplicate designs (hence \(k \mid r\)). This partition can be generated sequentially, allowing us to obtain an \(r\)-equireplicate design by taking the union of a prescribed number of subsets from it. Algorithm 3 implements this procedure and Theorem 4 shows that, for any strictly increasing sequence of natural numbers \(\eta(\cdot)\), if \(n > 3\;\eta(k-1)\;\{\eta(k-1) - \eta(0)\}\) and \(r\in \{k,2k,\ldots,\phi(n)k \}\), then Algorithm 3 constructs \(r\)-equireplicate designs with linear computational complexity in the design size. In both the algorithm and theorem, we denote by \(\mathcal{C}_n = \{ a \in \mathbb{Z}_n| \; \text{gcd}(a,n) = 1 \}\) the set of coprimes of \(n\) and with \(\phi(n) = \mid \mathcal{C}_n \mid\) its cardinality, which is known as Euler’s totient function. Additional details are provided in 10.2, which also explores connections between our algorithm and factorizations of the complete \(k\)-uniform hypergraph \(K_n^{(k)}\).

5 Numerical Experiments↩︎

In this section, we empirically validate the theoretical results established in Section 3.3 regarding the asymptotic distribution of incomplete U-statistics of equireplicate designs in the finite order regime. We also investigate key aspects of our novel algorithms for constructing such equireplicate designs, introduced in Section 4, with particular focus on their variance-minimization property (when \(k=2\)) and linear computational complexity.

To this end, we conduct two sets of simulation studies on kernel-based testing methods: one on the two-sample test based on the unbiased MMD (uMMD) statistic, which is a second-order U-statistic (i.e., covering the \(k=2\) case), and the other on the independence test based on the HSIC [@gretton2007kernel], which is a fourth-order U-statistic (i.e., covering the \(k>2\) case). In addition, we illustrate our methodology on the widely used CIFAR-10 image classification dataset [@Krizhevsky09learningmultiple], which has frequently served as a benchmark for evaluating alternatives to the standard MMD two-sample test (see e.g., [@liu2020learning]). In this context, our novel approach provides a permutation-free version of the MMD test, that offers substantial computational advantages.

5.1 MMD experiments↩︎

In the first study, we perform a series of two-sample tests using incomplete versions of the uMMD test statistic—constructed via equireplicate designs—to assess departures from normality in their asymptotic distributions. Under the null \(H_0\)—i.e., when the two samples are drawn from the same distribution—the complete uMMD statistic (see eq. (6) in [@gretton2006kernel]) is a degenerate U-statistic and has a non-Gaussian limiting distribution. Under the alternative \(H_1\), it is non-degenerate and converges to a normal (see Theorem 8 in [@gretton2006kernel] for details).

More specifically, we simulate each time both samples \((X,Y)\) i.i.d. from a \(\mathcal{N}(0,1)\). We vary the common sample size \(n \in \{100, 200, 400,800,1600\}\). Each incomplete uMMD statistic is computed using a linear kernel and an \(r\)-equireplicate design generated by Algorithm 1. We select \(r \in \{1,\log(n),\log^2(n),\log^3(n), n/2, n-1\}\). We obtain the empirical distribution of the standardized uMMD statistic by repeating the experiment \(500\) times, standardizing the statistic in each replicate using its Monte Carlo estimate of the standard deviation. To quantify departures from normality, we compute the Kolmogorov-Smirnov (KS) distance between each empirical distribution and \(\mathcal{N}(0,1)\). We then repeat the entire procedure \(100\) times to produce a sampling distribution of KS distances. The left panel in Figure 4, reports \(95\%\) Monte Carlo confidence intervals (CI) for the KS distance (centered at the corresponding mean) and, for reference, shows \(Q_{0.975}\) and \(Q_{0.5}\), the \(97.5\%\) and \(50\%\) quantiles of the KS distance when the data are truly \(\mathcal{N}(0,1)\), respectively.

We observe that for \(r \in \{\log (n),\, \log^2 (n)\}\), the KS CIs lie below \(Q_{0.975}\) even at \(n=100\) and quickly contract toward \(Q_{0.5}\) as \(n\) increases. When \(r=\log^3(n)\), a larger sample size is required before the upper bound drops below \(Q_{0.975}\), but the decreasing trend is evident. In contrast, for larger values of \(r\) (e.g., \(r=n/2\)), the distance to normality remains roughly constant as \(n\) grows³. These findings validate Theorem 2 for \(k=2\), which predicts that for \(r=O(\log^q (n))\) (with fixed \(q\)), deviations from normality are small and diminish as \(n\) increases.

In practice, a Monte Carlo estimate of the uMMD standard deviation is not available. However, since the uMMD is degenerate under \(H_0\), we can use the variance estimator of Proposition 2. Figure 9 in 11 replicates Figure 4 but standardizes by \(s_2^2\) rather than by a Monte Carlo standard deviation. The results are virtually identical, thus using this estimator does not materially affect the KS CI as predicted by Corollary 2.

In the second study, we compare the variance of the uMMD statistic computed under equireplicate designs with that under random designs, to verify that with the equireplicate designs produced by Algorithm 1 we achieve minimum variance. To this end, we perform another series of two-sample tests, this time drawing the two samples from different distributions so that the uMMD statistic is non-degenerate.

More specifically, we simulate each time \(X \overset{iid}{\sim} \mathcal{N}(0,1)\) and \(Y \overset{iid}{\sim} \mathcal{N}(2,1)\). We vary \(n \in \{100, 200, 400,800,1200, 1600, 2000\}\) and select \(r \in \{1,\log(n),\log^2(n),\log^3(n)\}\). We follow the same procedure as in the first study to obtain Monte Carlo estimates of the variance of the uMMD statistic under both equireplicate and random designs. We then take their ratio (equireplicate/random) as a measure of relative efficiency and repeat the entire experiment \(100\) times to obtain a sampling distribution of variance ratios. The right panel of Figure 4 reports \(95\%\) Monte Carlo CI for these ratios, together with a reference line at \(1\), which corresponds to equal variance. We observe that the upper bounds of all CIs lie below \(1\), with a single exception at \(r=\log^{3}(n)\) when \(n=100\): in this case \(r \approx n-1\), so the statistic is effectively the complete uMMD and the variances coincide. As \(r\) increases, the efficiency gains diminish but appear to stabilize at a strictly positive level as \(n\) grows. These findings validate Theorem 3, confirming that the designs produced by Algorithm 1 achieve the minimum-variance property.

5.2 HSIC experiments↩︎

Following the same approach of the first MMD study, we conduct a series of independence tests using incomplete versions of the unbiased HSIC (uHSIC) statistic—constructed via equireplicate designs—to assess deviations from normality in their asymptotic distributions. Under \(H_0\)—i.e., when the two samples are independent—the complete uHSIC statistic (see the \(\mathrm{HSIC}_s(Z)\) definition in [@gretton2007kernel]) is degenerate of order \(1\) and has a non-Gaussian limiting distribution. Under \(H_1\), it is non-degenerate and converges to a normal (see Theorem 1 and 2 in [@gretton2007kernel] for details).

More specifically, we simulate each time both samples \((X,Y)\) i.i.d. from a \(\mathcal{N}(0,1)\). We choose \(n \in \{200, 400,800,1600\}\). The incomplete uHSIC statistic is computed from eq. (23) of [@schrab2025practical], accounting for kernel symmetrization, using a linear kernel and an \(r\)-equireplicate design constructed by Algorithm 3, with \(k=4\) and \(\eta(j)=2^j\) (see Remark 8). We select \(r \in \{1,\log(n),\log^2(n),\log^3(n)\}\)⁴. We obtain a sampling distribution of the KS distance between the empirical distribution of the standardized uHSIC statistic and \(\mathcal{N}(0,1)\) as described in the first MMD study. The left panel in Figure 5, reports \(95\%\) Monte Carlo CI for the KS distance and, for reference, shows \(Q_{0.975}\) as well as \(Q_{0.5}\).

We observe that for \(r \in \{\log(n),\, \log^2(n),\, \log^3(n)\}\), the KS CI lie below \(Q_{0.975}\) even at \(n=200\) and contract rapidly toward \(Q_{0.5}\) as \(n\) increases. The case \(r=\log^3(n)\) shows a slower decrease toward \(Q_{0.5}\), although the downward trend is evident. These findings support Theorem 2 for \(k>2\) (in particular \(k=4\)), which predicts that for \(r=O(\log^q(n))\) with fixed \(q\), deviations from normality remain small and vanish as \(n\) grows.

Following the same approach of the second MMD study, we compare the variance of the uHSIC statistic computed under equireplicate designs with that under random designs, to verify that incomplete U-statistics based on the equireplicate designs produced by Algorithm 3 do not lose efficiency with respect to the ones based on random designs. To this end, we perform another series of independence tests, this time generating dependent samples so that the uHSIC statistic is non-degenerate.

More specifically, we simulate each time \(X \overset{iid}{\sim} \mathcal{N}(0,1)\) and \(Y = 0.5 \sin{(X)} + \sqrt{3/4} E\) with \(E \overset{iid}{\sim} \mathcal{N}(0,1)\). We vary \(n \in \{200, 400,800, 1600\}\) and select the replication parameter \(r \in \{1,\log(n),\log^2(n),\log^3(n)\}\). We follow the same procedure as in the second MMD study to obtain a sampling distribution of the ratio between the Monte Carlo variance estimates of the uHSIC statistic under equireplicate and random designs. The right panel of Figure 5 reports \(95\%\) Monte Carlo CI for these ratios (equireplicate/random), together with a reference line at \(1\), which corresponds to equal variance.

There is no evidence of efficiency loss for equireplicate designs relative to random designs, since for every value of \(r\) the confidence intervals intersect values below \(1\). In contrast, we observe efficiency gains for \(r \in \{1, \log(n), \log^2(n)\}\), particularly at moderate sample sizes. Consistent with the second MMD study, these gains diminish as \(r\) increases. Based on these results, future research could investigate whether the designs produced by Algorithm 3 are indeed more efficient than random designs.

5.3 Real data example: CIFAR-10 dataset↩︎

We compare our permutation-free (PF) version of the MMD two-sample test against its standard permutation-based (PB) counterpart in terms of power, type I error, and runtime. Both methods rely on the same incomplete uMMD test statistic—built using equireplicate designs—to test for distributional differences between two balanced stratified samples of CIFAR-10 images. This dataset contains \(60000\) color images of size \(32 \times 32\) pixels, evenly distributed across \(10\) mutually exclusive classes—airplane, automobile, bird, cat, deer, dog, frog, horse, ship, and truck—with 6000 images per class.

More specifically, we sample without replacement \(X = \{X_{\text{cat}},\, X_{\text{deer}},\, X_{\text{ship}},\, X_{\text{truck}}\}\) and \(Y = \{Y_{\text{airplane}},\, Y_{\text{automobile}},\, Y_{\text{dog}},\, Y_{\text{horse}}\}\), where each class contributes equally, i.e. \(|X_{\text{class}}| = |Y_{\text{class}}| = n/4\). In total, we use 8 classes from the CIFAR-10 dataset—4 animals (cat, deer, dog, horse) and 4 vehicles (airplane, automobile, ship, truck)—split evenly between \(X\) and \(Y\), which yields a final population of \(N=48000\) images. We vary the common sample size \(n \in \{100, 200, 400,800,1600\}\). For both methods, each incomplete uMMD statistic is computed using a Gaussian kernel with bandwidth selected via the standard median heuristic, and an \(r\)-equireplicate design generated by Algorithm 1. For the PF test, we use the variance estimator of Proposition 2. We select \(r \in \{\log(n),\log^2(n)\}\), fix \(B=1000\) permutations for the PB test and choose \(\alpha = 0.05\). We repeat the experiment \(500\) times to obtain Monte Carlo estimates of power, type I error (with \(X\) and \(Y\) drawn from the same distribution), and the empirical distribution of runtimes. The left panel of Figure 6 shows \(95\%\) Monte Carlo CI for power, while the right panel reports the corresponding CI for type I error. Table 1 summarizes mean runtimes (in seconds) along with the half-width of the \(95\%\) CI.

We find no evidence of a loss of power for the PF MMD test relative to the PB version. For both methods, increasing \(r\) from \(\log(n)\) to \(\log^2(n)\) yields substantial power gains without sacrificing type I error control (see right panel of Figure 6). Across all \(n\) and \(r\), both tests maintain level \(\alpha\), though they are slightly conservative for smaller sample sizes—especially when \(r=\log^2(n)\)—with PF less conservative than PB. From \(n=800\) onward, both exhibit size \(\alpha\) behavior. These results support the validity of Corollary 2, of which the PF MMD test is a special case. In terms of computation, the speedup is substantial: the PF test runs about \(1000\times\) faster than PB, while retaining both power and validity. More broadly, these findings confirm the expected computational gains are proportional to \(B\), the number of permutations.

6 Conclusion↩︎

In this work, we introduced a novel characterization of the dependence structure of a U-statistic via its dependency graph. This perspective allowed us to derive a new Berry–Esseen bound that applies to all incomplete U-statistics based on deterministic designs, establishing conditions for Gaussian limiting distributions even in degenerate cases and when the order diverges. We further developed efficient algorithms for constructing incomplete U-statistics using equireplicate designs, a subclass of deterministic designs that, for second-order U-statistics, achieve minimum variance. All theoretical results have been validated through extensive numerical experiments. Finally, applying our framework to kernel-based testing, we proposed a permutation-free version of the MMD two-sample test, which—as shown by our real data example—delivers substantial computational gains while preserving both power and type I error control.

An important direction for future work is to extend our results to incomplete U-statistics with random designs, particularly for the case \(k>2\). In addition, future research may investigate whether \(\Delta(\mathcal{D}_n)\) can grow at larger polynomial rates and still maintain asymptotic normality. In a complementary direction, it would be interesting to move beyond the classical incomplete U-statistic framework by exploiting the full information contained in a partition of \(\mathcal{B}_2\) into disjoint equireplicate designs. Indeed, appropriate aggregation of the resulting statistics could allow one to retain asymptotic normality while leveraging all available pairwise interactions. In the specific case of the MMD two-sample test, this type of aggregation may lead to power improvements and could be compared with the current state-of-the-art permutation-free MMD test of [@shekhar2022permutation], which also relies on \(O(n^2)\) terms. Moreover, inspired by the work of [@janson2021clt_mdependent] on \(m\)-dependent processes, we conjecture that a Lindeberg-type condition may suffice to obtain the conclusions of Theorem 2 under only a finite second moment assumption on the kernel. Further work is also needed in the infinite-order regime: both to evaluate the applicability of existing variance estimators within our framework (see Remark 6) and to develop new ones tailored to the equireplicate design construction of Algorithm 3. More broadly, an interesting challenge is to extend our framework to dependent U-statistics [@dehling2002dependent], which frequently arise in applications, including the analysis of time-series and network data.

Acknowledgments and Funding↩︎

This work was supported in part by NSF grant SES-2150615. We also acknowledge the support of Purdue University, where both authors were affiliated at the beginning of this project. In preparing this work, OpenAI models were used to proofread specific paragraphs and refine code. The authors subsequently reviewed and edited all AI-assisted content, and take full responsibility for the integrity of the manuscript.

Disclosure Statement↩︎

The authors report there are no competing interests to declare.

7 Introduction↩︎

Related Work (extended version). In his seminal work, [@blom1976some] introduces incomplete U-statistics, analyzes their finite-sample and asymptotic variance, and establishes conditions for asymptotic normality in the non-degenerate case. He is also the first to suggest methods for constructing minimum variance designs, such as Latin squares and Graeco-Latin squares. [@brown1978reduced] derive asymptotic properties for incomplete U-statistics of second order, under an equireplicate design structure. They are the first to incorporate graph-theoretic language in their proofs. However, they do not consider the case of an odd replication parameter when the sample size \(n\) is even, and the parameter is not allowed to grow with \(n\). [@weber1981incomplete] and [@oneil1993asymptotic] show that in the degenerate case, the limiting behavior can be either standard or non-standard, depending on the choice of the design. [@lee1982incomplete] studies the problem of choosing minimum variance designs for incomplete U-statistics. [@janson1984asymptotic] provides a comprehensive treatment of the asymptotic distribution of incomplete U-statistics of order \(k\), considering both random and deterministic designs. However, he does not investigate specifically the class of equireplicate designs. In general, all previously mentioned works—being theoretical in nature—lack a practical procedure for constructing minimum variance designs. On the other hand, [@rempala2003incomplete] and [@kong2021design] propose asymptotically efficient incomplete U-statistics constructions, but they do not provide a finite sample analysis. Existing works addressing convergence rates to limiting distributions of incomplete U-statistics are: [@rinott1997coupling], which examine a Markov-type dependence framework with applications to incomplete U-statistics; [@chen2019randomized], which consider randomized incomplete U-statistics in high-dimensional settings; [@sturma2024testing] which extend [@chen2019randomized] to a mixed degenerate setting with applications to testing a null hypothesis defined by equality and inequality constraints; [@kim2024dimension], which propose a new class of incomplete U-statistics that enable valid inference regardless of how the dimension scales with the sample size; [@shao2025u] that develop higher-order approximations for the sampling distribution of studentized non-degenerate incomplete U-statistics; and [@leung2026berry], which analyzes incomplete U-statistics of random designs—generated via Bernoulli sampling—across different regimes relating the sample size \(n\) and the expected design size. Among these prior works, none has provided a comprehensive framework encompassing: (i) finite-sample results on the distance to normality that account for both degenerate and non-degenerate cases; (ii) asymptotic analyses allowing the order \(k\) to grow with \(n\) and the design size to increase superlinearly in \(n\); and (iii) efficient algorithms for constructing minimum-variance designs when \(k=2\).

8 Background and Notation↩︎

8.1 Minimum variance designs↩︎

Being equireplicate is sufficient for minimum variance designs only in the special case \(k=2\) (see Theorem 1 on page 195 of [@lee1990u]). For \(k>2\), additional structural constraints become necessary, as the intersection sizes among the subsets in the design can no longer be adequately controlled. This is why, in the combinatorial design literature, researchers impose additional conditions, beyond the equireplicate one, in order to obtain more balanced designs. For instance, any BIBD must satisfy the condition \(r(k-1) = \lambda (n-1)\), where the parameter \(\lambda\) represents the number of design elements in which each distinct pair of elements from \(V\) appears together. This additional constraint ensures balance not only at the level of individual elements of \(V\), but also among pairs of elements. However, explicit constructions of BIBDs are known only for specific cases e.g., when \(k=3\), we have Steiner triple systems (see ch. 12 in [@wallis2016introduction]) and certain specific families of designs described in [@sprott1954note]. Constructing BIBDs for general values of \(k\), on the other hand, remains a notoriously challenging problem. Another approach is to construct a design such that \(f_c = 0\) for all \(c \in \{2, \ldots, k\}\) i.e., to require that the intersection between any two distinct elements of \(D\) contains at most one element. If an equireplicate design satisfies this property, then by Theorem 2, page 196 in [@lee1990u], it attains minimum variance. In the main text, we refer to these type of designs as equireplicate and linear designs (see Remark 3 for further details). One way to practically build these designs is shown in Example 7 of [@lee1982incomplete]. Another particularly interesting way, was recently introduced in [@shao2025u]-Section 3.1-where the authors provide a novel approach to avoid unwanted overlaps among the blocks of the design.

8.2 Adjacency matrix of the line graph of an hypergraph↩︎

In this paragraph, we introduce the adjacency matrix of \(L(H)\), the line graph of the hypergraph \(H=(V,E)\). In particular, we denote it by \(A_{L(H)}\in \{0,1 \}^{|E| \times |E|}\) and define it as: \[(A_{L(H)})_{i,j} = \begin{cases} 1, & \text{if } |e_i \cap e_j| \neq 0, \\ 0, & \text{otherwise}. \end{cases}\]

To give a concrete example, consider \(K_{n}^{(k)}\) and its line graph \(L(K_{n}^{(k)})\). Then \(A_{L(K_{n}^{(k)})} \in \{0,1 \}^{\binom{n}{k} \times \binom{n}{k}}\) and the matrix indicates whether any given pairwise intersection among subsets of size \(k\) of \(V\) is empty or not. In the same way, we can consider the hypergraph of the design \(\mathcal{D} = (V,D)\), and denote its line graph with \(L(\mathcal{D})\). Then, its adjacency matrix \(A_{L(\mathcal{D})} \in \{0,1 \}^{m \times m}\) is a sub-matrix of \(A_{L(K_{n}^{(k)})}\) that indicates whether any given two elements of the design \(D\) share at least one index or not.

8.3 Equireplicate designs and \(k\)-uniform, \(r\)-regular hypergraphs↩︎

The above holds true because, as discussed in 2.3, the hyperedge set of a \(k\)-uniform hypergraph, being a subset of \(\mathcal{B}_k\) by definition, identifies uniquely a design \(D\) on the vertex set \(V\). If the \(k\)-uniform hypergraph is also \(r\)-regular, then it means that each vertex \(v \in V\) appears exactly \(r\) times in the collection of all hyperedges. Thus, the corresponding design \(D\) must be \(r\)-equireplicate since the vertices are the indices and occur in the same number of blocks i.e., elements of the design. The converse is true for the reverse reasoning. Because of this equivalence, equireplicate designs are also known as regular designs. Moreover, equation (4 ) can now be interpreted as an extension of the classical handshaking lemma—that relates the number of edges in a graph with the sum of the degrees—for \(k\)-uniform, \(r\)-regular hypergraphs.

9 Normal Approximations of Incomplete U-statistics↩︎

9.1 Proof of Proposition 1↩︎

9.2 Proof of Lemma 1↩︎

9.3 Proof of Theorem 1↩︎

9.4 Statement and proof of Theorem 5↩︎

9.5 Variance of incomplete U-statistics of deterministic designs↩︎

The dependency graph of an incomplete U-statistics i.e., \(L(\mathcal{D})\), encodes the presence or absence of dependence relationships between pairs of random variables in \(\left\{ h(S), \;S \in D \right\}\). Clearly, linear types of dependencies are also represented in \(L(\mathcal{D})\). Thus, for example, if there is no edge between two random variables \(h(S_1)\) and \(h(S_2)\), with \(S_1,S_2 \in D\), then \(\operatorname{Cov}[h((S_1),h(S_2)] = 0\). Obviously, the opposite direction does not hold, as the absence of a linear dependence does not imply independence. Following this line of reasoning, it is straightforward to conclude that an upper bound on \(\Delta(L(\mathcal{D}))\) is also an upper bound on the number of covariance terms in \(\operatorname{Var}U_{n,D}^{(k)}\). However, a U-statistic can be degenerate, which implies that \(\operatorname{Cov}[h((S_1),h(S_2)] = 0\) even if an edge connects \(h(S_1)\) and \(h(S_1)\) in \(L(\mathcal{D})\). Thus, Lemma 1 allows us to obtain only an upper bound for \(\operatorname{Var}U_{n,D}^{(k)}\). The lower bound represents a scenario of extreme degeneracy, where \(\sigma^2_{k-1} = 0\) which, by inequality (3 ), implies that \(\operatorname{Cov}[h((S_1),h(S_2)] = 0\) for all \(S_1,S_2 \in D\), with \(S_1 \neq S_2\). In this situation, excluding the case in which \(\operatorname{Var}U_{n,D}^{(k)} = 0\), only \(\operatorname{Var}[h(S)] = \sigma^{2}_{k} > 0\) and the number of these variance terms corresponds to the number of vertices of \(L(\mathcal{D})\). The next Lemma formalizes these results.

In the proof provided in the next paragraph, we begin by stacking the random variables in the set \(\left\{ h(S), \;S \in D \right\}\) to form the vector \(\boldsymbol{h(S)}\). We then leverage the interpretation of \(A_{L(\mathcal{D})}\)—the adjacency matrix of the line graph \(L(\mathcal{D})\) (see Section 8.2)—as an unweighted analogue of the variance-covariance matrix \(\Sigma = \operatorname{Var}[ \boldsymbol{h(S)}]\). This correspondence allows us to apply the upper bound established in Lemma 1. The upper bound is relatively loose for deterministic designs that are far from being equireplicate. However, the bound becomes tighter as the gap between the maximum degree \(\Delta(\mathcal{D})\) and the average degree \(\Bar{d}(\mathcal{D})\) decreases, for a fixed design size. Moreover, as long as \(k = o(\sqrt{n})\), the bound is asymptotically tight for designs \(D^{\diamond}\) that are both equireplicate and linear. When \(k \asymp \sqrt{n}\), the bound stays asymptotically tight for this class of designs, but now up to a constant (see the proof below for further details).

To conclude, we underline that the results derived in Lemma 2 and 4 encompass both degenerate and non-degenerate cases of incomplete U-statistics based on deterministic designs. Moreover, they imply that \(\operatorname{Var}U_{n,D}^{(k)}\) is \(O(\frac{1}{n})\) when \(k\) is fixed and \(\Delta(\mathcal{D})\) and \(\Bar{d}(\mathcal{D})\) grow at the same rate. This is akin to the variance of a classical mean estimator, even if the random variables in the set \(\left\{ h(S), \;S \in D \right\}\) are dependent.

9.6 Proof of Lemma 2↩︎

9.7 Statement and proof of Proposition 3↩︎

9.8 Statement and proof of Lemma 3↩︎

Proof sketch: the proof, provided below, leverages the extension of the classical handshaking lemma, discussed in Section 2.3, to \(k\)-uniform hypergraphs. Indeed, assuming that both designs have the same block size \(k\) and cardinality, by equations 4 and 6 , implies that \(\Bar{d}(\mathcal{D}) = r\). But then, intuitively, there must be an index which appears strictly more than \(r\) times in the non-equireplicate deterministic design, therefore forcing \(\Delta(\mathcal{D}) > r\).

9.9 Statement and proof of Lemma 4↩︎

9.10 Statement and proof of Corollary 3↩︎

9.11 Proof of Theorem 2↩︎

9.12 Proof of Proposition 2↩︎

9.13 Statement and proof of Corollary 4↩︎

10 Efficient Construction of Equireplicate Designs↩︎

10.1 Construction of r-equireplicate designs when k = 2↩︎

To facilitate our construction of equireplicate designs, we first introduce the new concept of an equireplicate partition, which partitions \(\mathcal{B}_2\) into disjoint equireplicate designs.

Note that if \(G_1,\ldots, G_{(n-1)/r}\) is an \(r\)-equireplicate partition, then the union of any \(q\) of the subsets yields a \(qr\)-equireplicate design. For \(n\) even, this implies that from a \(1\)-equireplicate partition, we can produce an \(r\)-equireplicate design for any \(r \in \{1,2,\ldots\, n-1\}\). For \(n\) odd, it is not possible to have a \(1\)-equireplicate design as the size would be \(n/2\), which is not an integer. Indeed, the only possible \(r\)-equireplicate designs in this case are for \(r\) even. Thus, for \(n\) odd, we can consider a \(2\)-equireplicate partition that enables the construction of an \(r\)-equireplicate design for any \(r \in \{2,4,6,\ldots\, n-1\}\), still by unioning the subsets of the partition.

Based on the above discussion, in Theorem 6 we construct a \(1\)-equireplicate partition for \(n\) even, and in Theorem 7 the construction a \(2\)-equireplicate partition for \(n\) odd. To aid understanding of our constructions, we provide a visual representation of the \(1\)-equireplicate partition in Figure 7 for \(n=6\) and of the \(2\)-equireplicate partition in Figure 8 when \(n=7\).

10.1.1 Statement and proof of Theorem 6↩︎

10.1.2 Statement and proof of Theorem 7↩︎

10.1.3 Proof of Theorem 3↩︎

10.2 Construction of r-equireplicate designs when k > 2↩︎

In principle, we can extend definition 2 to consider an \(r\)-equireplicate partition of \(\mathcal{B}_k\) and aim at building \(r\)-equireplicate designs by unioning the elements of either a \(1\)-equireplicate partition if \(k \mid n\) or of a \(k\)-equireplicate partition if \(k \nmid n\). However, when \(k \mid n\) and \(k>2\), even if the existence of a \(1\)-equireplicate partition is guaranteed by Baranyai’s theorem [@baranyai1975factorization], there is no known general sequential construction that would allow us to build an \(r\)-equireplicate design for any \(n\) and \(r \in \{1,\ldots,\binom{n-1}{k-1} \}\). Indeed, there exist only some efficient algorithmic constructions for \(k=3\) and \(k=4\) cases [@yan2022construction], but they impose additional divisibility conditions on \(n\). The situation becomes even more challenging when \(k \nmid n\). In this case, the problem is equivalent to constructing a \(k\)-factorization of the complete \(k\)-uniform hypergraph \(K_n^{(k)}\). However, even the existence of such a factorization is not guaranteed—let alone an efficient algorithm for its construction—as this remains an open problem in combinatorics (see e.g., [@bailey2010hamiltonian; @petecki2014cyclic] that discusses cyclic decompositions of \(K_n^{(k)}\)).

However, if we settle for a more attainable goal, we can still construct a \(r\)-equireplicate design in some settings, even when \(k>2\), via a partial equireplicate partition.

In the next theorem, we show that we can construct a partial \(k\)-equireplicate partition under some conditions on \(n\) and \(k\), and for any strictly increasing natural number valued sequence of choice. We denote with \(\mathcal{C}_n = \{ a \in \mathbb{Z}_n| \; \text{gcd}(a,n) = 1 \}\) the set of coprimes of \(n\) and with \(\phi(n) = \mid \mathcal{C}_n \mid\) its cardinality, which is known as Euler’s totient function.

10.2.1 Proof of Theorem 8↩︎

10.2.2 Proof of Theorem 4↩︎

11 Numerical Experiments↩︎

---

1. In 8.3, we discuss an analog of Observation 1 that applies specifically to equireplicate designs and \(k\)-uniform, \(r\)-regular hypergraphs.↩︎

2. If \(k \nmid n\), \(D^{\perp}_n\) would be of size \(\lfloor n\rfloor/k\) where \(\lfloor n\rfloor\) is the largest integer, smaller than \(n\), such that \(k \mid \lfloor n\rfloor\). To avoid unnecessary complications, we restrict ourselves to the case \(k \mid n\) in the main text.↩︎

3. Note that, as expected, when \(r=n-1\) the KS CI indicate a clear departure from normality. This occurs because we recover the complete uMMD, which under the null is non-Gaussian.↩︎

4. Note that \(r\) must be a multiple of \(k=4\). If it is not the case, we select the nearest multiple. Moreover, when \(r=1\), we use the design \(D_n^{\perp}\), described in Proposition 2.↩︎