U-Statistics: Theory and Applications

• U-statistics are symmetrized averages that compute a kernel over all distinct subsets of data to yield unbiased estimators, unifying measures like sample means and variances.
• They rely on the Hoeffding decomposition, splitting the statistic into linear and degenerate parts, with asymptotic behavior governed by the first nonvanishing projection.
• Modern approaches extend U-statistics to robust, computationally efficient methods — including incomplete and median-of-means variants — with applications in resampling, high-dimensional, sequential, and quantum settings.

U-statistics, often abbreviated “u-stats”, are symmetrized sample averages of a kernel over tuples of observations, introduced in Hoeffding’s sense as unbiased estimators of regular parameters such as expectations of pairwise or higher-order functions. For a symmetric measurable kernel $h:\mathcal X^m\to\mathbb R$, the classical target is $\theta=\mathbb E[h(X_1,\dots,X_m)]$, and the associated statistic averages $h$ over all distinct $m$-tuples from the sample. This construction underlies a large part of nonparametric inference because it unifies sample means, sample variance, Gini-type functionals, rank-based procedures, and many resampling-based estimators within a single projection-based framework (Joly et al., 2015, Fuchs et al., 2013).

1. Definition and canonical examples

For iid observations $X_1,\dots,X_n$, the order-$m$ U-statistic with symmetric kernel $h$ is

$$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$

where $I_n^m$ is the set of ordered $m$-tuples of distinct indices. In the degree-two case emphasized in several later developments,

$\theta=\mathbb E[h(X_1,\dots,X_m)]$0

The kernel is called degenerate in the order-two setting when $\theta=\mathbb E[h(X_1,\dots,X_m)]$1 for all $\theta=\mathbb E[h(X_1,\dots,X_m)]$2, and nondegenerate otherwise (Sharipov et al., 2015).

This general form specializes to many familiar statistics. The sample variance can be written as a degree-two U-statistic with kernel $\theta=\mathbb E[h(X_1,\dots,X_m)]$3, and Gini’s mean difference corresponds to $\theta=\mathbb E[h(X_1,\dots,X_m)]$4 (Sharipov et al., 2015, Wendler, 2010). In resampling-based classification error estimation, the complete “all learning–testing splits” estimator is itself a U-statistic once the learning/test asymmetry is removed by symmetrization of the kernel $\theta=\mathbb E[h(X_1,\dots,X_m)]$5 into $\theta=\mathbb E[h(X_1,\dots,X_m)]$6; the resulting statistic estimates the difference of unconditional error rates of two deterministic learning algorithms (Fuchs et al., 2013).

The complete U-statistic has two classical properties that remain central in later work. First, it is unbiased for the corresponding regular parameter. Second, when a parameter is representable through a symmetric kernel in the admissible nonparametric class considered by Hoeffding’s theory, the complete U-statistic has minimum variance among unbiased estimators; this same logic is used in modern work both for error-rate estimation and for unbiased estimation of the variance of that U-statistic itself (Fuchs et al., 2013).

2. Hoeffding decomposition and degeneracy regimes

The structural backbone of U-statistics is the Hoeffding decomposition. For degree two,

$\theta=\mathbb E[h(X_1,\dots,X_m)]$7

with

$\theta=\mathbb E[h(X_1,\dots,X_m)]$8

The term involving $\theta=\mathbb E[h(X_1,\dots,X_m)]$9 is the linear part, while $h$0 is degenerate. For nondegenerate kernels, asymptotic normality is governed by the linear projection; the degenerate remainder must be shown negligible after the relevant normalization (Sharipov et al., 2015).

In the general symmetric-kernel setting, degeneracy is indexed by the first nonvanishing Hoeffding projection. If $h$1 is $h$2-degenerate of order $h$3, then the first $h$4 projections vanish, and the leading stochastic scale becomes $h$5 rather than the generic $h$6. For bounded $h$7-canonical kernels, the Arcones–Giné benchmark recalled in later robust work yields deviation rates of order $h$8; in the canonical order-two case this is essentially $h$9 (Joly et al., 2015).

This dichotomy becomes sharper in continuous-monitoring asymptotics. In the nondegenerate regime, degree-two U-statistics admit Gaussian partial-sum approximations based on the first projection. In the degenerate regime, the leading approximation is not Gaussian mean behavior but a centered quadratic Gaussian chaos generated by the Hilbert–Schmidt spectrum of the kernel operator. This distinction drives the different time-uniform rates $m$0 and $m$1 in asymptotic confidence sequences (Cai et al., 14 May 2026).

A closely related noncommutative analogue appears in quantum U-statistics. There the kernel is a symmetric selfadjoint operator $m$2 acting on $m$3 tensor factors, the Hoeffding decomposition is formulated in $m$4, and the order of degeneracy again determines the normalization. If the first nonzero Hoeffding component is of order $m$5, then $m$6, and the normalized limit is expressed through Hermite-polynomial functionals of canonical variables in a CCR algebra rather than ordinary Gaussian chaoses (Guta et al., 2010).

