A Berry-Esseen theorem for incomplete U-statistics with Bernoulli sampling

Suppose $X_{1},\dots,X_{n}\in M$ are independent and identically distributed random variables taking values in a given measurable space $(M,\mathcal{M})$. Hoeffding (1948) first introduced U-statistics, which, for a given degree $m\in\mathbb{N}$, are statistics of the form

where ${I_{n,m}}\equiv\{(i_{1},\dots,i_{m}):1\leq i_{1}<\cdots<i_{m}\leq n\}$ and $h:M^{m}\rightarrow\mathbb{R}$ is a given symmetric function in $m$ arguments, i.e.

for any permutation $(\pi_{1},\dots,\pi_{m})$ of the indices $(1,\dots,m)$. It generalizes the sample mean and encompasses many examples such as the sample variance; see Koroljuk and Borovskich (1994) for a comprehensive treatment of the classical theory on U-statistics. To simplify notation, for each $m$-tuple ${\bf i}=(i_{1},\dots,i_{m})\in I_{n,m}$, we will henceforth adopt the shorthands ${\bf X}_{\bf i}\equiv(X_{i_{1}},\dots,X_{i_{m}})$ and $h({\bf X}_{\bf i})\equiv h(X_{i_{1}},\dots,X_{i_{m}})$, as well as using $\mathcal{X}\equiv(X_{1},\dots,X_{n})$ to represent all the raw data inserted into $U_{n}$. Throughout this paper, we will assume

and that the single-argument projection function $g:M\rightarrow\mathbb{R}$ associated with the kernel $h$, defined by $g(x_{1})\equiv\mathbb{E}[h(x_{1},X_{1},\dots,X_{m-1})]$, is non-degenerate, i.e.

note that (1.3) necessarily implies $\sigma_{h}^{2}>0$ by Jensen’s inequality.

$U_{n}$ is complete, in a sense that it is an average of the kernel $h$ evaluated at all possible ordered subsamples of size $m$ from $\{1,\dots,n\}$. Blom (1976) suggested alternative constructions dubbed “incomplete U-statistics” which only average over an appropriate subset of ${I_{n,m}}$, with the intuition that the strong dependence among the summands in (1.1) should allow one to only use a well-chosen subset of them without losing too much statistical efficiency. In this work, we study one particular such construction, known as the incomplete U-statistic with Bernoulli sampling, which has the form

where the $Z_{\bf i}$’s are ${n\choose m}$ independent and identically distributed Bernoulli random variables with success probability

i.e., $P(Z_{\bf i}=1)=p=N{n\choose m}^{-1}$ and $P(Z_{\bf i}=0)=1-p=1-N{n\choose m}^{-1}$, and $\hat{N}$ is defined as

We call $N$ the computational budget because it is the expected number of kernel evaluations $\hat{N}$ that $U_{n,N}^{\prime}$ winds up summing over. In the sequel, we will also use $\mathcal{Z}\equiv(Z_{\bf i})_{{\bf i}\in{I_{n,m}}}$ to represent all these Bernoulli samplers, which are also assumed to be independent of the raw data $\mathcal{X}$.

It is an established result (Janson, 1984, Corollary 1) that, if

the following weak convergence holds for any degree $m\geq 2$ under (1.2)-(1.3):

where

Note that the asymptotic regime in (1.6) captures a common computationally efficient use case where the number of kernel evaluations being summed over by $U_{n,N}^{\prime}$ is in the same order as the raw sample size $n$, whereas the full-blown U-statistic $U_{n}$ entails a lot more summations, of the order $O(n^{m})$. The asymptotic normality of incomplete U-statistics such as (1.7) has found a resurgence of interest among the machine learning community, since Mentch and Hooker (2016) pointed out that it provides a computationally feasible framework for quantifying the uncertainty of ensemble predictions; see also Zhou et al. (2021), Peng et al. (2022), DiCiccio and Romano (2022), Xu et al. (2022) for some related recent works. While the convergent result (1.7) itself has been known for a long time, a Berry-Esseen (B-E) bound of a natural rate that characterizes the associated normal approximating accuracy is thus far missing from the literature.