3. Dependence, censoring, and nonstandard sampling schemes

Classical iid theory is only one branch of the subject. For dependent data, the central difficulty is that both the linear projection and the degenerate remainder must be controlled without destroying the dependence structure. One recent bootstrap construction for nondegenerate degree-two U-statistics from stationary absolutely regular or strongly mixing sequences avoids the usual “resample blocks of observations, then recompute the full U-statistic” strategy. Instead, it computes U-statistics within circular or nonoverlapping blocks, resamples those block summaries with replacement, and averages them. The method is asymptotically valid because the linear part coincides with that of the plug-in block bootstrap, while the degenerate remainder is shown negligible under mixing and moment assumptions; it also reduces the per-replication cost from recomputing an $m$7 statistic to averaging $m$8 precomputed block summaries (Sharipov et al., 2015).

A broader dependent-data theory exists for empirical U-processes, U-quantile processes, and generalized linear statistics under strong mixing or $m$9 near epoch dependence on an absolutely regular process. In that setting, the empirical U-process admits an almost sure Gaussian approximation, the U-quantile process follows from a generalized Bahadur representation, and functional CLTs and laws of the iterated logarithm become available for statistics such as Gini’s mean difference and robust scale estimators (Wendler, 2010).

Other nonstandard sampling schemes require more specialized constructions. For left-truncated and right-censored survival data, inverse-probability weighting modifies the kernel contribution by $X_1,\dots,X_n$0, where $X_1,\dots,X_n$1, $X_1,\dots,X_n$2, and $X_1,\dots,X_n$3 is a left-truncation-adjusted product-limit estimator for the censoring distribution. The resulting weighted U-statistics are $X_1,\dots,X_n$4-consistent, and their asymptotic variance splits into a first-order projection term plus an additional martingale term induced by estimation of the censoring weights (Sudheesh et al., 2021).

There are also settings in which the “sample points” entering the kernel are themselves dependent functionals. U-statistics of overlapping sample spacings are not standard Hoeffding U-statistics because normalized overlapping spacings form a dependent sequence; under the uniform null they admit an exponential representation through overlapping gamma blocks divided by a common normalization. The asymptotic null variance therefore includes both a long-run covariance term and a correction term generated by the common denominator (Singh et al., 2021). In weakly dependent time series, another nonclassical example is the order-two U-statistic built from blockwise local empirical moments $X_1,\dots,X_n$5, where both the triangular-array law and the interblock dependence must be handled simultaneously (Dehling et al., 2023).

4. Robustness, incomplete computation, and variance reduction

A central modern theme is robustness to heavy tails. Classical U-statistics concentrate well for bounded kernels, but can become unstable when $X_1,\dots,X_n$6 is heavy-tailed. A median-of-means analogue addresses this by partitioning the sample into $X_1,\dots,X_n$7 blocks, computing decoupled blockwise U-statistics across distinct blocks, and taking their median. For a symmetric kernel $X_1,\dots,X_n$8 that is $X_1,\dots,X_n$9-degenerate of order $m$0, finite variance alone suffices for

$m$1

with probability at least $m$2; for canonical kernels and finite $m$3-th moment with $m$4, the rate becomes $m$5 (Joly et al., 2015).

Computational considerations motivate incomplete U-statistics. A design-based construction, ICUDO, first divides the data into homogeneous groups, then uses an orthogonal array $m$6 to select group combinations, and finally samples one tuple within each chosen group-combination. Its mean-squared error satisfies

$m$7

where $m$8 contains only Hoeffding projection variances of orders exceeding the orthogonal-array strength $m$9. This means the leading $h$0 inflation term does not involve the low-order projections that dominate the variance of the full U-statistic. Under stronger smoothness, the remainder sharpens to $h$1, and in some cases asymptotic efficiency is possible even when $h$2 grows faster than $h$3 rather than faster than $h$4 (Kong et al., 2020).

A different computational role appears in importance-weighted variational inference. If a base objective or gradient estimator already requires $h$5 samples and one has a total budget $h$6, then averaging that base estimator over all overlapping subsets of size $h$7 yields a complete U-statistic estimator. Classical U-statistic variance theory implies

$h$8

so the complete U-statistic never has larger variance than averaging over $h$9 disjoint batches. Practical incomplete versions based on multiple random permutations of the $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$0 samples recover a fraction $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$1 of the complete-U variance reduction after only $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$2 structured passes (Burroni et al., 2023).

Label-efficient inference leads to another augmentation. In active inference for nondegenerate U-statistics with expensive labels, an augmented inverse probability weighted estimator combines a plug-in term based on predictions $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$3 with an IPW correction involving the queried labels. Its asymptotic variance is driven not by raw residuals $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$4, but by the first Hoeffding projection discrepancy $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$5. The variance-optimal sampling rule under a labeling budget is therefore proportional to