In this paper, we offer the latter to fill such a gap. Capitalizing on the recent framework of Stein’s method for the so-called “Studentized nonlinear statistics” (Leung et al., 2024) and an exponential lower tail bound for non-negative U-statistics (Leung and Shao, 2023, Lemma 4.3), we establish B-E bounds for $U_{n,N}^{\prime}$ under minimalistic moment assumptions on the kernel $h$, which has the rate $1/\sqrt{\min(n,N)}$ when $n$ and $N$ are of the same order, the optimal rate we believe one could possibly achieve. In passing, we mention that the statistical community has also been developing B-E bounds for high-dimensional incomplete U-statistics (Song et al., 2019, Chen and Kato, 2019); our focus here is instead to provide a bound of the best rate for the one-dimensional convergent result in (1.7).

We first review the nature of the normal convergence in (1.7). To facilitate the discussion, we define the U-statistics

whose non-negative kernels $h^{2}$ and $|h|^{3}$ are induced by taking powers of the absolute value of the original kernel $h$. One can easily check that, as in the recent works on B-E bounds for incomplete U-statistics (Chen and Kato, 2019, Song et al., 2019, Peng et al., 2022, Sturma et al., 2022), $U_{n,N}^{\prime}$ admits the decomposition

where

Hence, up to multiplicative scaling factors, $U^{\prime}_{n,N}$ can be understood as a sum of the complete U-statistic $U_{n}$ and the term $B_{n}$. Under the assumptions (1.2) and (1.3), it is well-known that as $n$ tends to infinity, $\sqrt{n}U_{n}$ converges weakly to the limit $\mathcal{N}(0,m^{2}\sigma_{g}^{2})$ (Koroljuk and Borovskich, 1994, Ch.4). On the other hand, when conditioned on the raw data $\mathcal{X}$, $\sqrt{n}B_{n}$ should be approximately distributed as $\mathcal{N}(0,\alpha U_{h^{2}})$, for being a sum of independent mean-zero terms solely driven by the randomness of $\mathcal{Z}$, with variance $\frac{np}{N^{2}}\sum_{{\bf i}\in I_{n,m}}h^{2}(X_{\bf i})=\alpha U_{h^{2}}$; since

marginally, $\sqrt{n}B_{n}$ is expected to be approximately distributed as $\mathcal{N}(0,\alpha\sigma_{h}^{2})$ and also approximately independent of $U_{n}$. Altogether, in light of $N/\hat{N}\approx 1$, $\frac{\sqrt{n}U_{n,N}^{\prime}}{\sigma}$, being a sum of two approximately independent and normal random variables, is expected to be asymptotically normal with mean zero and variance

under the regime of (1.6), because $p\equiv N/{n\choose m}\longrightarrow 0$. Indeed, this is stated by the weak convergence in (1.7).

An extension of the preceding discussion also sheds light on what one should expect from a B-E bound that assesses the normal approximating accuracy of (1.7). Firstly, the complete U-statistic $U_{n}$ has a well-established B-E bound (Chen and Shao, 2007, Peng et al., 2022) that reads as

Moreover, ignoring the technical difficulty that may arise from division by zero when $U_{h^{2}}=0$ for the time being, one can write

so $\sqrt{N}B_{n}/\sqrt{U_{h^{2}}}$ is a sum of conditionally independent random variables $(Z_{\bf i}-p)h({\bf X}_{\bf i})/\sqrt{U_{h^{2}}N(1-p)}$ indexed by ${\bf i}\in{I_{n,m}}$, each having the conditional mean

and the conditional second and third absolute moments

where we have used the calculation $\mathbb{E}[|Z_{\bf i}|^{3}]=p(1-p)(1-2p+2p^{2})$. Since the second moment calculation above suggests that

the classical B-E bound (Berry, 1941, Esseen, 1942, Shevtsova, 2010) for a sum of independent variables with unit variance shall imply that

As discussed, since the asymptotic normality of $U_{n,N}^{\prime}$ under the regime (1.6) can be viewed as the distribution resulting from summing the two approximately independent and normal random variables $U_{n}$ and $B_{n}$, the two B-E bounds in (2.5) and (2.6) intuitively suggest that a “correct” B-E bound for the marginal cumulative distribution of $U_{n,N}^{\prime}$ should be of the rate $1/\sqrt{\min(N,n)}$ when $n\asymp N$, i.e. when there exists a positive constant $C>0$ such that

as a non-asymptotic analog of the regime (1.6).

This paper materializes the latter intuition into a Berry-Esseen theorem, which can be stated in terms of projection functions with respect to the non-negative kernel $h^{2}$ of the U-statistic $U_{h^{2}}$ in (2.1). For each $r=1,\dots,m$, first define the function $\tilde{\Psi}_{r}:M^{r}\rightarrow\mathbb{R}$ in $r$ arguments by

which is centered in the sense that $\mathbb{E}[\tilde{\Psi}_{r}]\equiv\mathbb{E}[\tilde{\Psi}_{r}(X_{1},\dots,X_{r})]=0$; specifically for $\tilde{\Psi}_{1}$, we also let $\Psi_{1}:M\rightarrow\mathbb{R}_{\geq 0}$ be its “uncentered version”

which must also be non-negative. Next, we can recursively define the Hoeffding projection kernels $\pi_{r}(h^{2}):M^{r}\rightarrow\mathbb{R}$ of $h^{2}$ as

and

These projection functions give rise to the famous Hoeffding decomposition (De la Pena and Giné, 1999, p.137) of the U-statistic $U_{h^{2}}$ that, upon re-centering and rescaling with $\sigma_{h}^{2}$, can be written as

where

are the leading terms in the decomposition. We note that each of the $\eta_{i}$’s has mean $\mathbb{E}[\eta_{i}]=m/n$, giving rise to the re-centering of $\sum_{i=1}^{n}\eta_{i}$ by $m$ in (2.9). We can now state our main B-E bound for $U_{n,N}^{\prime}$; recall $\mathbb{E}[\pi_{r}(h^{2})]\equiv\mathbb{E}[\pi_{r}(h^{2})(X_{1},\dots,X_{r})]$ and $\mathbb{E}[\Psi_{1}^{3/2}]\equiv\mathbb{E}[\Psi_{1}^{3/2}(X_{1})]=\mathbb{E}[(\Psi_{1}(X_{1}))^{3/2}]$, as per our convention in Section 1.1.

It is easy to see that each term on the right hand side of our B-E bound in Theorem 2.1 decays no slower than the rate $1/\sqrt{\min(N,n)}$ as $n,N\longrightarrow\infty$, when $n\asymp N$ as in (2.7); in particular, since the expression

is seen to be of order $O(n^{-2/3})$, the term $\mathcal{R}_{n,m,N,h}$ decays at the rate $n^{-1/2}$ under (2.7), whether under the minimalistic finite third moment assumption $\mathbb{E}[|h|^{3}]<\infty$ or the stronger finite fourth moment assumption $\mathbb{E}[h^{4}]<\infty$. Our result can also be contrasted with Peng et al. (2022, Theorem 5), which to our best knowledge is the first attempt in proving a B-E bound for the convergence in (1.7). However, Peng et al. (2022) had to assume six finite moments on the kernel $h$ and their bound exhibits a slower rate than $1/\sqrt{\min(N,n)}$ under (2.7). To further appreciate Theorem 2.1, the various parts of our bound could be interpreted as follows:

• •

  $\big{(}\frac{m}{n}(\frac{\sigma_{h}^{2}}{m\sigma_{g}^{2}}-1)\big{)}^{1/2}+\frac{\mathbb{E}[|g|^{3}]}{\sqrt{n}\sigma_{g}^{3}}$ stems from the complete U-statistic B-E bound in (2.5), which account for the “part-complete-U-statistic” behavior of $U_{n,N}^{\prime}$ coming from $U_{n}$; recall the decomposition in (2.2).
  