$$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$6

which generalizes the active-mean-estimation rule only when $$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$7 (Wang et al., 12 May 2026).

5. High-dimensional, adaptive, and sequential inference

Recent work has shifted from fixed-dimensional CLTs to pathwise approximations in high dimension. For vector-valued, order-two, nondegenerate U-statistics

$$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$8

a strong Gaussian approximation now couples the entire sequential process

$$U_n(h)=\frac{(n-m)!}{n!}\sum_{(i_1,\ldots,i_m)\in I_n^m} h(X_{i_1},\ldots,X_{i_m}),$$9

to a Gaussian partial-sum process $I_n^m$0 in Euclidean norm, with explicit error rates under mild $I_n^m$1 moments on the Hájek projection and only second moments on the degenerate remainder. This yields Brownian bridge approximations for U-statistic-based change-point processes and a self-normalized relevant testing procedure with a fully pivotal limit (Li et al., 11 Mar 2026).

A related but distinct high-dimensional program constructs, for each integer order $I_n^m$2, a U-statistic that unbiasedly estimates the $I_n^m$3-norm of a high-dimensional feature collection $I_n^m$4. Under the null, U-statistics of different finite orders are asymptotically normal and asymptotically independent, and they are also asymptotically independent of the corresponding maximum-type statistic, whose limit is extreme value. This separation allows adaptive combination of $I_n^m$5-values across dense-alternative and sparse-alternative regimes (He et al., 2018).

Sequential inference under continuous monitoring requires a different calibration. In the nondegenerate degree-two case, a leave-one-out jackknife estimator

$I_n^m$6

estimates the variance of the first projection strongly enough to support asymptotic anytime-valid confidence sequences

$I_n^m$7

In the degenerate case, the leading approximation is a centered quadratic Gaussian chaos, and the SAGE boundary calibrates its time-uniform excursions through the positive spectrum of the kernel operator (Cai et al., 14 May 2026).

High-dimensional change-point testing also exploits order-two distance U-statistics in a different way. Two processes based on $I_n^m$8-distance averages, one comparing within-segment pairwise distances and the other comparing cross-segment distances to a pooled baseline, share the same Brownian-bridge null behavior but separate under different alternatives. The first is especially sensitive to scale changes, the second to location changes, and their maximum yields a flexible detector in high dimensions (Boniece et al., 2022).

6. Applications and nonclassical generalizations

Applications of U-statistics are unusually broad because the kernel can encode estimation, testing, or risk. In binary classification, averaging the loss difference over all learning–testing splits of fixed size gives a complete U-statistic estimator of the unconditional error-rate difference between two deterministic algorithms. Its variance is itself a regular parameter of degree at most $I_n^m$9, which yields an unbiased U-statistic variance estimator whenever the total sample size satisfies $m$0; this in turn supports asymptotically exact tests of equality of error rates (Fuchs et al., 2013).

In clustering under heavy-tailed dissimilarities $m$1, the empirical risk

$m$2

may fail to concentrate uniformly over candidate partitions. Replacing it by the median-of-means U-statistic $m$3 restores finite-variance guarantees, and under a margin-type low-noise condition the excess clustering risk obeys faster rates of order $m$4 (Joly et al., 2015).

Goodness-of-fit testing by spacings provides another representative use. For overlapping $m$5-spacings, the order-two statistic

$m$6

is asymptotically normal under the uniform null and under $m$7-local alternatives. Within the class of such second-order tests, the kernel $m$8, corresponding to Gini’s mean squared difference of overlapping spacings, is asymptotically locally most powerful; its efficacy coincides with that of the Greenwood test based on overlapping spacings (Singh et al., 2021).

The concept also extends beyond classical probability spaces. Quantum U-statistics average a symmetric selfadjoint kernel over all $m$9-subsets of $\theta=\mathbb E[h(X_1,\dots,X_m)]$00 identical quantum systems prepared in $\theta=\mathbb E[h(X_1,\dots,X_m)]$01. Their normalized limits are Hermite-polynomial functionals of canonical CCR variables, with asymptotic normality in the nondegenerate case and quadratic-form limits for degenerate order-two kernels. The framework supports hypothesis testing beyond simple-versus-simple discrimination and asymptotic analysis of $\theta=\mathbb E[h(X_1,\dots,X_m)]$02-body interaction Hamiltonians in quantum metrology (Guta et al., 2010).

A plausible unifying implication is that U-statistics serve less as a single estimator than as a reusable inferential architecture: once a target can be expressed as the expectation of a symmetric kernel, the same projection logic, degeneracy taxonomy, and resampling or concentration machinery become available across iid, dependent, high-dimensional, sequential, censored, and even noncommutative settings.