• •

  $\frac{\mathbb{E}[|h|^{3}](1-2p+2p^{2})}{\sigma_{h}^{3}(N(1-p))^{1/2}}+\exp(\frac{-[n/m]\sigma_{h}^{6}}{24(\mathbb{E}[|h|^{3}])^{2}})$ is analogous to the B-E bound in (2.6) for the conditional normal convergence of $\sqrt{N}B_{n}/\sqrt{U_{h^{2}}}$, which accounts for the “part-independent-sum” behavior of $U_{n,N}^{\prime}$ coming from $B_{n}$ as per the decomposition in (2.2). However, because Theorem 2.1 really concerns the convergence of the marginal distribution for $U_{n,N}^{\prime}$, the data-dependent ratio $U_{|h|^{3}}/U_{h^{2}}^{3/2}$ in (2.6) is now replaced by $\mathbb{E}[|h|^{3}]/\sigma_{h}^{3}$, a ratio of population moments that can be roughly thought of as its “average”; note that $\mathbb{E}[U_{|h|^{3}}]=\mathbb{E}[|h|^{3}]$ and $\mathbb{E}[U_{h^{2}}]=\sigma_{h}^{2}$. The extra term $\exp\big{(}\frac{-[n/m]\sigma_{h}^{6}}{24(\mathbb{E}[|h|^{3}])^{2}}\big{)}$ captures the fact that $U_{|h|^{3}}/U_{h^{2}}^{3/2}$ can never be too large due to the exponentially decaying lower tail behavior of $U_{h^{2}}$ as a non-negative U-statistic; see (3.20) below.

  The factor $\frac{1-2p+2p^{2}}{(1-p)^{1/2}}$ should also be understood with some care, because $p$ depends on $N$ and $n$ as per (1.5). Indeed, as $N$ approaches ${n\choose m}$, $p$ approaches $1$, which makes $\frac{1-2p+2p^{2}}{(1-p)^{1/2}}$ blows up and the whole bound become vacuous. However, our B-E bound in Theorem 2.1 is really geared for $n$ and $N$ having the same order as in (2.7), in which case $\lim_{n\rightarrow\infty}\frac{1-2p+2p^{2}}{(1-p)^{1/2}}=1$ as $\lim_{n\rightarrow\infty}p=0$. In fact, it is well-known that $U_{n,N}^{\prime}$ has different limiting distributions than the one considered in (1.7), when the relative growth rates of $n$ and $N$ are different from the prescription in (1.6) (Janson, 1984, Mentch and Hooker, 2016).
  

• •

  $\mathcal{R}_{n,m,N,h}$ stems from various “remainder terms” emerging in the proof of Theorem 2.1. For instance, one such remainder term is $D_{2}$ defined in (3.6) below, which corresponds to the error incurred by estimating $\sigma_{h}^{2}$ with $U_{h^{2}}$, as is implicit in (2.4). While the bound on $\mathcal{R}_{n,m,N,h}$ under the minimalistic third moment assumption $\mathbb{E}[|h|^{3}]<\infty$ appears to have a relatively complicated form as exhibited by (2.11), it comes from standard moment inequalities on the centered U-statistic $U_{h^{2}}-\sigma_{h}^{2}$; see Ferger (1996) and Koroljuk and Borovskich (1994, Chapter 2) for instance.

The proof of Theorem 2.1 adopts a technique called variable censoring, recently introduced in Leung et al. (2024) as a refinement of the variable truncation technique for proving B-E bounds with Stein’s method. For a given real-valued random variable $Y$, its censored version with respect to two constants $a,b\in\mathbb{R}\cup\{-\infty,\infty\}$, where $a\leq b$, is defined as

i.e. once $Y$ goes beyond the range $[a,b]$, its censored version $\bar{Y}$ can only be either $a$ or $b$ depending on which direction the value of $Y$ goes; that $a$ or/and $b$ are allowed to be infinity means that $Y$ may only be censored in one direction, or not be censored at all. This technique has the following apparent but very useful properties for establishing bounds:

We now explain the structure of Theorem 2.1’s proof. Throughout, we will use $\Omega$ to denote the entire underlying probability space. Moreover, we will define the following censored version of $\hat{N}$,

and the event

that only depends on the Bernoulli samplers in $\mathcal{Z}$; note that $\hat{N}_{c}=\hat{N}$ on $\mathcal{E}_{\mathcal{Z}}$. The first step of the proof is to establish, for any $z\in\mathbb{R}$, the bound

where

and

Following other existing works on B-E bounds for incomplete U-statistics, this can be developed with iterative arguments similar to ones pioneered by Chen and Kato (2019), and will be proved at the end of this section. Apparently, the remaining steps boil down to bounding each term on the right hand side of (2.14). The first term can be readily handled by the B-E bound for complete U-statistics in (2.5), where as the second term can be bounded by Bernstein’s inequality (van der Vaart and Wellner, 1996, Lemma 2.2.9, p.102) as in

where the last inequality comes from the fact that

and $\mathbb{E}[|h|^{3}]\geq\sigma_{h}^{3}$. The proof of Theorem 2.1 can then be completed by the following two lemmas on the other two terms in (2.14):

Collecting (2.5), (2.14) and (2.18), Lemmas 2.3 and 2.4 immediately render Theorem 2.1. The critical result is Lemma 2.3, the B-E bound for the term $T$ with the random threshold $\mathfrak{z}_{\mathcal{X}}$, which also depends on the data $\mathcal{X}$ via $U_{n}$ as suggested in (2.16). Most of the rest of this paper is devoted to proving Lemma 2.3, while the straightforward proof of Lemma 2.4 is included in Appendix A. We note that in (2.15), $\big{(}\frac{N-\hat{N}_{c}}{\hat{N}_{c}}\big{)}\frac{\sqrt{N}U_{n}}{\sigma_{h}\sqrt{1-p}}$ is a lower-order term while $\frac{N^{3/2}B_{n}}{\hat{N}_{c}\sigma_{h}}$ is the main term that drives the distribution of $T$; hence, it may be tempting to draw a parallel between Lemma 2.3 and the B-E bound for the conditional distribution $P(\sqrt{N}B_{n}/\sqrt{U_{h^{2}}}\leq z\mid\mathcal{X})$ in (2.6). Nevertheless, while (2.6) is a direct corollary of the classical result by Berry (1941) and Esseen (1942), the normalizer $\hat{N}_{c}\sigma_{h}$ for $N^{3/2}B_{n}$ is not the latter’s conditional standard deviation $N\sqrt{U_{h^{2}}}$ given $\mathcal{X}$. Moreover, unlike standard B-E bound problems, one has to handle a random threshold. These attributes make proving Lemma 2.3 non-trivially challenging.

Turns out, two ingredients that have only been recently formalized in the literature come in handy: the framework of proving B-E bounds for the so-called “Studentized nonlinear statistics” with Stein’s method (Leung et al., 2024), and an exponential lower tail bound for the U-statistic $U_{h^{2}}$ in (2.1). We note that analyzing the tail behavior of $U_{h^{2}}$ is quite an important aspect, as will become apparent later in our proof of Lemma 2.3 in Section 3. For another technical comparison with Peng et al. (2022, Theorem 5), we note that the Chebyshev’s inequality used therein (Peng et al., 2022, p.285) for controlling the tail of $U_{h^{2}}$ is too crude to produce a bound of the desired rate $1/\sqrt{\min(N,n)}$. In contrast, we exploit the non-negativity of $U_{h^{2}}$ and only control its lower tail, which decays exponentially in $n$ as per the latest result in Leung and Shao (2023, Lemma 4.3).

Before ending this section, we establish (2.14) as promised: