# U-Statistic Reduction: Higher-Order Accurate Risk Control and Statistical-
Computational Trade-Off, with Application to Network Method-of-Moments

Meijia Shao
Department of Statistics
The Ohio State University
and
Dong Xia
Department of Mathematics
The Hong Kong University of Science and Technology
and
Yuan Zhang
Department of Statistics
The Ohio State University

###### Abstract

U-statistics play central roles in many statistical learning tools but face
the haunting issue of scalability. Significant efforts have been devoted into
accelerating computation by U-statistic reduction. However, existing results
almost exclusively focus on power analysis, while little work addresses risk
control accuracy – comparatively, the latter requires distinct and much more
challenging techniques. In this paper, we establish the first statistical
inference procedure with provably higher-order accurate risk control for
incomplete U-statistics. The sharpness of our new result enables us to reveal
how risk control accuracy also trades off with speed for the first time in
literature, which complements the well-known variance-speed trade-off. Our
proposed general framework converts the long-standing challenge of formulating
accurate statistical inference procedures for many different designs into a
surprisingly routine task. This paper covers non-degenerate and degenerate
U-statistics, and network moments. We conducted comprehensive numerical
studies and observed results that validate our theory’s sharpness. Our method
also demonstrates effectiveness on real-world data applications.

Keywords: Nonparametrics, statistical learning, Edgeworth expansion, fast
computation.

##  1 Introduction

###  1.1 Background: U-statistics in statistical learning

A _U-statistic_ , associated with $X_{1},\ldots,X_{n}$ randomly sampled from a
general probability space and a degree-$r$ permutation-invariant kernel
function $h(x_{1},\ldots,x_{r})$ such that
$h(x_{1},\ldots,x_{r})=h(x_{\pi(1)},\ldots,x_{\pi(r)})$ for any bijection
$\pi:[1:r]\leftrightarrow[1:r]$, is defined as

$U_{n}:=\binom{n}{r}^{-1}\sum_{1\leq i_{1}<\ldots<i_{r}\leq
n}h(X_{i_{1}},\ldots,X_{i_{r}})=:\binom{n}{r}^{-1}\sum_{I_{r}\in{\cal
C}_{n}^{r}}h(X_{I_{r}}),$ (1)

where ${\cal C}_{n}^{k}:=\\{(i_{1},\ldots,i_{k}):1\leq i_{1}<\ldots<i_{k}\leq
n\\}$ and $X_{I_{k}}:=(X_{i_{1}},\ldots,X_{i_{k}})$ for $k\in[1:r]$.
U-statistics play central roles in many contemporary statistical learning
methods.

######  Example 1.1 (Example 1 of Kong and Zheng (2021)).

Test the symmetry of the distribution of $X\in\mathbb{R}$ by
$h(x_{1},x_{2},x_{3}):={\rm sign}(2x_{1}-x_{2}-x_{3})+{\rm
sign}(2x_{2}-x_{3}-x_{1})+{\rm sign}(2x_{3}-x_{1}-x_{2})$.

######  Example 1.2 (Bergsma-Dassios sign covariance (Bergsma and Dassios,
2014; Moon and Chen, 2022)).

To test the independence of $X\in\mathbb{S}_{X}$ and $Y\in\mathbb{S}_{Y}$,
where $\mathbb{S}_{X}$ and $\mathbb{S}_{Y}$ are Banach spaces equipped with
metrics $\rho_{X}$ and $\rho_{Y}$, define
$h\big{(}(x_{1},y_{1}),\ldots(x_{4},y_{4})\big{)}:=s_{X}(x_{i_{1}},\ldots,x_{i_{4}})s_{Y}(y_{i_{1}},\ldots,y_{i_{4}})$,
where $s_{X}(t_{1},\ldots,t_{4}):={\rm
sign}\big{\\{}\rho_{X}(t_{1},t_{2})+\rho_{X}(t_{3},t_{4})-\rho_{X}(t_{1},t_{3})-\rho_{X}(t_{2},t_{4})\big{\\}}$
and similarly define $s_{Y}$.

######  Example 1.3 (Treatment effect measurement (Rosenbaum, 2011; Zhao,
2019)).

Let $Y_{1},\ldots,Y_{n}$ denote the observed treated-minus-control matched
pair differences. Given integers $1\leq\underline{r}\leq\overline{r}\leq r$,
consider any $r$ observations $Y_{I_{r}}:=(Y_{i_{1}},\ldots,Y_{i_{r}})$,
sorted by their absolute values: $|Y_{I_{r}(1)}|\leq\ldots\leq|Y_{I_{r}(r)}|$.
Define
$h(Y_{I_{r}}):=\sum_{\ell=\underline{r}}^{\overline{r}}\mathbbm{1}_{[Y_{I_{r}(\ell)}>0]}$.
Here $I_{r}(\ell)$ denotes the $\ell$-th element of $I_{r}$.

###  1.2 Computational challenge: why we need U-statistic reduction

One major obstacle in the practical use of U-statistics lies in its heavy
computational cost. It takes $O(n^{r})$ time even to evaluate the point
estimator $U_{n}$. For example, $r=2$ for Maximum Mean Descrepancy (MMD,
assuming equal sample sizes) Gretton et al. (2012); Schrab et al. (2021) and
energy distance Székely et al. (2007); $r=3$ for testing distribution symmetry
(Example 1.1); $r=4$ for distance covariance (dCov) (Székely et al., 2007; Yao
et al., 2018) and sign covariance (Example 1.2); and the new U-statistic in
Rosenbaum (2011) has a user-selected degree that can go up to 20-ish (Tables 3
and 4 in Zhao (2019)). Efforts to alleviate the computational challenge can be
summarized into two mainstream directions. The first direction explores
shortcuts to fast-compute the complete U-statistic $U_{n}$. For example,
Rosenbaum (2011) shows that its proposed U-statistic can be evaluated in
$O(n\log n)$ time when the treatment effects are continuously distributed (no
ties); also, Huo and Székely (2016) and Chaudhuri and Hu (2019) show that dCov
can be computed in $O(n\log n)$ time when $X_{i}\in\mathbb{R}$; a similar
acceleration is available for sign covariance with scalar inputs (Even-Zohar
and Leng, 2021). Unfortunately, most such shortcuts require scalar input. For
example, the Bergsma-Dassios sign covariance (Example 1.2) with manifold-
valued functional trajectories as input (see Moon and Chen (2022), Simulation
II) cannot be accelerated by Heller and Heller (2016) or Even-Zohar and Leng
(2021) 111To see why, notice that the key step in Even-Zohar and Leng (2021)
discusses the joint order relationships between the coordinates of
$(x_{1},y_{1})$ and $(x_{2},y_{2})$. This enables stripping the absolute value
operations in $s_{X}(x_{[1:4]})s_{Y}(y_{[1:4]})$, where
$\rho_{X}(u,v)=\rho_{Y}(u,v):=|u-v|$. This trick is clearly inapplicable to
non-scalar inputs.. Moreover, with non-scalar $X_{i}$ inputs, it may even be
quite expensive to evaluate just one $h(X_{I_{r}})$ term. For example,
Chakraborty and Zhang (2021) computed the dissimilarity between two observed
time series of earthquake waves (Dau et al., 2018) – they directly computed an
$\ell_{2}$ pairwise distance, but we may want to synchronize them before
computing the distance (Srivastava et al., 2011; Strait et al., 2019).

The above discussion naturally leads us to the second acceleration strategy,
namely, _U-statistic reduction_ , by averaging over a much smaller set of
$r$-tuples in (1). Let

$\displaystyle{\cal J}_{n,\alpha}:=\big{(}I_{r}^{(1)},\ldots,I_{r}^{(|{\cal
J}_{n,\alpha}|)}\big{)}$ (2)

be a sequence of elements in ${\cal C}_{n}^{r}$ satisfying $|{\cal
J}_{n,\alpha}|\asymp n^{\alpha}$ for some $\alpha\in(1,r)$. We shall treat
${\cal J}_{n,\alpha}$ almost like a subset of ${\cal C}_{n}^{r}$, except that
${\cal J}_{n,\alpha}$ allows duplications. The _reduced U-statistic_ , also
called _incomplete U-statistic_ (Blom, 1976; Chen and Kato, 2019a), with
design ${\cal J}_{n,\alpha}$, is

$\displaystyle U_{J}:=$ $\displaystyle~{}|{\cal
J}_{n,\alpha}|^{-1}\sum_{I_{r}\in{\cal J}_{n,\alpha}}h(X_{I_{r}}).$ (3)

###  1.3 Two aspects of computational-statistical trade-off

There are two kinds of prices we must pay for computation reduction. First,
this reduction inflates $\mathrm{Var}(U_{J})$, which further determines: (i)
the confidence interval radius; and (ii) the minimum separation condition
$|\mu_{H_{a}}-\mu_{H_{0}}|$ for consistently testing $H_{0}:\mu=\mu_{H_{0}}$
versus $H_{a}:\mu=\mu_{H_{a}}$. Here $\mu:=\mathbb{E}U_{n}$. This aspect of
computational-statistical trade-off is easy to quantify thus well-understood.
The overwhelming majority of existing literature on U-statistic reduction
regards this aspect. The pioneering work Blom (1976) provides the earliest
systematic study, followed up by many works aiming at designing ${\cal
J}_{n,\alpha}$ smartly to minimize $\mathrm{Var}(U_{J})$ under a given
computational budget $O(n^{\alpha})$ (Lee, 1979, 1982, 2019; Rempala and
Wesolowski, 2003; Clémençon et al., 2016; Kong and Zheng, 2021; Dürre and
Paindaveine, 2021).

The second kind of price for speeding-up is the depreciation of _risk control
accuracy_ in statistical inference. The risk control accuracy refers to: (i)
$|\mathbb{P}(\textrm{true }\mu\in\textrm{CI})-(1-\beta)|$ (confidence
interval); and (ii) $|\mathbb{P}(\textrm{actual type I error rate})-\beta|$
(hypothesis testing), where the target confidence and significance levels are
$1-\beta$ and $\beta$, respectively. Technically speaking, it is governed by
the accuracy of approximating the distribution of a properly-designed
studentization of $U_{J}$. Most existing works only prove asymptotic normality
of $U_{J}$ (Brown and Kildea, 1978; Janson, 1984; Chen and Peng, 2021); while
to our best knowledge, Chen and Kato (2019a) is the only existing work that
provides finite-sample error bound result in this regard. However, revealing
the computational-statistical trade-off in risk control accuracy requires more
accurate characterization of the sampling distribution of studentized
incomplete U-statistic that is not yet available in existing literature222For
example, the “computational-statistical trade-off” in Remark 3.1 of Chen and
Kato (2019a) refers to the $\mathrm{Var}(U_{J})$ aspect; while their Berry-
Esseen bound on distribution approximation (Theorem 3.3) does not characterize
how risk control accuracy depreciates as computation reduces in the most
interesting regime $N\gg n$.. Our paper aims to fill this blank.

###  1.4 Our contributions

In this article, we present the first comprehensive study on _risk control
accuracy_ in statistical inference for incomplete U-statistics. Our results
cover non-degenerate, degenerate and noisy U-statistics. For non-degenerate
U-statistics, we established the first higher-order accurate distribution
approximation to incomplete U-statistics with general designs. This leads to
Cornish-Fisher CI’s and tests with higher-order accurate risk controls. Our
approach only requires the design to satisfy two natural, weak and easy-to-
verify assumptions. This converts the challenging task of analyzing many
different designs into routine calculations. Our method strictly complies with
the $O(n^{\alpha})$ computational budget and permits convenient
parallelization. For degenerate U-statistics, Weber (1981) and Chen and Kato
(2019a) observed a very interesting phenomenon that the incompleteness of
reduced U-statistics, apart from speeding up, also reinstates normality.
Pushing their results forward, we established the first finite-sample higher-
order accurate distribution approximation for the studentization of an
incomplete degenerate U-statistic. Thirdly, we applied our theory and methods
to incomplete network U-statistics, commonly known as _network moments_
(Bickel et al., 2011; Zhang and Xia, 2022).

Across all three scenarios, our method’s accuracy significantly improves over
the best existing results. The sharpness of our error bounds allows us to
reveal the trade-off between speed and the risk control accuracy of incomplete
U-statistics, for the first time in literature. For practitioners, our method
provides fast and easy-to-implement solutions with guidance on how much
computation would be required if the user aims at a certain risk control
accuracy goal under _finite sample size_. Our analysis discovers interesting
differences between the three scenarios. They are reflected in the variation
of practical guidance. For non-degenerate U-statistics, achieving higher-order
risk control accuracy only requires $\alpha>1$; and we can reduce the
computation cost from $O(n^{r})$ to $O(n^{2})$ almost “for free”: it does not
deteriorate risk control accuracy at all and only inflates
$\mathrm{Var}(U_{J})$ imperceptibly. In sharp contrast, for most network
moments (except the _edge_ motif), higher-order accurate inference requires
$\alpha>2$ even for densest networks; in presence of network sparsity, the
guidance will further vary accordingly.

The theory part of this paper is majorly different from complete U-statistic
literature Helmers (1991); Maesono (1997); Putter and van Zwet (1998) and
features several innovations. The incompleteness nature of our studied
U-statistic not only introduces new and complicated leading terms, but
moreover breaches the symmetry of remainder terms, making existing bounds
inapplicable. We addressed this using a different technique. A second example
is that the formulation of a succinct and weak Assumption 2, a highlighted
contribution of this paper, was also distilled from our theoretical
exploration. In the proof of Lemma 2.3, a key supporting result accompanying
Assumption 2, we handled the challenging dependency in some random sampling
schemes. On nonparametric network analysis, our Theorem 4.2 significantly
expands the classical analysis of triangles and V-shapes under sparse Erdös-
Renyi model in Gao and Lafferty (2017) to all motifs under general sparse
graph models. The proof features our novel martingale analysis.

In summary, what our paper caters to is not just one specific application or
data structure, but the infrastructural question of _risk control accuracy_ in
U-statistic reduction. The comprehensive theory & methods toolbox we developed
here provides an urgently needed (but currently missing) service to many
U-statistic-based learning methods who want to properly control their
inference risks while scaling up.

##  2 Reduction of non-degenerate noiseless U-statistics

###  2.1 Method and theory for general designs

Recall that ${\cal J}_{n,\alpha}\subseteq{\cal C}_{n}^{r}$ denotes the index
set for the reduced U-statistic $U_{J}$ defined in (3). Following the
convention (Lee, 2019; Chen and Kato, 2019a), we call ${\cal J}_{n,\alpha}$
the _design_. Throughout this subsection, ${\cal J}_{n,\alpha}$ is given and
fixed. For any $I_{k}\in{\cal C}_{n}^{k}$, define
$a_{n,\alpha;k}(I_{k}):=\big{|}\\{\widetilde{I}_{r}\in{\cal
J}_{n,\alpha}:I_{k}\subseteq\widetilde{I}_{r}\\}\big{|}$ to be the number of
elements in ${\cal J}_{n,\alpha}$ that contain $I_{k}$ as a subset. An
immediate fact is $\sum_{I_{k}\in{\cal
C}_{n}^{k}}a_{n,\alpha;k}(I_{k})={r\choose k}|{\cal J}_{n,\alpha}|$ for any
$k\in[1:r]$. By Han and Qian (2018), we have

$\displaystyle U_{J}=$ $\displaystyle~{}\mu+|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}(i)g_{1}(X_{i})+|{\cal
J}_{n,\alpha}|^{-1}\sum_{k=2}^{r}\sum_{I_{k}\in{\cal
C}_{n}^{k}}a_{n,\alpha;k}(I_{k})g_{k}(X_{I_{k}}),$ (4)

where $\mu:=\mathbb{E}[h(X_{[1:r]})]$, and define projection terms:
$g_{1}(X_{1}):=\mathbb{E}[h(X_{[1:r]})|X_{1}]-\mu$, for all $k\in[2:r]$,
$g_{k}(X_{[1:k]}):=\mathbb{E}[h(X_{[1:r]})|X_{[1:k]}]-\mu-\sum_{k^{\prime}=1}^{k-1}\sum_{I_{k^{\prime}}\in{\cal
C}_{[1:k]}^{k^{\prime}}}g_{k^{\prime}}(X_{I_{k^{\prime}}})$. All $g_{k}$ terms
are mutually uncorrelated. We call the U-statistic _non-degenerate_ if
$\xi_{1}^{2}:=\mathrm{Var}(g_{1}(X_{1}))\geq{\rm constant}>0$ and call it
_degenerate_ if $\xi_{1}=0$. Here, we focus on the non-degenerate case and
relegate the degenerate case to Section 3. In both Sections 2 and 3, we
observe $h(X_{I_{r}})$ without observational error – we call such U-statistics
_noiseless_. In Section 4, we study network moments as _noisy_ U-statistics,
where we only observe noise-contaminated $h(X_{I_{r}})$’s.

Next, we estimate $\mathrm{Var}(U_{J})$ for studentizing $U_{J}$. This
seemingly straightforward task (Han and Qian, 2018) now requires careful
thoughts due to incompleteness. Specifically, we need to choose between
estimating $\sigma_{J}^{2}:=\mathrm{Var}(U_{J})$ or
$\sigma_{J;1}^{2}:=\mathrm{Var}(|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}(i)g_{1}(X_{i}))$. This
question does not exist in the complete U-statistic setting where the two
quantities differ by $O(n^{-2})$ (Maesono, 1997; Zhang and Xia, 2022) , but
now, we have $|\sigma_{J}^{2}-\sigma_{J;1}^{2}|\asymp n^{-\alpha}$, which
cannot be directly ignored when $\alpha\in(1,2]$. In this paper, we choose to
estimate $\sigma_{J;1}$ for two reasons: formulation cleanness, and
computation convenience. Fortunately, the bias caused by using $\sigma_{J;1}$
instead of $\sigma_{J}$ to studentize $U_{J}$ will be properly accounted for
by the Edgeworth correction terms, see Remark 2.3. Now since
$\sigma_{J;1}^{2}=|{\cal
J}_{n,\alpha}|^{-2}\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}\xi_{1}^{2}$,
it suffices to estimate $\xi_{1}^{2}:=\mathrm{Var}(g_{1}(X_{1}))$ by

$\displaystyle\widetilde{\xi}_{1}^{2}:=$
$\displaystyle~{}n^{-\alpha}\sum_{i=1}^{n}\sum_{d=1}^{n^{\alpha-1}}h(X_{[i:d:(i+(r-1)d)]})h(X_{[i:(-d):(i-(r-1)d)]})-\widetilde{\mu}^{2},$
(5)

where
$\widetilde{\mu}^{2}:=n^{-\alpha}\sum_{i=1}^{n}\sum_{d=1}^{n^{\alpha-1}}h(X_{[i:d:(i+(r-1)d)]})h(X_{[(i+rd):d:(i+(2r-1)d)]}).$
Computing $\widetilde{\xi}_{1}^{2}$ costs $O(n^{\alpha})$. Denote the
shorthand $M_{\alpha}:=|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\asymp n^{\alpha-1}$.
Let $\\{B_{n}\\}_{n\geq 1}$ and $\\{b_{n}\\}_{n\geq 1}$ be sequences of random
variables and deterministic values, respectively. We write
$B_{n}={\widetilde{O}_{p}}(b_{n})$ if $\mathbb{P}(B_{n}\geq C\cdot
b_{n})=O(n^{-1})$ for some absolute constant $C>0$ and when $n$ is large
enough. The plug-in bias of $\widetilde{\xi}_{1}^{2}$ is characterized by the
following lemma.

######  Lemma 2.1.

Set $\alpha\in(1,2)$, we have

$\displaystyle\widetilde{\xi}_{1}^{2}-\xi_{1}^{2}=$
$\displaystyle~{}\frac{1}{n}\sum_{i=1}^{n}\big{\\{}g_{1}^{2}(X_{i})-\xi_{1}^{2}\big{\\}}+\frac{1}{n{M_{\alpha}}}\sum_{i=1}^{n}\sum_{d=1}^{{M_{\alpha}}}\sum_{\ell=1}^{r-1}g_{1}(X_{i})\Big{\\{}g_{2}(X_{i},X_{i+\ell
d})+g_{2}(X_{i},X_{i-\ell d})\Big{\\}}$
$\displaystyle+{\widetilde{O}_{p}}(n^{-\alpha/2}\log^{1/2}n).$ (6)

With Lemma 2.1, we are ready to formulate and decompose the studentization of
$U_{J}$.

$\displaystyle T_{J}:=$ $\displaystyle~{}\frac{U_{J}-\mu}{|{\cal
J}_{n,\alpha}|^{-1}\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{1/2}\cdot\widetilde{\xi}_{1}}=({\cal
T}_{1}+{\cal T}_{2})(1+{\cal T}_{3})^{-1/2}$ $\displaystyle=$
$\displaystyle~{}{\cal T}_{1}+{\cal T}_{2}-\frac{1}{2}{\cal T}_{1}{\cal
T}_{3}+{\widetilde{O}_{p}}\big{(}n^{-\alpha/2}\log n\big{)},$ (7)

where we define the following shorthand

$\displaystyle{\cal T}_{1}:=$
$\displaystyle~{}\frac{\sum_{i=1}^{n}a_{n,\alpha;1}(i)g_{1}(X_{i})}{\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{1/2}\xi_{1}},\quad{\cal
T}_{2}:=\frac{\sum_{k=2}^{r}\sum_{I_{k}\in{\cal
C}_{n}^{k}}a_{n,\alpha;k}(I_{k})g_{k}(X_{I_{k}})}{\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{1/2}\xi_{1}},$
$\displaystyle{\cal T}_{3}:=$
$\displaystyle~{}\sum_{i=1}^{n}\frac{g_{1}^{2}(X_{i})-\xi_{1}^{2}}{n\xi_{1}^{2}}+\frac{1}{n{M_{\alpha}}\xi_{1}^{2}}\sum_{i=1}^{n}\sum_{d=1}^{{M_{\alpha}}}\sum_{\ell=1}^{r-1}g_{1}(X_{i})\Big{\\{}g_{2}(X_{i},X_{i+\ell
d})+g_{2}(X_{i},X_{i-\ell d})\Big{\\}}.$ (8)

Next, we formulate the sampling distribution of $T_{J}$ with a fixed design
${\cal J}_{n,\alpha}$. We aim at addressing general ${\cal J}_{n,\alpha}$’s,
but of course it cannot be completely arbitrary. Fortunately, it turns out
that we only need two weak and natural regularity conditions on ${\cal
J}_{n,\alpha}$.

######  Assumption 1.

The design of ${\cal J}_{n,\alpha}$ is data-oblivious, namely,

$\displaystyle{\cal J}_{n,\alpha}\perp(X_{1},\ldots,X_{n}).$ (9)

For a deterministic ${\cal J}_{n,\alpha}$, (9) means ${\cal J}_{n,\alpha}$ is
designed without consulting the data $X_{[1:n]}$.

Assumption 1 is a very common assumption in U-statistic literature. We will
briefly discuss a recent but seminal paper Kong and Zheng (2021) on data-aware
designs in Section 7, but the quick gist is that being data-aware usually
requires (possibly expensive) additional computation. A second motivation for
data-oblivious designs comes from the analysis of network U-statistics, where
$X_{[1:n]}$ are latent (Zhang and Xia, 2022).

Our second assumption is motivated by the natural idea of avoiding bad designs
such as ${\cal J}_{n,\alpha}:=\big{(}[1:r],\ldots,[1:r]\big{)}$. Existing
literature suggests that a good design should be balanced, i.e.,
$g_{k}(X_{I_{k}})$’s with different $I_{k}$’s should have similar frequencies
in (4). Our next assumption formalizes this intuition.

######  Assumption 2.

Set $\alpha\in(1,r)$. It holds simultaneously for all $k\in[1:r]$ and
$I_{k}\in{\cal C}_{n}^{k}$ that

$\displaystyle
a_{n,\alpha;k}(I_{k})\in\begin{cases}[C_{2},C_{3}]n^{\alpha-k},&\textrm{ if
}k<\alpha,\\\ [0,C_{3}]\log n,&\textrm{ if }k=\alpha,\\\ [0,C_{3}],&\textrm{
if }k>\alpha,\end{cases}$ (10)

where $C_{2},C_{3}:0<C_{2}<C_{3}$ are universal constants. Moreover, for all
$k\in[1:r]$

$\displaystyle\sum_{I_{k}\in{\cal C}_{n}^{k}}a_{n,\alpha;k}^{2}(I_{k})\asymp$
$\displaystyle~{}n^{\alpha}\cdot n^{(\alpha-k)\vee 0}.$ (11)

One might expect to see more assumptions, but surprisingly, Assumptions 1 and
2 are all we need. This is one of the highlighted contributions of this paper.
These two assumptions are satisfied by our proposed deterministic “no-waste”
design in Section 2.2.1 and many popular randomized designs such as Chen and
Kato (2019a) with high probability. We shall show in Section 2.2.2 that four
typical random designs satisfy Assumption 2 with high probability.

Aside from Assumptions 1 and 2, another common assumption in U-statistic
literature is Cramér’s condition (Helmers, 1991; Maesono, 1997).

$\limsup_{t\to\infty}\big{|}\mathbb{E}[e^{{\mathbbm{i}}t\xi_{1}^{-1}\cdot
g_{1}(X_{1})}]\big{|}<1.$ (12)

The assumption (12) is undesirably restrictive and violated by important
applications, e.g., Example 1.1 with a discrete $X_{1}$ distribution. To waive
(12), inspired by Lahiri (1993) and Shao et al. (2022), we add to $T_{J}$ an
_artificial_ smoothing term $\delta_{J}\sim
N(0,\sigma_{\delta}^{2}=C_{\delta_{J}}\log n\cdot n^{-\alpha/2})$ independent
of $T_{J}$ with a large enough constant $C_{\delta_{J}}>0$. We will show that
$\delta_{J}$ waives (12) by smoothing the CDF of $T_{J}$ without altering the
distribution approximation formula and preserves the error bound333“Preserves
the error bound”: we can achieve the same error bound in approximating
$F_{T_{J}+\delta_{J}}$, possibly without Cramér’s condition, as that if we
assume Cramér’s condition and approximate $F_{T_{J}}$; where we use the same
Edgeworth expansion formula across both cases.. Now we are ready to present
our main result of this subsection. Let
$\xi_{k}^{2}:=\mathrm{Var}(g_{k}(X_{[1:k]}))$. Define the population Edgeworth
expansion formula for $T_{J}$ to be

$\displaystyle G_{{\cal
J}_{n,\alpha}}(u):=\Phi(u)+\phi(u)\cdot\Bigg{[}\Bigg{(}-\dfrac{\sum_{i=1}^{n}a_{n,\alpha;1}^{3}(i)(u^{2}-1)}{6\xi_{1}^{3}\\{\sum_{i^{\prime}=1}^{n}a_{n,\alpha;1}^{2}(i^{\prime})\\}^{3/2}}+\dfrac{\sum_{i=1}^{n}a_{n,\alpha;1}(i)u^{2}}{2\\{\sum_{i^{\prime}=1}^{n}a_{n,\alpha;1}^{2}(i^{\prime})\\}^{1/2}n\xi_{1}^{3}}\Bigg{)}\mathbb{E}[g_{1}^{3}(X_{1})]$
$\displaystyle+\Bigg{(}-\dfrac{\sum_{1\leq i<j\leq
n}a_{n,\alpha;1}(i)a_{n,\alpha;1}(j)a_{n,\alpha;2}(\\{i,j\\})(u^{2}-1)}{\\{\sum_{i^{\prime}=1}^{n}a_{n,\alpha;1}^{2}(i^{\prime})\\}^{3/2}\xi_{1}^{3}}$
$\displaystyle+\dfrac{(r-1)\sum_{i=1}^{n}a_{n,\alpha;1}(i)u^{2}}{\\{\sum_{i^{\prime}=1}^{n}a_{n,\alpha;1}^{2}(i^{\prime})\\}^{1/2}n\xi_{1}^{3}}\Bigg{)}\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]-\sum_{\ell=1}^{\lfloor\frac{\alpha/2}{\alpha-1}\rfloor}\Bigg{\\{}\frac{H_{2\ell-1}(u)}{(2\ell)!\\{\sum_{i^{\prime}=1}^{n}a_{n,\alpha;1}^{2}(i^{\prime})\\}^{\ell}\xi_{1}^{2\ell}}$
$\displaystyle\times\sum_{k_{1},\ldots,k_{\ell}\in[2:r]}\sum_{\begin{subarray}{c}\big{(}I_{k_{1}}^{(1)},\ldots,I_{k_{\ell}}^{(\ell)}\big{)}\in{\cal
C}_{n}^{k_{1}}\otimes\cdots\otimes{\cal C}_{n}^{k_{\ell}}\\\
\textrm{$I_{k_{1}}^{(1)},\ldots,I_{k_{\ell}}^{(\ell)}$ mutually
disjoint}\end{subarray}}a_{n,\alpha;k_{1}}^{2}(I_{k_{1}}^{(1)})\cdots
a_{n,\alpha;k_{\ell}}^{2}(I_{k_{\ell}}^{(\ell)})\xi_{k_{1}}^{2}\cdots\xi_{k_{\ell}}^{2}\Bigg{\\}}\Bigg{]},$
(13)

where
$H_{k}(u):=(-1)^{k}e^{u^{2}/2}{\mathrm{d}}^{k}/{\mathrm{d}}u^{k}(e^{-u^{2}/2})$
is the $k$th Hermite polynomial (Slater (1960), page 99). In practice, we use
the empirical version of (13) with estimated coefficient components. We stress
that under many ${\cal J}_{n,\alpha}$ designs, the last term in (13) can be
greatly simplified, thus we still comply with the $O(n^{\alpha})$ computation
(see Sections 2.2.1, 2.2.2). Define

$\displaystyle\widetilde{\mathbb{E}}[g_{1}^{3}(X_{1})]:=\frac{1}{n}\sum_{i=1}^{n}h(X_{[i:(i+r-1)]})$
$\displaystyle
h(X_{\\{i,[(i+r):(i+2r-2)]\\}})h(X_{\\{i,[(i+2r-1):(i+3r-3)]\\}})-\widetilde{\mu}^{3}$
(14)
$\displaystyle\widetilde{\mathbb{E}}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]:=$
$\displaystyle~{}\frac{1}{n}\sum_{i=1}^{n}h(X_{[(i-r+1):i]})h(X_{[i:(i+r-1)]})h(X_{[(i+r-1):(i+2r-2)]})$
$\displaystyle-\widetilde{\mu}^{3}-2U_{J}\cdot\widetilde{\xi}_{1}^{2}$ (15)
$\displaystyle\widetilde{\xi}_{k}^{2}:=\frac{1}{n^{\alpha}}\sum_{i=1}^{n}\sum_{d=1}^{n^{\alpha-1}}h(X_{[i:d:(i+(r-1)d)]})$
$\displaystyle
h(X_{[(i+(k-1)d):(-d):(i-(r-k)d)]})-\widetilde{\mu}^{2}-\sum_{k^{\prime}=1}^{k-1}\binom{k}{k^{\prime}}\widetilde{\xi}_{k^{\prime}}^{2},$
(16)

for $k\in[2:r]$. Let $\widetilde{G}_{{\cal J}_{n,\alpha}}(u)$ be the empirical
version of $G_{{\cal J}_{n,\alpha}}(u)$ with coefficients estimated by (5),
(14), (15) and (16). We will also write $G_{J}(u)$ in short for $G_{{\cal
J}_{n,\alpha}}(u)$, especially in the proof. Now we present our main theorem
on the distribution approximation error.

######  Theorem 2.1.

Set $\alpha\in(1,2]$. If $\xi_{1}>0$ and ${\cal J}_{n,\alpha}$ satisfies
Assumptions 1 and 2, then

$\displaystyle\big{\|}F_{T_{J}+\delta_{J}|{\cal J}_{n,\alpha}}(u)-G_{{\cal
J}_{n,\alpha}}(u)\big{\|}_{\infty}=$
$\displaystyle~{}O\big{(}n^{-\alpha/2}\log^{1/2}n\big{)}$ (17)
$\displaystyle\big{\|}F_{T_{J}+\delta_{J}|{\cal
J}_{n,\alpha}}(u)-\widetilde{G}_{{\cal J}_{n,\alpha}}(u)\big{\|}_{\infty}=$
$\displaystyle~{}{\widetilde{O}_{p}}(n^{-\alpha/2}\log^{1/2}n).$ (18)

If ${\cal J}_{n,\alpha}$ is deterministic, then $F_{T_{J}+\delta_{J}|{\cal
J}_{n,\alpha}}(u)$ should be understood as $F_{T_{J}+\delta_{J}}(u)$.

######  Remark 2.1.

Theorem 2.1 provides the practical guidance that for non-degenerate
U-statistics, paying more than $O(n^{2})$ computation (i.e. ,$\alpha>2$) will
not further merit risk control accuracy, since the error bound at $\alpha=2$
already matches the risk control accuracy of the complete U-statistic
(Helmers, 1991; Maesono, 1997). Also, increasing $\alpha$ beyond 2 only brings
$O(n^{-2})$ improvement to $\mathrm{Var}(U_{J})$ (Lee, 2019). Considering the
computational cost grows exponentially in $\alpha$, it is therefore not
worthwhile to set $\alpha$ above 2 in this scenario.

######  Remark 2.2.

Remark 3.1 in Chen and Kato (2019a) points out that as $\alpha$ decreases
towards 1, then $\sigma_{J;1}$ becomes a poorer approximation to $\sigma_{J}$;
when $\alpha=1$, their difference no longer diminishes with $n\to\infty$. This
is called _phase change_ in Weber (1981); Chen and Kato (2019a). While Weber
(1981); Chen and Kato (2019a) exclusively focused on $\mathrm{Var}(U_{J})$, we
reveal an interesting second aspect of what happens as $\alpha\to 1$: the
Edgeworth expansion becomes lengthier, and the risk control accuracy also
depreciates. If we do not incorporate an increasing number of bias-correction
terms in the Edgeworth expansion, it will depreciate the risk control accuracy
even faster (the $n^{-\alpha/2}$ term in Theorem 2.1 will be replaced by
$n^{-(\alpha-1)}$, which is the Berry-Esseen bound of the normal approximation
to $T_{J}$ that uses an estimated $\sigma_{J;1}$ to studentize).

Next, we invert the Edgeworth expansion to formulate the Cornish-Fisher
confidence interval with higher-order accurate confidence level control. This
inversion, however, is a quite nontrivial generalization of the convention
(Hall, 1983, 2013), because the Edgeworth expansion of $T_{J}$ contains not
only a $O(n^{-1/2})$ term, but also $O\big{(}n^{-(\alpha-1)\ell}\big{)}$,
$\ell\in\\{1,\ldots,\lfloor\alpha/(2(\alpha-1))\rfloor\\}$ terms. To ease
notation, we summarize (13) into the following prototype formula

$\displaystyle G_{{\cal J}_{n,\alpha}}(u):=$
$\displaystyle~{}\Phi(u)+\phi(u)\Bigg{\\{}\frac{\Gamma_{0}(u)}{\sqrt{n}}+\sum_{\ell=1}^{\lfloor\frac{\alpha/2}{\alpha-1}\rfloor}\frac{\Gamma_{\ell}(u)}{{M_{\alpha}}^{\ell}}\Bigg{\\}}$
(19)

Recall that we set ${M_{\alpha}}:=|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\asymp n^{\alpha-1}$ due
to Assumption 2. The other components $\Gamma_{0}(\cdot)$ and
$\Gamma_{\ell}(\cdot)$’s can be easily defined by comparing (13) to (19).
Given any user-selected significance level $\beta\in(0,1/2)$, define the
Cornish-Fisher expansion terms $\Psi_{0}(z_{\beta})$ and
$\Psi_{\ell}(z_{\beta})$ for
$\ell=1,\ldots,\lfloor(\alpha/2)/(\alpha-1)\rfloor$ recursively as follows.
First, set

$\displaystyle\Psi_{0}(z_{\beta}):=$
$\displaystyle~{}\Gamma_{0}(z_{\beta}),\quad\Psi_{1}(z_{\beta}):=-\Gamma_{1}(z_{\beta}).$
(20)

Then for each $k\in\big{[}2:\lfloor(\alpha/2)/(\alpha-1)\rfloor\big{]}$ and
with the already-defined $\Psi_{0}(z_{\beta}),\ldots,\Psi_{k-1}(z_{\beta})$,
we define $\Psi_{k}(z_{\beta})$ by solving the following recursive equation

$\displaystyle\Psi_{k}(z_{\beta})=-\sum_{\ell^{\prime}=2}^{k}\Bigg{\\{}\sum_{\begin{subarray}{c}j_{1},\ldots,j_{\ell^{\prime}}:\\\
1\leq\\{j_{1},\ldots,j_{\ell^{\prime}}\\}\leq k-\ell^{\prime}+1\\\
j_{1}+\cdots+j_{\ell^{\prime}}=k\end{subarray}}\Psi_{j_{1}}(z_{\beta})\cdots\Psi_{j_{\ell^{\prime}}}(z_{\beta})\cdot\frac{\phi^{(\ell^{\prime}-1)}(z_{\beta})}{(\ell^{\prime})!}\Bigg{\\}}$
$\displaystyle-\sum_{\begin{subarray}{c}k_{1},k_{2}:k_{1}+k_{2}=k\\\
k_{1}=0,\ldots,k-1\\\
\end{subarray}}\Bigg{[}\Bigg{\\{}\sum_{\ell^{\prime}=1}^{k_{1}}\sum_{\begin{subarray}{c}j_{1},\ldots,j_{\ell^{\prime}}:\\\
1\leq\\{j_{1},\ldots,j_{\ell^{\prime}}\\}\leq k-\ell^{\prime}+1\\\
j_{1}+\cdots+j_{\ell^{\prime}}=k\end{subarray}}\Psi_{j_{1}}(z_{\beta})\cdots\Psi_{j_{\ell^{\prime}}}(z_{\beta})\cdot\frac{\phi^{(\ell^{\prime}-1)}(z_{\beta})}{(\ell^{\prime})!}\Bigg{\\}}$
$\displaystyle\cdot\Bigg{\\{}\Gamma_{k_{2}}(z_{\beta})+\sum_{\ell^{\prime}=1}^{k_{2}-1}\sum_{\ell^{\prime\prime}=1}^{k-\ell^{\prime}}\Bigg{\\{}\sum_{\begin{subarray}{c}j_{1},\ldots,j_{\ell^{\prime\prime}}:\\\
1\leq\\{j_{1},\ldots,j_{\ell^{\prime\prime}}\\}\leq
k-\ell^{\prime}-\ell^{\prime\prime}+1\\\
j_{1}+\cdots+j_{\ell^{\prime\prime}}=k-\ell^{\prime}\end{subarray}}\Psi_{j_{1}}(z_{\beta})\cdots\Psi_{j_{\ell^{\prime\prime}}}(z_{\beta})\cdot\frac{\Gamma_{\ell^{\prime}}^{(\ell^{\prime\prime})}(z_{\beta})}{(\ell^{\prime\prime})!}\Bigg{\\}}\Bigg{\\}}\Bigg{]}.$
(21)

Our population Cornish-Fisher expansion that approximates the lower
$\beta$-quantile of $T_{J}$ is

$\displaystyle G_{{\cal J}_{n,\alpha}}^{-1}(z_{\beta})=:$
$\displaystyle~{}z_{\beta}-\frac{\Psi_{0}(z_{\beta})}{\sqrt{n}}+\sum_{\ell=1}^{\lfloor\frac{\alpha/2}{\alpha-1}\rfloor}\frac{\Psi_{\ell}(z_{\beta})}{{M_{\alpha}}^{\ell}}.$
(22)

Next, to define the empirical C-F expansion, we first re-express the empirical
Edgeworth expansion $\widetilde{G}_{{\cal J}_{n,\alpha}}$ in the format of
(19), which defines the empirical versions of $\Gamma_{0},\Gamma_{1},\ldots$,
denoted by $\widetilde{\Gamma}_{0},\widetilde{\Gamma}_{1},\ldots$. Replacing
$\Gamma$’s in (20) and (21) by $\widetilde{\Gamma}$’s, we obtain the empirical
$\widetilde{\Psi}$’s and thus the empirical C-F expansion
$\widetilde{G}_{{\cal J}_{n,\alpha}}^{-1}(\cdot)$. We omit further details due
to page limit.

######  Theorem 2.2.

Under the conditions of Theorem 2.1, for any given $\beta\in(0,1)$, the
Cornish-Fisher expansion (22) satisfies

$\displaystyle F_{T_{J}+\delta_{J}|{\cal J}_{n,\alpha}}\big{(}G_{{\cal
J}_{n,\alpha}}^{-1}(z_{\beta})\big{)}=$
$\displaystyle~{}\beta+O\big{(}n^{-\alpha/2}\big{)}$ (23) $\displaystyle
F_{T_{J}+\delta_{J}|{\cal J}_{n,\alpha}}\big{(}\widetilde{G}_{{\cal
J}_{n,\alpha}}^{-1}(z_{\beta})\big{)}=$
$\displaystyle~{}\beta+{\widetilde{O}_{p}}\big{(}n^{-\alpha/2}\log^{1/2}n\big{)}.$
(24)

Theorems 2.1 and 2.2 enable inference methods with higher-order accurate risk
control.

######  Corollary 2.1.

Under the conditions of Theorem 2.1, we have

  1. (a).

Test the hypotheses

$H_{0}:\mu=\mu_{0};\quad\textrm{vs.}\quad H_{a}:\mu\neq\mu_{0}\\\ $

using the following Empirical p-value, denoted by $\mathfrak{p}$

$\displaystyle\mathfrak{p}:=$
$\displaystyle~{}2\min\Big{\\{}\widetilde{G}_{{\cal J}_{n,\alpha}}(T_{J}^{(\rm
obs)}+\delta_{J}),1-\widetilde{G}_{{\cal J}_{n,\alpha}}(T_{J}^{(\rm
obs)}+\delta_{J})\Big{\\}},$ (25)

where $T_{J}^{(\rm obs)}:=(U_{J}-\mu_{0})/\big{\\{}|{\cal
J}_{n,\alpha}|^{-1}\\{\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\\}^{1/2}\widetilde{\xi}_{1}\big{\\}}$.
This test enjoys higher-order accurate type-I error control:

$\displaystyle\mathbb{P}_{H_{0}}\big{(}\mathfrak{p}<\beta\big{|}{\cal
J}_{n,\alpha}\big{)}=$ $\displaystyle~{}\beta+O(n^{-\alpha/2}\log^{1/2}n).$

  2. (b).

The Cornish-Fisher confidence interval ${\cal I}_{\beta}$ defined by

$\displaystyle{\cal I}_{\beta}:=\Big{(}$ $\displaystyle
U_{J}-(\widetilde{G}_{{\cal
J}_{n,\alpha}}^{-1}(z_{1-\beta/2})-\delta_{J})\cdot|{\cal
J}_{n,\alpha}|^{-1}\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{1/2}\cdot\widetilde{\xi}_{1},$
$\displaystyle U_{J}-(\widetilde{G}_{{\cal
J}_{n,\alpha}}^{-1}(z_{\beta/2})-\delta_{J})\cdot|{\cal
J}_{n,\alpha}|^{-1}\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{1/2}\cdot\widetilde{\xi}_{1}\Big{)}$

enjoys a higher-order accurate control of actual coverage probability around
$1-\beta$:

$\displaystyle\mathbb{P}\big{(}\mu\in{\cal I}_{\beta}\big{|}{\cal
J}_{n,\alpha}\big{)}=$ $\displaystyle~{}1-\beta+O(n^{-\alpha/2}\log^{1/2}n).$
(26)

In both parts, $\mathbb{P}(\cdot|{\cal J}_{n,\alpha})$ are taken with respect
to the randomness of both $X_{[1:n]}$ and $\delta_{J}$.

As mentioned in Section 1.3, reducing the U-statistic inflates
$\mathrm{Var}(U_{J})$. However, we studentize $U_{J}$ by
$\widetilde{\sigma}_{J;1}$, which only captures the leading term in
$\mathrm{Var}(U_{J})$, whose order does _not_ vary with $\alpha$. Readers
naturally wonder where the variance inflation is reflected in our statistical
inference procedure. Here, we use our CI formula as an example to clarify.

######  Remark 2.3.

The radius of our Cornish-Fisher CI is
$O\big{(}n^{-1/2}+n^{-(\alpha-1/2)}\big{)}$444 To see this, notice that
$\Gamma_{0}(-u)=\Gamma_{0}(u)$, while $\Gamma_{\ell}(-u)=-\Gamma_{\ell}(u)$
for all $\ell\geq 1$. Also notice that $(\widetilde{G}_{{\cal
J}_{n,\alpha}}^{-1}(z_{1-\beta/2})-\delta_{J})\asymp 1+n^{-(\alpha-1)}$ and
$|{\cal
J}_{n,\alpha}|^{-1}\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{1/2}\asymp
n^{-1/2}$. . Studentizing $U_{J}$ with $\widetilde{\sigma}_{J}$ will also
yield a CI radius of
$\\{O(n^{-1}+n^{-\alpha})\\}^{1/2}=O\big{(}n^{-1/2}+n^{-(\alpha-1/2)}\big{)}$.
In other words, using $\widetilde{\sigma}_{J}$ or $\widetilde{\sigma}_{J;1}$
to studentize $U_{J}$ lead to different pivots as intermediate steps, but
eventually, they will produce essentially the same Cornish-Fisher CI’s.

We conclude this subsection with the observation that any test, based on an
asymptotically normal $T_{J}$, that correctly estimates the leading term of
$\mathrm{Var}(U_{J})$, namely, $\sigma_{J}^{2}$, is asymptotically power-
optimal (see how Theorem 3.5 of Banerjee and Ma (2017) establishes the
asymptotic power-optimality of their test). This of course also applies to our
method. We reiterate that power-optimality is conceptually distinct from risk
control accuracy (see Section 1.3).

###  2.2 Applications: statistical inference with higher-order accurate risk
control for deterministic and randomized designs

In this subsection, we apply our general framework in Section 2.1 to analyzing
several popular designs. In Section 2.2.1, we propose and analyze a novel
variance-optimal deterministic reduction scheme. In Section 2.2.2, we present
the first provably higher-order accurate inference for randomized incomplete
U-statistics (Lee, 2019; Chen and Kato, 2019a).

####  2.2.1 A novel variance-optimal deterministic design

As aforementioned, most past works on U-statistic reduction focus on designing
a deterministic ${\cal J}_{n,\alpha}$ that minimizes $\mathrm{Var}(U_{J})$.
This task demands a design ${\cal J}_{n,\alpha}$, under which all
$a_{n,\alpha;k}(I_{k})$’s take as similar values as possible; specifically,
for $k>\alpha$, we desire $a_{n,\alpha;k}(I_{k})\in\\{0,1\\}$. Despite the
significant efforts in existing literature, it remains a challenging open
problem to construct such ${\cal J}_{n,\alpha}$ for a general $\alpha>1$.
However, our analysis reveals that increasing $\alpha$ above 2 brings little
merit to both inference power and risk control accuracy. Therefore, we can
naturally consider the much more doable problem with $\alpha\in(1,2]$.
Slightly abusing notation, set ${M_{\alpha}}:=\lfloor
n^{\alpha-1}\rfloor$555This ${M_{\alpha}}$ will later play the same role as
the ${M_{\alpha}}$ defined in (22), so we do not set up two symbols.. Define
an $r$-tuple ${\cal A}_{i,d}^{n,r}$ as follows.

$\displaystyle{\cal A}_{i,d}^{n,r}:=$
$\displaystyle\Big{\\{}(i_{1},\ldots,i_{r}):1\leq i_{k}\leq n,i_{k}\equiv
j_{k}\ (\mathrm{mod}\ n),j_{k}=i+(2^{k-1}-1)d,\forall k\in[1:r]\Big{\\}}.$
(27)

Our proposed deterministic design ${\cal J}_{n,\alpha}$ is

$\displaystyle{\cal J}_{n,\alpha}:=$
$\displaystyle\bigcup_{d=b_{1}{M_{\alpha}}}^{b_{2}{M_{\alpha}}}\bigcup_{i=1}^{n}{\cal
A}_{i,d}^{n,r},$ (28)

where $b_{1},b_{2}:0<b_{1}<b_{2}$ are constants. It is not difficult to check
that $|{\cal J}_{n,\alpha}|=O(n^{\alpha})$.

######  Lemma 2.2.

Suppose $n\gg r$. Set $\alpha\in(1,2)$ and
$b_{1}/b_{2}\in\big{(}(2^{r-1}-1)/(2^{r-1}),1\big{)}$. Our design (27) and
(28) satisfies $a_{n,\alpha;1}(i)\asymp{M_{\alpha}}$ for all $i\in[1:n]$ and
$a_{n,\alpha;k}(I_{k})\in\\{0,1\\}$, for all $k\geq 2$ and all $I_{k}\in{\cal
C}_{n}^{k}$. Moreover, $\sum_{I_{k}\in{\cal
C}_{n}^{k}}a_{n,\alpha;k}^{2}(I_{k})\asymp n^{\alpha+(\alpha-k)\vee 0}$ for
all $k\in[1:r]$.

Lemma 2.2 implies that our design achieves the optimal variance among all
data-oblivious reduction schemes under $O(n^{\alpha})$ computational budget
limit. As an immediate corollary of Lemma 2.2, our deterministic ${\cal
J}_{n,\alpha}$ also satisfies Assumption 2, so Theorem 2.1 and Corollary 2.1
hold. It has a handy Edgeworth expansion simplification

$\displaystyle\Gamma_{0}(u):=$
$\displaystyle~{}\frac{2u^{2}+1}{6\xi_{1}^{3}}\mathbb{E}[g_{1}^{3}(X_{1})]+\frac{(r-1)(u^{2}+1)}{2\xi_{1}^{3}}\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]$
(29) $\displaystyle\Gamma_{\ell}(u):=$
$\displaystyle~{}-\Bigg{\\{}\frac{\sum_{k=2}^{r}\xi_{k}^{2}\binom{r}{k}}{(b_{2}-b_{1})r^{2}\xi_{1}^{2}}\Bigg{\\}}^{\ell}\cdot\frac{H_{2\ell-1}(u)}{(2\ell)!}=-\Bigg{\\{}\frac{\sigma_{h}^{2}-r\xi_{1}^{2}}{(b_{2}-b_{1})r^{2}\xi_{1}^{2}}\Bigg{\\}}^{\ell}\cdot\frac{H_{2\ell-1}(u)}{(2\ell)!}$
(30)

for all $\ell=1,\ldots,\lfloor\alpha/\\{2(\alpha-1)\\}\rfloor$, where
$\sigma_{h}^{2}:=\mathrm{Var}(h(X_{[1:r]}))$. We can estimate $\sigma_{h}^{2}$
by

$\displaystyle\widetilde{\sigma}_{h}^{2}:=$
$\displaystyle~{}\frac{1}{n{M_{\alpha}}}\sum_{i=1}^{n}\sum_{d=1}^{M_{\alpha}}h^{2}(X_{[i:d:(i+(r-1)d)]})-\widetilde{\mu}^{2},$
(31)

where $\widetilde{\mu}$ is inherited from (5). Plugging (29), (30) and our
setting ${M_{\alpha}}=\lfloor n^{\alpha-1}\rfloor$ into (19), we obtain the
population Edgeworth expansion $G_{J}(u)$ for $T_{J}+\delta_{J}$ under our
${\cal J}_{n,\alpha}$. Using the plug-in estimators (14), (15) and (31), we
obtain the empirical Edgeworth expansion $\widetilde{G}_{J}(u)$. Finally,
using (20), (21), and (22), we can formulate the Cornish-Fisher CI.
Theoretical guarantees for our proposed design are straight corollaries of our
main theorems in Section 2.1 and thus omitted here due to page limit.

Interestingly, our method, while being an acceleration tool itself, can also
service some other acceleration tools and make them even faster. One such
“customer” is the _divide-and-conquer acceleration_ by parallel computing in
Chen and Peng (2021). We have $K$ parallel computing servers who will return
summary statistics to the main server for aggregation. In Chen and Peng
(2021), each individual server still computes a _complete_ U-statistic, thus
readers already see that our method can significantly improve on this point.
Due to page limit, we sink the details of coupling our method with Chen and
Peng (2021) to Algorithm 1 in Supplementary Material. We stress that Algorithm
1 can also be viewed as a parallel-computing strategy for implementing our
inference method on several servers, while the inference formula and
accompanying theory remain unchanged. We compare our method coupled with Chen
and Peng (2021) versus the vanilla Chen and Peng (2021) in Table 1, where for
cleanness, we equate all split sizes, set $K\asymp n^{\tau^{\prime}}$ as in
Chen and Peng (2021) and equate the orders of the second leading terms in the
two approaches’ variances by setting $\alpha=2-\tau^{\prime}$. Table 1 shows
that our method significantly accelerates Chen and Peng (2021), with little to
no variance inflation (depending on $\tau^{\prime}$ vs $\alpha$) and higher
accuracy in risk control.

Table 1: Merits of our method coupled with Chen and Peng (2021). Set $\alpha\in(1,2)$. | Vanilla Chen and Peng (2021) | Our method + Chen and Peng (2021)
---|---|---
Time cost on each server | $O(n^{r}/K^{r-1})$ | $O\big{(}(n/K)^{\alpha-1}\cdot K\big{)}$
Variance of aggregated U-stat. | $r^{2}\xi_{1}^{2}/n+O(Kn^{-2})$ | $r^{2}\xi_{1}^{2}/n+O(n^{-\alpha})$
CDF approximation error |  $o(n^{-1/2})$666This requires $K=O(n^{\tau^{\prime}})$ for $\tau^{\prime}\in(0,1/4)$. See Theorem 3.3-(i) in Chen and Peng (2021). | $O(n^{-\alpha/2})$
Inference error777Inference error may be represented by the discrepancy between the nominal and the actual CI coverage probabilities or that between the nominal significance level and the actual type I error rate. |  $o_{p}(1)$ 888 The vanilla Chen and Peng (2021) _standardizes_ $U_{J}$. However, standardization-based inference needs a bias-correction that de-facto reproduces the result of the studentization-based inference in order to achieve higher-order accuracy. See Hall (2013), Section 3.10.2. | ${\widetilde{O}_{p}}(n^{-\alpha/2}\log^{1/2}n)$

####  2.2.2 Analysis of randomized incomplete U-statistics

The general framework we propose in Section 2.1 is particularly convenient for
analyzing randomized designs. Here, we showcase the analysis for a few example
designs in literature:

  1. (J1)

Sample $n^{\alpha}$ size-$r$ subsets from ${\cal C}_{n}^{r}$ at random, with
replacement.

  2. (J2)

Similar to (J1) but sample without replacement999Notice that sampling ${\cal
J}_{n,\alpha}:|{\cal J}_{n,\alpha}|=O(n^{\alpha})$ without replacement could
also be done with $O(n^{\alpha})$ budget, in both time and memory, by a
lexicographic indexing of elements in ${\cal C}_{n}^{r}$..

  3. (J3)

For $i=1,\ldots,n$, sample $n^{\alpha-1}$ size-$r$ subsets from ${\cal
C}_{n}^{r}$ containing $i$, with replacement.

  4. (J4)

Similar to (J3), but for each $i$, sample without replacement.

These sampling schemes are natural to come up with, and there are many more
similar randomized designs in existing literature. However, no available
theory and methods yet exist to address the problem of higher-order accurate
risk control for any of them. The conventional analysis typically starts with
re-expressing $U_{J}$ as follows (Chen and Kato, 2019a).

$\displaystyle U_{J}-\mu:=$
$\displaystyle~{}\underbrace{(U_{n}-\mu)}_{\textrm{(Part
I)}}+\underbrace{|{\cal J}_{n,\alpha}|^{-1}\sum_{I_{r}\in{\cal
J}_{n,\alpha}}\big{\\{}h(X_{I_{r}})-U_{n}\big{\\}}}_{\textrm{(Part
II)}}=:(U_{n}-\mu)+V_{J},$ (32)

where part I is a rescaled complete U-statistic (see definition in eq. (1))
and part II captures the randomness in ${\cal J}_{n,\alpha}$. Then one can
normal-approximate both parts, respectively, to establish the asymptotic
normality of $U_{J}$ via careful conditioning and convolution, see pages 9–20
in Chen and Kato (2019b). While (32) is very useful for analyzing degenerate
U-statistics (see our Section 3), it is not a sharp tool in the non-degenerate
case, where both parts will contribute non-ignorable terms to the Edgeworth
expansion of studentized $U_{J}$. Therefore, our analysis here takes a very
different route: the key is to apply our general framework in Section 2.1. To
this end, we first verify that our considered randomized design ${\cal
J}_{n,\alpha}$ satisfies Assumption 2 with high probability. (Assumption 1 is
easily verified.)

######  Lemma 2.3.

Let ${\cal J}_{n,\alpha}$ be constructed by one of (J1)–(J4). For any given
constant $C_{1}>0$, there exist constants $C_{2},C_{3}:C_{3}>C_{2}>0$
depending on $C_{1}$ and the design ${\cal J}_{n,\alpha}$, such that
Assumption 2 holds with probability at least $1-n^{-C_{1}}$.

Next, we showcase the formula simplification for (J1) and (J3) by taking
another layer of expectation with respect to the randomness of ${\cal
J}_{n,\alpha}$, denoted by $\mathbb{E}_{J}$, as follows.

######  Corollary 2.2.

Set $\alpha\in(1,2)$. For randomized deign (J1), we have

$\displaystyle\mathbb{E}_{J}[\Gamma_{0}(u)]:=$
$\displaystyle~{}\frac{2u^{2}+1}{6\xi_{1}^{3}}\mathbb{E}[g_{1}^{3}(X_{1})]+\frac{(r-1)(u^{2}+1)}{2\xi_{1}^{3}}\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]$
(33) $\displaystyle\mathbb{E}_{J}[\Gamma_{\ell}(u)]:=$
$\displaystyle~{}-\frac{H_{2\ell-1}(u)}{(2\ell)!}\Big{\\{}\frac{\sum_{k=2}^{r}\binom{r}{k}\xi_{k}^{2}}{r^{2}\xi_{1}^{2}}\Big{\\}}^{\ell}.$
(34)

For randomized deign (J3), we have

$\displaystyle\mathbb{E}_{J}[$
$\displaystyle\Gamma_{0}(u)]:=\frac{2u^{2}+1}{6\xi_{1}^{3}}\mathbb{E}[g_{1}^{3}(X_{1})]$
$\displaystyle+\frac{(r-1)\big{\\{}(r^{3}+2r^{2}-2)u^{2}+r^{3}-2r^{2}+2\big{\\}}}{2r^{3}\xi_{1}^{3}}\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]$
(35) $\displaystyle\mathbb{E}_{J}[$ $\displaystyle\Gamma_{\ell}(u)]\textrm{ is
the same as the $\mathbb{E}_{J}[\Gamma_{\ell}(u)]$ under \ref{jna-
choice-1:SRS-w-replace}}.$ (36)

In both examples, set

$\displaystyle G_{J}(u):=$
$\displaystyle~{}\Phi(u)+\phi(u)\Big{\\{}\frac{\mathbb{E}_{J}[\Gamma_{0}(u)]}{\sqrt{n}}+\sum_{\ell=1}^{\lfloor\frac{\alpha/2}{\alpha-1}\rfloor}\frac{\mathbb{E}_{J}[\Gamma_{\ell}(u)]}{{M_{\alpha}}^{\ell}}\Big{\\}},$
(37)

where

$\displaystyle{M_{\alpha}}:=$
$\displaystyle~{}n^{\alpha-1}\cdot\big{\\{}1+1/(rn^{\alpha-1})\big{\\}}\Big{/}\begin{cases}1+\dfrac{n^{\alpha-2}\xi_{2}^{2}\cdot
r(r-1)}{\sum_{k=2}^{r}\binom{r}{k}\xi_{k}^{2}},&\textrm{under \ref{jna-
choice-1:SRS-w-replace}}\\\ 1+\dfrac{n^{\alpha-2}\xi_{2}^{2}\cdot
r^{2}(r-1)}{2\sum_{k=2}^{r}\binom{r}{k}\xi_{k}^{2}},&\textrm{under \ref{jna-
choice-3:stratified-w-replace}}\end{cases}.$ (38)

Then we have

$\displaystyle\big{\|}F_{T_{J}+\delta_{J}}(u)-G_{J}(u)\big{\|}_{\infty}=$
$\displaystyle~{}O(n^{-\alpha/2}\log n).$ (39)

We can naturally define the empirical version $\widetilde{G}_{J}(u)$ with
coefficient estimated by (5), (14), (15) and (31) and use it for downstream
analysis, accompanied by theoretical guarantees exactly similar to Corollary
2.1. We skip the repetitive details.

##  3 Reduction of degenerate noiseless U-statistics

We call a U-statistic _degenerate up to degree- $(k_{0}-1)$_ (Chen and Kato,
2019a), if for some $k_{0}\in[2:r]$, we have $\xi_{1}=\cdots=\xi_{k_{0}-1}=0$
and $\xi_{k_{0}}^{2}>0$. Unlike the non-degenerate scenario, the incomplete
U-statistics with deterministic and randomized designs are mainly driven by
very different sources of stochastic variations, so they need to be treated
separately.

We first analyze $U_{J}$ with the deterministic design ${\cal J}_{n,\alpha}$
in Section 2.2.1, starting with variance estimation. Since this ${\cal
J}_{n,\alpha}$ has the convenient property that $\sum_{I_{k}\in{\cal
C}_{n}^{k}}a_{n,\alpha;k}^{2}(I_{k})=\binom{r}{k}|{\cal J}_{n,\alpha}|$ for
all $k\in[1:r]$, we have $\mathrm{Var}(U_{J})=|{\cal
J}_{n,\alpha}|^{-1}\sum_{k=k_{0}}^{r}\binom{r}{k}\xi_{k}^{2}=|{\cal
J}_{n,\alpha}|^{-1}\sigma_{h}^{2}$, where
$\sigma_{h}^{2}:=\mathrm{Var}(h(X_{[1:r]}))=\sum_{k=k_{0}}^{r}\binom{r}{k}\xi_{k}^{2}$
can be accurately estimated in $O(n^{\alpha})$ time by
$\widetilde{\sigma}_{h}^{2}:=|{\cal J}_{n,\alpha}|^{-1}\sum_{I_{r}\in{\cal
J}_{n,\alpha}}\Big{\\{}h\big{(}X_{I_{r}}\big{)}-U_{J}\Big{\\}}^{2}$. We
studentize $U_{J}$ as $T_{J}:=(U_{J}-\mu)/\big{(}|{\cal
J}_{n,\alpha}|^{-1/2}\cdot\widetilde{\sigma}_{h}\big{)}$.

######  Theorem 3.1.

Suppose $h$ is degenerate up to degree $(k_{0}-1)$ for some $k_{0}\in[2:r]$.
Set $\alpha\in(1,2)$. If ${\cal J}_{n,\alpha}$ is constructed as that in
Section 2.2.1, then as $n\to\infty$, $T_{J}\stackrel{{\scriptstyle
d}}{{\to}}N(0,1)$.

The downstream inference using Theorem 3.1 is completely standard. We skip it.

Next, we analyze $U_{J}$ with a randomized ${\cal J}_{n,\alpha}$. For
simplicity, we focus on two particular designs (J1) and (J3), where ${\cal
J}_{n,\alpha}$ elements are independent. The decomposition (32) from Chen and
Kato (2019a) now becomes particularly useful. In the degenerate scenario,
$V_{J}$ is the dominating term in (32) whenever $\alpha<k_{0}$. By Major
(2007), we know
$U_{n}-\mu={\widetilde{O}_{p}}\big{(}n^{-k_{0}/2}\log^{k_{0}/2}n\big{)}$.
Rewrite $V_{J}$ as $V_{J}=|{\cal J}_{n,\alpha}|^{-1}\sum_{j=1}^{|{\cal
J}_{n,\alpha}|}{\cal X}_{j}$, where ${\cal X}_{j}:=h(X_{I_{r}^{(j)}})-U_{n}$
by recalling (2). Conditioning on $X_{[1:n]}$, we can view $V_{J}$ as a sample
mean of independent mean-zero random variables. We have $|{\cal
J}_{n,\alpha}|\cdot\mathrm{Var}(V_{J}|X_{[1:n]})=\binom{n}{r}^{-1}\sum_{I_{r}\in{\cal
C}_{n}^{r}}\big{\\{}h(X_{I_{r}})-U_{n}\big{\\}}^{2}=\sigma_{h}^{2}+{\widetilde{O}_{p}}\big{(}n^{-k_{0}/2}\log^{k_{0}/2}n\big{)}$.
Estimating $\sigma_{h}^{2}$ by $\widetilde{\sigma}_{h}^{2}$, the
studentization is

$\displaystyle T_{J}:=(U_{J}-\mu)/\big{(}|{\cal
J}_{n,\alpha}|^{-1/2}\cdot\widetilde{\sigma}_{h}\big{)}={\widetilde{O}_{p}}\big{(}n^{(\alpha-
k_{0})/2}\log^{k_{0}/2}n\big{)}+V_{J}/(|{\cal
J}_{n,\alpha}|^{-1/2}\cdot\widetilde{\sigma}_{h}).$ (40)

To characterize the distribution of $T_{J}$, we study the second term on the
RHS of (40). We formulate an Edgeworth expansion that refines the normal
approximation (Chen and Kato, 2019a). Define

$\displaystyle G_{V|X_{[1:n]}}(u):=$
$\displaystyle~{}\Phi(u)+\frac{\mathbb{E}[{\cal
X}_{1}^{3}|X_{[1:n]}](2u^{2}+1)\phi(u)}{\big{\\{}\mathrm{Var}({\cal
X}_{1}|X_{[1:n]})\\}^{3/2}\cdot 6|{\cal J}_{n,\alpha}|^{1/2}}$ (41)
$\displaystyle\widetilde{G}_{V|X_{[1:n]}}(u):=$
$\displaystyle~{}\Phi(u)+\frac{\widetilde{\nu}_{h}^{3}(2u^{2}+1)\phi(u)}{\widetilde{\sigma}_{h}^{3}\cdot
6|{\cal J}_{n,\alpha}|^{1/2}},$ (42)

and define $\widetilde{\nu}_{h}^{3}:=|{\cal
J}_{n,\alpha}|^{-1}\sum_{I_{r}\in{\cal
J}_{n,\alpha}}\big{\\{}h(X_{I_{r}})-U_{J}\big{\\}}^{3}$ as the empirical
version of $\mathbb{E}[{\cal X}_{1}^{3}|X_{[1:n]}]$. The next theorem bounds
the distribution approximation errors of $G_{V|X_{[1:n]}}(u)$ and
$\widetilde{G}_{V|X_{[1:n]}}(u)$.

######  Theorem 3.2.

Suppose the U-statistic is degenerate up to degree $(k_{0}-1)$ for some
$k_{0}\in[2:r]$. Set $\alpha\in(1,k_{0})$. If ${\cal J}_{n,\alpha}$ is follows
(J1) or (J3) in Section 2.2.2, then

$\displaystyle\Big{\|}F_{T_{J}+\delta_{V}|X_{[1:n]}}(u)-G_{V|X_{[1:n]}}(u)\Big{\|}_{\infty}=$
$\displaystyle~{}{\widetilde{O}_{p}}\Big{(}n^{-\alpha}\log^{1/2}n+n^{(\alpha-
k_{0})/2}\log^{k_{0}/2}n\Big{)}$ (43)
$\displaystyle\Big{\|}F_{T_{J}+\delta_{V}|X_{[1:n]}}(u)-\widetilde{G}_{V|X_{[1:n]}}(u)\Big{\|}_{\infty}=$
$\displaystyle~{}{\widetilde{O}_{p}}\Big{(}n^{-\alpha}\log^{1/2}n+n^{(\alpha-
k_{0})/2}\log^{k_{0}/2}n\Big{)},$ (44)

where we add to $T_{J}$ an artificial Gaussian smoother independent of all
data $\delta_{V}\sim N(0,\sigma^{2}=Cn^{-\alpha}\log n)$ with a large enough
constant $C$.

To construct the Cornish-Fisher CI, we invert
$\widetilde{G}_{V|X_{[1:n]}}(u)$.

######  Corollary 3.1.

Define

$\displaystyle\widetilde{G}_{V|X_{[1:n]}}^{-1}(z_{\beta}):=$
$\displaystyle~{}z_{\beta}-\frac{\widetilde{\nu}_{h}^{3}(2z_{\beta}^{2}+1)}{\widetilde{\sigma}_{h}^{3}\cdot
6|{\cal J}_{n,\alpha}|^{1/2}}.$ (45)

and the two-sided Cornish-Fisher confidence interval for estimating $\mu$ to
be

$\displaystyle{\cal I}_{V}:=\Big{(}$ $\displaystyle
U_{J}-\big{\\{}\widetilde{G}_{V|X_{[1:n]}}^{-1}(1-z_{\beta}/2)-\delta_{V}\big{\\}}|{\cal
J}_{n,\alpha}|^{-1/2}\cdot\widetilde{\sigma}_{h},$ $\displaystyle
U_{J}-\big{\\{}\widetilde{G}_{V|X_{[1:n]}}^{-1}(z_{\beta}/2)-\delta_{V}\big{\\}}|{\cal
J}_{n,\alpha}|^{-1/2}\cdot\widetilde{\sigma}_{h}\Big{)},$ (46)

Then we have

$\displaystyle\mathbb{P}(\mu\in{\cal I}_{V})=$
$\displaystyle~{}1-\beta+O\big{(}n^{-\alpha}\log^{1/2}n+n^{(\alpha-
k_{0})/2}\log^{k_{0}/2}n\big{)}.$ (47)

Theorem 3.2 and Corollary 3.1 provide the first finite-sample characterization
of the trade-off between risk control accuracy and computational speed in the
reduction of degenerate U-statistics. We remark three quick interpretations of
these results. First, when $\alpha<k_{0}/2$, our method achieves a sharper
error bound than normal approximation. Second, our method’s error rate is
optimized to $O\big{(}n^{-k_{0}/3}\log^{k_{0}/2}n\big{)}$ at $\alpha=k_{0}/3$.
Third, for all choices of $\alpha<r$, our method’s error bound improves over
the $O\big{(}n^{(\alpha-k_{0})/4}+n^{-(\alpha\wedge 1)/6}\big{)}$ Berry-Esseen
bound (Chen and Kato, 2019a) for one-dimensional U-statistics.101010The first
term in the error bound of Theorem 3.3 in Chen and Kato (2019a) needs a slight
revision, changing “$\log^{k+3}d$” into “$\log^{k+3}(d+C_{0})$” for some
constant $C_{0}>0$. To see this, notice that when (the “$d$” in their paper)
$d=1$, Equation (C.19) in the proof of Theorem 3.3 (Chen and Kato, 2019b)
needs this modification to hold.

We conclude this subsection by an interesting comparison between the guidance
on selecting $\alpha$ in the non-degenerate and degenerate scenarios. For non-
degenerate U-statistics, the reason not to choose a large $\alpha$ is purely
computation (recall Remark 2.1); if computation resources permit, we can set
$\alpha=r$ to minimize $\mathrm{Var}(U_{J})$. The degenerate scenario,
however, is very different. As $\alpha$ approaches $k_{0}$, the risk control
accuracy of both our method and normal approximation will drop, because the
_non-Gaussian_ part $U_{n}-\mu$ becomes increasingly non-ignorable in (32).
This blessing reinstatement of normality by incompleteness was discovered by
Weber (1981); Chen and Kato (2019a), and our paper provides a much-refined
finite-sample description of this phenomenon. Approximating the distribution
of the non-Gaussian $U_{n}-\mu$ is quite challenging, and no finite-sample
result yet exists to our best knowledge (Van der Vaart, 2000; Giraudo, 2021).

##  4 Reduction of network moments as noisy U-statistics

Sections 2 and 3 comprehensively addressed noiseless U-statistics. Now, we
turn our attention to an important class of noisy U-statistics, called
_network moments_ , for which the data-oblivious reduction schemes is
particularly useful. Consider an undirected network with binary edges, where
each node is associated with a latent position $X_{i}\sim$ Uniform$[0,1]$, and
the edge probability between $(i,j)$ is $W_{i,j}:=\rho_{n}\cdot
f(X_{i},X_{j})$ (Aldous, 1981; Hoover, 1979), where $\rho_{n}$ is a global
sparsity parameter and $f(\cdot,\cdot)$ is a latent function, called _graphon
function_ , that encodes all network structural information. Then observed
data are $A_{i,j}|W\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}$
Bernoulli$(W_{i,j})$. A network moment, indexed by a motif $R$ of $r$ nodes
and $s$ edges, is the count of $R$ in the network. Specifically, let
$h(A_{[1:r]}):=\sum_{\pi\in{\rm Permu.([1:r])}}\mathbbm{1}_{A_{[\pi(1:r)]}\geq
R}$, in which we treat $R$ as its own adjacency matrix, and $\geq$ is element-
wise. Define the empirical network moment $\widehat{U}_{n}$ by

$\displaystyle\widehat{U}_{n}:=$
$\displaystyle~{}\binom{n}{r}^{-1}\sum_{I_{r}\in{\cal
C}_{n}^{r}}h(A_{I_{r}}).$ (48)

Denote $\mu:=\mathbb{E}[U_{n}]$. The goal is to perform statistical inference
for $\mu$.

The higher-order approximation to the distribution of the studentization of
the complete network moment $\widehat{U}_{n}$ has been established in Zhang
and Xia (2022). However, similar to the noiseless scenario, the expensive
$O(n^{r})$ computation complexity is also a major obstacle in statistical
inference. This leads us to consider U-statistic reduction. Notice that in
(48), both $f$ and $X_{i}$’s are not only unobserved, but even inestimable,
due to an identifiability issue Gao et al. (2015). Therefore, data-aware
reduction techniques such as (Kong and Zheng, 2021) are inapplicable, leaving
data-oblivious reduction the natural choice. The incomplete network
U-statistic with design ${\cal J}_{n,\alpha}$ is

$\displaystyle\widehat{U}_{J}:=\big{|}{\cal
J}_{n,\alpha}\big{|}^{-1}\sum_{I_{r}\in{\cal
J}_{n,\alpha}}h(A_{I_{r}})=\big{|}{\cal
J}_{n,\alpha}\big{|}^{-1}\sum_{I_{r}\in{\cal
C}_{n}^{r}}a_{n,\alpha;r}(I_{r})h(A_{I_{r}}).$ (49)

Like Zhang and Xia (2022), we decompose the stochastic variations in
$\widehat{U}_{J}$ into two parts:
$\widehat{U}_{J}-\mu=(U_{J}-\mu)+(\widehat{U}_{J}-U_{J})$, where we define
$U_{J}:=\mathbb{E}[\widehat{U}_{J}|X_{[1:n]}]=\mu+|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}(i)g_{1}(X_{i})+|{\cal
J}_{n,\alpha}|^{-1}\sum_{k=2}^{r}\sum_{I_{k}\in{\cal
C}_{n}^{k}}a_{n,\alpha;k}(I_{k})g_{k}(X_{I_{k}})$. We call a network moment
$\widehat{U}_{J}$ non-degenerate if the noiseless part in its decomposition
$U_{J}$ is non-degenerate. In this paper, we focus on non-degenerate network
moments, leaving the degenerate case to an interesting future work.

To design a variance estimator for studentizing $\widehat{U}_{J}$, notice that
$\mathrm{Var}(U_{J}-\mu|{\cal J}_{n,\alpha})\asymp\rho_{n}^{2s}n^{-1}$ while
$\mathrm{Var}(\widehat{U}_{J}-U_{J}|{\cal
J}_{n,\alpha})\asymp\rho_{n}^{2s-1}n^{-2}+\rho_{n}^{s}n^{-\alpha}$ (see proof
of Theorem 4.1). When $\rho_{n}^{s}n^{\alpha-1}\to\infty$, the randomness in
$U_{J}$ dominates that in $A|W$111111Notice that this condition is always
satisfied for _complete_ network U-statistics, see Zhang and Xia (2022): where
$\alpha=r$ (i.e. $U_{J}=U_{n}$), with the presence of Condition (ii) in their
Lemma 3.1.. We now assume $\rho_{n}^{s}n^{\alpha-1}\to\infty$ and relegate the
case $\rho_{n}^{s}n^{\alpha-1}\to 0$ to Theorem 4.2121212For simplicity, we
only study the case where $\rho_{n}^{s}n^{\alpha-1}\to 0$, leaving the
intermediate case to future work.. With this assumption, we can approach
$\mathrm{Var}(\widehat{U}_{J})$ via $\xi_{1}^{2}$. Define
$\widehat{a}_{i}:=\sum_{I_{r}\in{\cal J}_{n,\alpha},i\in
I_{r}}h(A_{I_{r}})/a_{n,\alpha;1}(i)$, and estimate $\xi_{1}^{2}$ by

$\displaystyle\widehat{\xi}_{1}^{2}:=$ $\displaystyle~{}r^{-1}|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}(i)\big{(}\widehat{a}_{i}-\widehat{U}_{J}\big{)}^{2}.$
(50)

There are several good alternative choices for $\widehat{\xi}_{1}^{2}$, e.g.
$\big{\\{}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\big{\\}}^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)(\widehat{a}_{i}-\widehat{U}_{J})^{2}$
and $n^{-1}\sum_{i=1}^{n}(\widehat{a}_{i}-\widehat{U}_{J})^{2}$; we choose
(50) for technical convenience. Studentize $\widehat{U}_{J}$ as

$\displaystyle\widehat{T}_{J}:=\frac{\widehat{U}_{J}-\mu}{\\{\sum_{i=1}^{n}a_{n,\alpha;1}^{2}(i)\\}^{1/2}\cdot\widehat{\xi}_{1}/|{\cal
J}_{n,\alpha}|}+\delta_{J},$ (51)

where, like before, we add an artificial noise $\delta_{J}$ to completely
waive Cramér’s condition. For a motif $R$, define the population Edgeworth
expansion for $\widehat{T}_{J}$ as follows.

$\displaystyle G_{J}(x):=$
$\displaystyle~{}\Phi(x)+\frac{\varphi(x)}{\sqrt{n}\cdot\xi_{1}^{3}}\Big{\\{}\frac{2x^{2}+1}{6}\cdot\mathbb{E}[g_{1}^{3}(X_{1})]$
$\displaystyle+\frac{r-1}{2}\cdot(x^{2}+1)\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]\Big{\\}}.$
(52)

To formulate the empirical Edgeworth expansion (EEE), we also need to estimate
$\mathbb{E}[g_{1}^{3}(X_{1})]$ and
$\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]$. Limited by the
$O(n^{\alpha})$ computational budget, we need to estimate these quantities
very differently than that in Zhang and Xia (2022). Define

$\displaystyle\widehat{\mathbb{E}}[g_{1}^{3}(X_{1})]:=$
$\displaystyle~{}r^{-1}|{\cal
J}_{n,\alpha}|^{-1}\sum_{i=1}^{n}a_{n,\alpha;1}(i)\big{(}\widehat{a}_{i}-\widehat{U}_{J}\big{)}^{3},$
(53) $\displaystyle\widehat{\mathbb{E}}[g_{1}(X_{1})g_{1}(X_{2})$
$\displaystyle g_{2}(X_{1},X_{2})]$ $\displaystyle:=$ $\displaystyle~{}|{\cal
J}_{n,\alpha}|^{-1}\sum_{I_{r}\in{\cal
J}_{n,\alpha}}\big{\\{}(h(A_{I_{r}})-\widehat{U}_{J})\cdot\big{(}\widehat{a}_{I_{r}(1)}-\widehat{U}_{J}\big{)}\cdot\big{(}\widehat{a}_{I_{r}(2)}-\widehat{U}_{J}\big{)}\big{\\}},$
(54)

where $I_{r}(\ell)$ indicates the $\ell$th element of $I_{r}$, and define the
EEE, denoted by $\widehat{G}_{J}(u)$, to be the empirical version of
$G_{J}(u)$ with $\xi_{1}$, $\mathbb{E}[g_{1}^{3}(X_{1})]$ and
$\mathbb{E}[g_{1}(X_{1})g_{1}(X_{2})g_{2}(X_{1},X_{2})]$ replaced by their
estimators. Define

$\displaystyle\widetilde{M}_{\alpha}(\rho_{n},n,\alpha;R):=$
$\displaystyle~{}\begin{cases}\rho_{n}^{-1}n^{-(\alpha-1)},&(r=2)\\\
\Big{\\{}\rho_{n}^{-s/2}n^{-(\alpha-1)/2}+\max_{p\in[3:(\lceil\alpha\rceil-1)]}\rho_{n}^{-v_{p}/2}n^{-(p-1)/2}\Big{\\}}\log^{1/2}n,&(r\geq
3)\end{cases}$
$\displaystyle+n^{-1}\cdot\log^{3/2}n+\begin{cases}(\rho_{n}\cdot
n)^{-1}\cdot\log^{1/2}n,&\textrm{ (Acyclic }R)\\\ \rho_{n}^{-r/2}\cdot
n^{-1}\cdot\log^{1/2}n,&\textrm{ (Cyclic }R)\end{cases},$ (55)

where $v_{p}$ is the maximum number of edges among all $p$-node subgraphs of
$R$. We have

######  Theorem 4.1.

Suppose the network is generated by the graphon model. Set $\alpha\in(1,r]$
satisfying $\widetilde{M}_{\alpha}(\rho_{n},n,\alpha;R)\to 0$. Also assume the
following conditions:

  1. 1.

$\rho_{n}^{-2s}\cdot\xi_{1}^{2}\geq{\rm constant}>0$.

  2. 2.

The design ${\cal J}_{n,\alpha}$ satisfies Assumptions 1 and 2.

Then we have

$\displaystyle\big{\|}F_{\widehat{T}_{J}(u)}-G_{J}(u)\big{\|}_{\infty}=$
$\displaystyle~{}O(\widetilde{\cal M}_{\alpha}(\rho_{n},n,\alpha;R)),$ (56)
$\displaystyle\big{\|}F_{\widehat{T}_{J}(u)}-\widehat{G}_{J}(u)\big{\|}_{\infty}=$
$\displaystyle~{}{\widetilde{O}_{p}}(\widetilde{\cal
M}_{\alpha}(\rho_{n},n,\alpha;R)).$ (57)

Theorem 4.1 provides the first rigorous description of the trade-off between
computational speed and the accuracy of risk control in reduced network
moments, reflected by the first term in $\widetilde{\cal
M}_{\alpha}(\rho_{n},n,\alpha;R)$. Theorem 4.1 sharply contrasts Theorem 2.1.
For noiseless U-statistics, we can achieve higher-order accurate risk control
by setting $\alpha>1$, and increasing $\alpha$ beyond 2 will not further
benefit this accuracy; but for network moments with $r\geq 3$ edges (i.e.
except $R=$ Edge), even in the densest network setting $\rho_{n}\asymp 1$,
higher-order accurate risk control still requires $\alpha>2$ regardless of the
motif’s shape. This distinction is caused by edge-wise observational errors. A
close inspection of the result of Theorem 4.1 also finds an interesting fact
that for completely dense networks, increasing computational budget beyond
$O(n^{3})$ will not further improve inference accuracy, analogous to its
sibling result for noiseless U-statistics. Comparing the results of Theorem
4.1 and Theorem 2.1 also raises another interesting question: if we
incorporate some $O(n^{-(\alpha-1)\ell})$ terms, like those in (13), into
(52), will that improve the error bound in Theorem 4.1? The answer is no,
because the first term in $\widetilde{\cal M}_{\alpha}(\rho_{n},n,\alpha;R)$
is _not_ caused by ignoring those terms in the Edgeworth expansion, but
instead due to edge-wise observational errors.

We conclude this section by studying the case $\rho_{n}^{s}n^{\alpha-1}\to 0$,
either due to network sparsity, or because we reduce the computation too much.
Now, the randomness in $\widehat{U}_{J}$ is dominated by the randomness in
$\widehat{U}_{J}-U_{J}$. Not surprisingly, the variance estimator for
$\widehat{U}_{J}$ needs some revision. With the new variance estimator
$(\widehat{\sigma}_{J}^{\prime})^{2}:=|{\cal
J}_{n,\alpha}|^{-1}\widehat{U}_{J}$, the studentization form is

$\displaystyle\widehat{T}_{J}^{\prime}:=$
$\displaystyle~{}\frac{\widehat{U}_{J}-\mu}{\widehat{\sigma}_{J}^{\prime}}.$
(58)

The next theorem establishes the asymptotic normality of
$\widehat{T}_{J}^{\prime}$. Define
$\widecheck{M}_{\alpha}(\rho_{n},n,\alpha;R)$ to almost equal $\widetilde{\cal
M}_{\alpha}(\rho_{n},n,\alpha;R)$, except for replacing
$\rho_{n}^{-s/2}n^{-(\alpha-1)/2}$ by $\rho_{n}^{-s/2}n^{-\alpha/2}$.

######  Theorem 4.2.

Assume the conditions of Theorem 4.1 hold, except replacing the condition
$\widetilde{\cal M}_{\alpha}(\rho_{n},n,\alpha;R)\to 0$ by
$\widecheck{M}_{\alpha}(\rho_{n},n,\alpha;R)\to 0$. Then we have

$\displaystyle\widehat{T}_{J}^{\prime}\stackrel{{\scriptstyle
d}}{{\to}}N(0,1).$ (59)

In Theorem 4.2, we still need the assumption
$\rho_{n}^{s}n^{\alpha}\to\infty$, which is not much weaker than
$\rho_{n}^{s}n^{\alpha-1}\to\infty$, but based on the empirical evidence from
our simulation result (Section 5.3, Figure 4), we conjecture that the
assumption $\rho_{n}^{s}n^{\alpha}\to\infty$ should be minimal131313As a
theoretical evidence: consider the simplest motif: $R=$ Edge. When
$\rho_{n}\cdot n^{\alpha}\to\lambda$ (constant), the limiting distribution of
$\widehat{T}_{J}^{\prime}$, after rescaling, is Poisson, not normal (Novak,
2019)..

##  5 Simulations

###  5.1 Simulation 1: non-degenerate U-statistics

Our first simulation assesses CDF approximation accuracy for noiseless non-
degenerate U-statistics. The goal is to accurately approximate
$F_{T_{J}+\delta_{J}}$, where we set a small variance with $C_{\delta}=0.008$
for $\delta_{J}$. We generated synthetic data with
$X_{1},\ldots,X_{n}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}$ PDF:
$(x+1)/2$, $x\in[-1,1]$; and the kernel function
$h(x_{1},x_{2},x_{3}):=\sin(x_{1}+x_{2}+x_{3})$. We experimented our proposed
deterministic design in Section 2.2.1 and a random design (J1) from Section
2.2.2. We compared our method to the following benchmarks: 1\. $N(0,1)$; 2.
Resample bootstrap (bootstrap iteration $B=200$ (Levin and Levina, 2019)); 3.
Subsample bootstrap (subsample size: $n^{1/2}$). We emulate the true sampling
distribution of $T_{J}+\delta_{J}$ by a Monte-Carlo approximation with samples
$n_{MC}:=10^{6}$. The performance measure is

$\displaystyle\sup_{u\in[-2,2];10u\in\mathbb{Z}}\big{|}\widehat{F}_{T_{J}+\delta_{J}}(u)-F_{T_{J}+\delta_{J}}(u)\big{|}.$
(60)

We varied $n\in\\{10,20,40,80\\}$ and set $\alpha=1.5$ (result for
$\alpha=1.7$ in Supplementary Material). For each $(n,\alpha)$ setting, we
repeated the experiment 30 times and recorded the mean and standard deviation
of the distribution approximation errors (60).

Figure 1:  Non-degenerate U-statistics: Row 1: deterministic ${\cal
J}_{n,\alpha}$ (Section 2.2); row 2: random ${\cal J}_{n,\alpha}$ (scheme
(J1)). Column 1: true CDF $=F_{T_{J}+\delta_{J}}(u)$, $n=80$; column 2: CDF
approximation error; column 3: CI coverage probability, dashed blue line
$=1-\beta=90\%$; column 4: CI length; column 5: computation time.

Column 1 of Figure 1 shows the true and estimated CDF curves for
$T_{J}+\delta_{J}$. We can see that our method’s estimated CDF curve almost
overlaps the true CDF curve; whereas other methods all exhibit noticeable
deviations at different levels. Then Column 2 of Figure 1 shows log-
transformed CDF approximation errors for different methods under different
$(n,\alpha)$ configurations. Our method shows clear advantage in accuracy
across all settings. Specifically, our method is the only method with an
empirical error rate faster than $n^{-1/2}$. This well-matches our theoretical
prediction on the higher-order accuracy of our method.

Next, we compare our Cornish-Fisher confidence interval to the three benchmark
methods in Simulation 1, plus the C-F CI constructed based on the complete
U-statistic, in the following aspects: coverage probability, CI length and
computation time. We fixed confidence level $1-\beta=90\%$ and focused on two-
sided confidence intervals for simplicity. Most simulation settings were
inherited from Simulation 1, but now we do not need a large $n_{MC}$ and can
afford larger $n$’s, so we varied $n\in\\{25,50,100,200,400\\}$. In each
experiment, to evaluate one empirical CI coverage probability, we generate
3000 CI’s for our method, $N(0,1)$ and complete U-statistics; and 500 CI’s for
resampling and subsampling bootstraps since they take longer to run; and the
experiment is repeated 100 times for all methods except the complete
U-statistic method (repeated 20 times).

Figure 1 shows the result for deterministic and random designs. Our method
shows clear advantage in accuracy of controlling the empirical coverage
probability around the nominal $90\%$ level, significantly improves over
normal approximation, especially for small $n$’s. As $n$ grows large, our
method’s speed advantage over bootstrap methods becomes clearer. Compared to
inference based on complete U-statistic, our method effectively reduces
computational complexity (much flatter computation time curve) without
noticeable loss in risk control accuracy. All methods (except subsampling
bootstrap) produce similar CI lengths. This echoes our earlier remarks that
the CI length reflects a different aspect of U-statistic reduction (inference
power, Section 1.3); and different approaches may perform similarly in this
aspect, if they are all asymptotically normal approximations.

###  5.2 Simulation 2: degenerate U-statistics

We generate $X_{1},\ldots,X_{n}\stackrel{{\scriptstyle\rm i.i.d.}}{{\sim}}$
Uniform$(0,1)$ and with the kernel
$h(x_{1},x_{2},x_{3}):=27\prod_{\ell=1}^{3}\big{\\{}\sin(2x_{\ell})-\sin^{2}1\big{\\}}$.
We set $n\in\\{10,20,\ldots,160\\}$. $\alpha=1.5$ and $n_{\rm MC}=10^{6}$, and
repeat each individual experiment 30 times. In Figure 2, plots 1 and 2 confirm
our method’s superiority in accurately approximating
$F_{\frac{U_{J}-U_{n}}{|{\cal
J}_{n,\alpha}|^{-1/2}\cdot\widetilde{\sigma}_{h}}+\delta_{J}\big{|}X_{[1:n]}}(u)$.
Next, we study the very interesting observation of our Theorem 3.2 and in
earlier literature Weber (1981); Chen and Kato (2019a) that the randomness and
incompleteness of a randomly designed ${\cal J}_{n,\alpha}$ reinstates
normality in $U_{J}$. To this end, we set
$\alpha\in\\{1,1.25,1.5,\ldots,2.75\\}$, $n\in\\{50,100,200,400\\}$, and set
$n_{\rm MC}=500$ due to the high computational cost at $\alpha=2.75$. Plot 3
in Figure 2 confirms the prediction of our Theorem 3.2 that the best $\alpha$
in this example should be close to $k_{0}/3=1$, despite setting $\alpha$
slightly above 1 noticeably improves the accuracy of CI level control,
possibly due to numerical reasons; whereas setting $\alpha$ close very close
to $k_{0}=3$ (e.g. $\alpha=2.75$) would significantly depreciate the accuracy
of our method, since the distribution of $T_{J}$ at $\alpha=k_{0}\geq 2$ is
non-Gaussian. The last two plots in Figure 2 confirm our method’s fine
performance in CI length (which reflects power) and computation speed.

Figure 2:  Degenerate U-statistics: Column 1: CDF curves; column 2: CDF
approximation error; columns 3–5: CI performance of our method with various
$\alpha$’s.

###  5.3 Simulation 3: network moments

We first compare the CDF approximation accuracy of different methods for dense
networks generated by the “BlockModel” in Zhang and Xia (2022) with
$\rho_{n}\asymp 1$. We experimented the four graphons (Edge, V-shape, Triangle
and Three-Star) in Zhang and Xia (2022), plus a new and more complex graphon,
which we call “Out-Triangle” (Biemann et al., 2012), defined by its edge set
$\\{(1,2),(2,3),(1,3),(1,4)\\}$. By Theorem 4.1, $R=$ Edge follows a different
guideline for choosing $\alpha$. Due to page limit, we sink the study of the
edge motif to Supplementary Material. For the other four motifs, we set
$\alpha=2.5$, $n\in\\{25,50,100,200\\}$, and inherit other simulation settings
from Section 5.1. Figure 15 shows that our method maintains stable advantage
over competing methods in CDF approximation accuracy. Next, Table 2 presents
detailed numerical comparison of the CI’s produced by different methods. We
see that our method achieves a clear acceleration over the complete
U-statistic method with only slightly depreciated confidence level control.
Our method’s speed is competitive even compared to the normal approximation.
It is significantly faster and more accurate than subsampling bootstrap, and
much faster than resampling bootstrap of comparable accuracy.

Figure 3:  Network moments: Row 1: CDF curves, true CDF
$=F_{\widehat{T}_{J}+\delta_{J}}(u)$, $n=100$. Our method: $\alpha=2.5$,
random ${\cal J}_{n,\alpha}$ (scheme (J1)). Row 2: CDF approximation
errors151515The Monte Carlo evaluation of the true CDF timed out for Out-
Triangle when $n=200$..  Table 2: Network moments: performance of $95\%$ CI’s,
$n=100$, BlockModel, Mean(SD)

Motif V-shape Triangle ThreeStar OutTriangle Our method ($\alpha=2.5$)
($n_{\rm MC}=10^{4})$ Coverage $=0.952(0.215)$ Length $=0.173(0.025)$ LogTime
$=-4.390(0.013)$ $0.949(0.219)$ $0.035(0.006)$ $-4.604(0.035)$ $0.945(0.228)$
$0.124(0.024)$ $-3.893(0.015)$ $0.940(0.238)$ $0.200(0.049)$ $-2.919(0.020)$
Norm. Approx. ($\alpha=2.5$) ($n_{\rm MC}=10^{4})$ $0.939(0.238)$
$0.173(0.025)$ $-4.549(0.072)$ $0.934(0.249)$ $0.035(0.006)$ $-4.869(0.019)$
$0.930(0.255)$ $0.124(0.024)$ $-4.062(0.018)$ $0.922(0.268)$ $0.200(0.049)$
$-3.058(0.010)$ Complete ($n_{\rm MC}=500)$ $0.953(0.211)$ $0.172(0.025)$
$-4.114(0.078)$ $0.951(0.215)$ $0.034(0.006)$ $-4.522(0.018)$ $0.948(0.223)$
$0.123(0.024)$ $-0.727(0.016)$ $0.942(0.234)$ $0.199(0.048)$ $0.272(0.026)$
Resample ($\alpha=2.5$) ($n_{\rm MC}=500)$ $0.947(0.224)$ $0.179(0.030)$
$0.756(0.010)$ $0.952(0.214)$ $0.038(0.007)$ $0.559(0.014)$ $0.947(0.223)$
$0.131(0.028)$ $1.266(0.012)$ $0.955(0.208)$ $0.228(0.057)$ $2.356(0.012)$
Subsample ($\alpha=2.5$) ($n_{\rm MC}=500)$ $0.964(0.187)$ $0.348(0.082)$
$-2.716(0.015)$ $0.923(0.267)$ $0.071(0.026)$ $-2.791(0.014)$ $0.959(0.198)$
$0.548(0.193)$ $-2.301(0.025)$ $0.929(0.256)$ $0.979(0.445)$ $-1.787(0.014)$

Our last experiment numerically pokes the asymptotic normality of
$\widehat{T}_{J}^{\prime}$ for sparse networks under different amounts of
computation reduction. Here, we set $R=$ Triangle and the graphon to be
BlockModel. We tested two sparsity regimes: $\rho_{n}\asymp n^{-1/2}$ and
$n^{-2/3}$. Our Theorem 4.2 predicts $\widehat{T}_{J}^{\prime}\to N(0,1)$ when
$\rho_{n}^{s}n^{\alpha-1}\to 0$ and $\rho_{n}^{s}n^{\alpha}\to\infty$,
reflected by our simulation settings $\alpha=2$ for $\rho_{n}\asymp n^{-1/2}$
and $\alpha=2.5$ for $\rho_{n}\asymp n^{-2/3}$, respectively. For each
$\rho_{n}$, we also tested smaller $\alpha$’s. We repeated each configuration
$n_{\rm MC}=500$ times and compared the sampling distribution of
$\widehat{T}_{J}^{\prime}$ to $N(0,1)$. The QQ plots in Figure 4.2 well-
matches our Theorem 4.2’s prediction. In plots 2 & 4, we clearly see the
expected convergence to $N(0,1)$. In contrast, in plots 1 & 3, where
$\rho_{n}^{s}n^{\alpha}\to 0$, we see non-Gaussian behaviors, suggesting that
the conditions of our Theorem 4.2 seem minimal.

Figure 4:  Network moments (sparse networks): QQ plot of sampling distribution
of $\widehat{T}_{J}^{\prime}$, Triangle, BlockModel. Columns 1 & 2:
$\rho_{n}\asymp n^{-1/2}$; columns 3 & 4: $\rho_{n}\asymp n^{-2/3}$. Due to
memory limit, we did not run $n\geq 5120$ for $\alpha=2.5$.

##  6 Data examples

###  6.1 Data example 1: Stock market data

The S&P 500 historical data ((Datahub), ) records the daily prices of 412
stocks from 11 sectors. Following Chakraborty and Zhang (2021), we computed
the _monthly logarithmic return rates_ of each stock from 1-March-2000 to
29-August-2022, yielding $n=138$ observations. The goal is to assess the
pairwise dependency between all sectors using an independence test. Denote the
stock returns from a sector $X$ containing $p$ stocks as
$S_{t}^{X}=(S_{it}^{X},...,S_{pt}^{X})$, $t\in[1:138]$; similarly define
$S_{t}^{Y}$. We measure dependency by dCov, rewritten as a U-statistic (Lemma
1 of Yao et al. (2018)).

$\displaystyle{\rm dCov}^{2}(X,Y):=$
$\displaystyle~{}\binom{n}{4}^{-1}\sum_{i<j<q<r}h(Z_{i},Z_{j},Z_{q},Z_{r}),$
(61)

where
$h(Z_{i},Z_{j},Z_{q},Z_{r}):=\sum_{s,t,u,v}^{i,j,q,r}(a_{st}b_{uv}+a_{st}b_{st}-a_{st}b_{su}-a_{st}b_{tv})\big{/}24$
161616The summation notation “$\sum_{s,t,u,v}^{i,j,q,r}$” means
$\sum_{(s,t,u,v)\in{\rm All.Permutations}(i,j,q,r)}$.,
$a_{ij}=\|S_{i}^{X}-S_{j}^{X}\|_{2}$, and
$b_{ij}=\|S_{i}^{Y}-S_{j}^{Y}\|_{2}$. We set $\alpha=1.5$ and used it to
evaluate the test statistic for $H_{0}:\mathbb{E}[{\rm dCov}^{2}]=0$ between
each sector pair. Observing that $U_{J}$’s are on the same order as
$\widehat{\xi}_{1}$, we used our formula for the non-degenerate case. As a
reference, on the diagonal, we also randomly split the stocks in each sector
into two sets and evaluated their dependency using the same method. Figure 5
shows that our inference procedure for reduced U-statistics well-preserves the
result of the complete U-statistic but computes much faster (Table 3). On the
diagonal, the sectors that exhibit strongest inner dependency include _CD, E,
F, I_ and _IT_. This is understandable since they tend to be more sensitive to
global economic changes. In contrast, members of _CmS, CnS_ and _U_ sectors
focus more on local markets, so their price fluctuations are less
synchronized. This understanding also applies to cross-sector relations, such
as the tight connection between the pairs (_CD_ , _I_) and (_I_ , _IT_),
whereas _U_ is comparatively less dependent on other sectors except _E_.

Figure 5:  Pairwise dependency test: heatmaps of test statistics. High values
(red): high detected dependency. Each cell reports mean(std.) of test
statistics over 30 repeated experiments, except the off-diagonal of complete
U-statistic method (no repetition needed).

###  6.2 Data example 2: UCR time series data (Earthquakes, Starlight)

In the second example, we analyze two UCR time series data sets (Dau et al.,
2018): _Earthquakes_ and _Starlight_. The former records $n=461$ earthquake
curves of length $T=512$. All curves are classified into $K=2$ types:
$n_{0}=368$ non-major and $n_{1}=93$ major ones. Following Chakraborty and
Zhang (2021) and Zhu and Shao (2021), we treat each earthquake curve as a data
point in a Hilbert space and aim at comparing the population distributions of
the curves of different types using _Maximum Mean Discrepancy (MMD)_ metric.
We measure the distance between two earthquake curves by comparing their SRVF
transforms (Srivastava et al., 2011) – a method that synchronizes the two
curves’ phases in presence of amplitude discrepancy. The computation of SRVF
for each curve pair is slow (Strait et al., 2019), so we reduced each curve
$\\{x_{t}\\}_{t=1}^{512}$ into $\big{\\{}\widetilde{x}_{t}:={\rm
Mean}\big{(}x_{[\\{t-(\ell-1)/2\\}:\\{t+(\ell-1)/2\\}]}\big{)}\big{\\}}_{t\in\\{4k+1,k\in[0:127]\\}}$.
Due to page limit, here we only present $\ell=7$; results for more window
sizes can be found in Supplementary Material. We set $\alpha=1.5$ and apply
our method to reduce the MMD, rewritten as a U-statistic according to Equation
(6) in Schrab et al. (2021) with the RBF kernel $k(x,y):=\exp(-{\rm
SRVF}(x,y)^{2}/5000)$. To apply this MMD formula, we subsampled the larger
class (small earthquakes) to equate the two sample sizes. We also computed the
averaged pairwise distance (using SRVF) within each cluster to illustrate its
inner cohesion. Row 1 of Figure 6 shows the result. We certainly do not have
access to the population MMD value, but deeming the empirical (complete) two-
sample MMD as a good approximation, our Cornish-Fisher confidence intervals
with randomized design (J1) shows good coverage. The same fine performance
shows up for the within-cluster average pairwise distances.

Next, we apply the same analysis to the much larger _Starlight_ data set that
contains 3 types of stars ($n_{1}=1329$, $n_{2}=2580$ and $n_{3}=5327$). Each
curve is recorded as a length $1024$ sequence. We down-sampled each curve
sequence to length $128$ without smoothing, because the starlight curves are
much smoother than the earthquake curves. Even so, evaluating a complete
U-statistic for comparing any two star types remains computationally
infeasible, due to the large sample sizes. Our method with $\alpha=1.5$ makes
it possible to implement a reduced version of Equation (6) in Schrab et al.
(2021) with the RBF kernel $k(x,y):=\exp(-{\rm SRVF}(x,y)^{2}/100)$. Due to
page limit, in row 2 of Figure 6, we only present the result for the
comparison between type 1 and type 2 stars. Compared to the _Earthquakes_
data, the MMD CI’s for _Starlight_ are much stabler, understandably due to the
latter’s much larger sample sizes of each cluster. Also, some MMD CI’s of the
_Earthquakes_ data contain 0 (will not reject $H_{0}$), while all CI’s for the
_Starlight_ data clearly support a two-sided alternative. This is not
surprising, given the much smaller within-cluster distance (column 4) and the
clearer between-cluster differences (compare columns 1 & 2) in the latter
data.

Figure 6:  Results of data example 2. Row 1: _Earthquakes_ ; row 2:
_Starlight_. Column 1 & 2: raw and pre-pocessed curves; column 3: 90% CI of
based on reduced between-cluster MMD; column 4: 90% CI of within-cluster
average pairwise distance (using SRVF (Srivastava et al., 2011)). Dashed line:
complete U-statistic (evaluations of complete U-statistics timed out ($>48$
hours) in most settings for _Starlight_ data).  Table 3: Time cost: our method
($\alpha=1.5$) vs. complete U-statistic

Time cost Stock Market ($r=4$) Earthquakes ($r=2$) (Unit = sec.) All Major
Non-major Maj. vs. Non-Maj. Our method 3.47 303.94 2471.70 1223.50 Complete U
8099.73 708.99 11199.92 17912.91 Time cost Starlight ($r=2$) (Unit = sec.)
Type 1 Type 2 Type 3 1 vs. 2 1 vs. 3 2 vs. 3 Our method 4512.95 12773.76
41282.26 19140.13 19149.33 50413.75 Complete U 48227.72 158233.7 (Time out)
(Time out) (Time out) (Time out)

##  7 Discussion

Our work assumes that the degree of degeneracy $k_{0}$ is known. This is a
common practice in the analysis of degenerate U-statistics and seemingly
cannot be circumvented by bootstraps: as Remark 3.2 in Chen and Kato (2019a)
pointed out, the knowledge of $k_{0}$ is critical for correctly setting how
many $r$-tuples to be randomly sampled. Fortunately, our proofs imply that
each $\xi_{k}$’s can be estimated by $\widehat{\xi}_{k}$ defined in Section 2
with concentration bound
$|\widehat{\xi}_{k}-\xi_{k}|={\widetilde{O}_{p}}(n^{-1/2}\log^{k/2}n)$.
Therefore, one can simply test $H_{0}(k):\xi_{k}=0$ versus
$H_{a}(k):\xi_{k}={\rm constant}>0$ by comparing $\widehat{\xi}_{k}$ to
$n^{-c_{0}}$, $c_{0}\in(0,1/2)$ to find $k_{0}$. Another approach, called the
“blocking trick”, by Ho and Shieh (2006), breaks the data $X_{[1:n]}$ into
$K_{0}$ blocks $X_{[\\{(k-1)(n/K_{0})+1\\}:(kn/K_{0})]},k\in[1:K_{0}]$,
computes a complete U-statistic on each block, and makes inference using a CLT
for the sample mean of $K_{0}$ i.i.d. terms. Its unification of degenerate and
non-degenerate cases and automatic adaptation to $k_{0}$, however, comes at
the price of heavily wasting data.

Our study throughout this paper exclusively focuses on _data-oblivious_
reduction schemes. Recently, Kong and Zheng (2021) proposed a _data-aware_
reduction scheme, based on their key observation that $X_{[1:r]}\approx
Y_{[1:r]}$ implies $h(X_{[1:r]})\approx h(Y_{[1:r]})$, thus by clustering
$X_{i}$’s, one can effectively reduce the U-statistic’s computation. While
their method shows very attractive performance, finite-sample higher-order
analysis for their method poses an interesting open challenge. There is also a
computational price for _being data-aware_. For example, in the setting
considered by Moon and Chen (2022), the clustering of all $X_{i}$’s in some
Banach space requires computing at least $O(n^{2})$ many potentially expensive
(like in our second data example) pairwise distances.

## References

  * Aldous (1981) Aldous, D. J. (1981).  Representations for partially exchangeable arrays of random variables.  Journal of Multivariate Analysis 11(4), 581–598.
  * Banerjee and Ma (2017) Banerjee, D. and Z. Ma (2017).  Optimal hypothesis testing for stochastic block models with growing degrees.  arXiv preprint arXiv:1705.05305.
  * Bergsma and Dassios (2014) Bergsma, W. and A. Dassios (2014).  A consistent test of independence based on a sign covariance related to Kendall’s tau.  Bernoulli 20(2), 1006–1028.
  * Bickel et al. (2011) Bickel, P. J., A. Chen, and E. Levina (2011).  The method of moments and degree distributions for network models.  The Annals of Statistics 39(5), 2280–2301.
  * Biemann et al. (2012) Biemann, C., S. Roos, and K. Weihe (2012).  Quantifying semantics using complex network analysis.  In Proceedings of COLING 2012, pp. 263–278.
  * Blom (1976) Blom, G. (1976).  Some properties of incomplete U-statistics.  Biometrika 63(3), 573–580.
  * Brown and Kildea (1978) Brown, B. and D. Kildea (1978).  Reduced U-statistics and the Hodges-Lehmann estimator.  The Annals of Statistics, 828–835.
  * Chakraborty and Zhang (2021) Chakraborty, S. and X. Zhang (2021).  A new framework for distance and kernel-based metrics in high dimensions.  Electronic Journal of Statistics 15(2), 5455–5522.
  * Chaudhuri and Hu (2019) Chaudhuri, A. and W. Hu (2019).  A fast algorithm for computing distance correlation.  Computational Statistics & Data Analysis 135, 15–24.
  * Chen and Peng (2021) Chen, S. X. and L. Peng (2021).  Distributed statistical inference for massive data.  The Annals of Statistics 49(5), 2851–2869.
  * Chen and Kato (2019a) Chen, X. and K. Kato (2019a).  Randomized incomplete U-statistics in high dimensions.  The Annals of Statistics 47(6), 3127–3156.
  * Chen and Kato (2019b) Chen, X. and K. Kato (2019b).  Supplementary Material to “Randomized incomplete U-statistics in high dimensions”.  The Annals of Statistics.
  * Clémençon et al. (2016) Clémençon, S., I. Colin, and A. Bellet (2016).  Scaling-up empirical risk minimization: optimization of incomplete U-statistics.  The Journal of Machine Learning Research 17(1), 2682–2717.
  * (14) (Datahub).  S&P 500 companies with financial information.  https://datahub.io/core/s-and-p-500-companies#data-cli.  Accessed 18-Nov-2022.
  * Dau et al. (2018) Dau, H. A., E. Keogh, K. Kamgar, C.-C. M. Yeh, Y. Zhu, S. Gharghabi, C. A. Ratanamahatana, Yanping, B. Hu, N. Begum, A. Bagnall, A. Mueen, G. Batista, and Hexagon-ML (2018, October).  The UCR time series classification archive.  https://www.cs.ucr.edu/~eamonn/time_series_data_2018/.
  * Dürre and Paindaveine (2021) Dürre, A. and D. Paindaveine (2021).  On the consistency of incomplete U-statistics under infinite second-order moments.  arXiv preprint arXiv:2112.14666.
  * Even-Zohar and Leng (2021) Even-Zohar, C. and C. Leng (2021).  Counting small permutation patterns.  In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 2288–2302. SIAM.
  * Gao and Lafferty (2017) Gao, C. and J. Lafferty (2017).  Testing network structure using relations between small subgraph probabilities.  arXiv preprint arXiv:1704.06742.
  * Gao et al. (2015) Gao, C., Y. Lu, and H. H. Zhou (2015).  Rate-optimal graphon estimation.  The Annals of Statistics 43(6), 2624–2652.
  * Giraudo (2021) Giraudo, D. (2021).  Limit theorems for U-statistics of Bernoulli data.  Latin American Journal of Probability and Mathematical Statistics 18, 793–828.
  * Gretton et al. (2012) Gretton, A., K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola (2012).  A kernel two-sample test.  The Journal of Machine Learning Research 13(1), 723–773.
  * Hall (1983) Hall, P. (1983).  Inverting an Edgeworth expansion.  The Annals of Statistics 11(2), 569–576.
  * Hall (2013) Hall, P. (2013).  The Bootstrap and Edgeworth Expansion.  Springer Science & Business Media.
  * Han and Qian (2018) Han, F. and T. Qian (2018).  On inference validity of weighted U-statistics under data heterogeneity.  Electronic Journal of Statistics 12(2), 2637–2708.
  * Heller and Heller (2016) Heller, Y. and R. Heller (2016).  Computing the bergsma dassios sign-covariance.  arXiv preprint arXiv:1605.08732.
  * Helmers (1991) Helmers, R. (1991).  On the Edgeworth expansion and the bootstrap approximation for a studentized U-statistic.  The Annals of Statistics 19(1), 470–484.
  * Ho and Shieh (2006) Ho, H.-C. and G. S. Shieh (2006).  Two-stage U-statistics for hypothesis testing.  Scandinavian Journal of Statistics 33(4), 861–873.
  * Hoover (1979) Hoover, D. N. (1979).  Relations on probability spaces and arrays of random variables.  Preprint, Institute for Advanced Study, Princeton, NJ 2.
  * Huo and Székely (2016) Huo, X. and G. J. Székely (2016).  Fast computing for distance covariance.  Technometrics 58(4), 435–447.
  * Janson (1984) Janson, S. (1984).  The asymptotic distributions of incomplete U-statistics.  Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 66(4), 495–505.
  * Kong and Zheng (2021) Kong, X. and W. Zheng (2021).  Design based incomplete U-statistics.  Statistica Sinica 31(3), 1593–1618.
  * Lahiri (1993) Lahiri, S. N. (1993).  Bootstrapping the studentized sample mean of lattice variables.  Journal of Multivariate Analysis 45(2), 247–256.
  * Lee (1979) Lee, A. J. (1979).  On the asymptotic distribution of certain incomplete U-statistics.  Technical report, North Carolina State University. Dept. of Statistics.
  * Lee (1982) Lee, A. J. (1982).  On incomplete U-statistics having minimum variance.  Australian Journal of Statistics 24(3), 275–282.
  * Lee (2019) Lee, A. J. (2019).  U-Statistics: Theory and Practice.  Routledge.
  * Levin and Levina (2019) Levin, K. and E. Levina (2019).  Bootstrapping networks with latent space structure.  arXiv preprint arXiv:1907.10821.
  * Maesono (1997) Maesono, Y. (1997).  Edgeworth expansions of a studentized U-statistic and a jackknife estimator of variance.  Journal of Statistical Planning and Inference 61(1), 61–84.
  * Major (2007) Major, P. (2007).  On a multivariate version of Bernstein’s inequality.  Electronic Journal of Probability 12, 966–988.
  * Moon and Chen (2022) Moon, H. and K. Chen (2022).  Interpoint-ranking sign covariance for the test of independence.  Biometrika 109(1), 165–179.
  * Novak (2019) Novak, S. (2019).  Poisson approximation.  Probability Surveys 16, 228–276.
  * Putter and van Zwet (1998) Putter, H. and W. R. van Zwet (1998).  Empirical edgeworth expansions for symmetric statistics.  The Annals of Statistics 26(4), 1540–1569.
  * Rempala and Wesolowski (2003) Rempala, G. and J. Wesolowski (2003).  Incomplete U-statistics of permanent design.  Journal of Nonparametric Statistics 15(2), 221–236.
  * Rosenbaum (2011) Rosenbaum, P. R. (2011).  A new U-statistic with superior design sensitivity in matched observational studies.  Biometrics 67(3), 1017–1027.
  * Schrab et al. (2021) Schrab, A., I. Kim, M. Albert, B. Laurent, B. Guedj, and A. Gretton (2021).  MMD aggregated two-sample test.  arXiv preprint arXiv:2110.15073.
  * Shao et al. (2022) Shao, M., D. Xia, Y. Zhang, Q. Wu, and S. Chen (2022).  Higher-order accurate two-sample network inference and network hashing.  arXiv preprint arXiv:2208.07573.
  * Slater (1960) Slater, L. J. (1960).  Confluent Hypergeometric Functions.  Cambridge University Press.
  * Srivastava et al. (2011) Srivastava, A., W. Wu, S. Kurtek, E. Klassen, and J. S. Marron (2011).  Registration of functional data using fisher-rao metric.  arXiv preprint arXiv:1103.3817.
  * Strait et al. (2019) Strait, J., O. Chkrebtii, and S. Kurtek (2019).  Automatic detection and uncertainty quantification of landmarks on elastic curves.  Journal of the American Statistical Association 114(527), 1002–1017.
  * Székely et al. (2007) Székely, G. J., M. L. Rizzo, and N. K. Bakirov (2007).  Measuring and testing dependence by correlation of distances.  The Annals of Statistics 35(6), 2769–2794.
  * Van der Vaart (2000) Van der Vaart, A. W. (2000).  Asymptotic Statistics, Volume 3.  Cambridge university press.
  * Weber (1981) Weber, N. (1981).  Incomplete degenerate U-statistics.  Scandinavian Journal of Statistics, 120–123.
  * Yao et al. (2018) Yao, S., X. Zhang, and X. Shao (2018).  Testing mutual independence in high dimension via distance covariance.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80(3), 455–480.
  * Zhang and Xia (2022) Zhang, Y. and D. Xia (2022).  Edgeworth expansions for network moments.  The Annals of Statistics 50(2), 726–753.
  * Zhao (2019) Zhao, Q. (2019).  On sensitivity value of pair-matched observational studies.  Journal of the American Statistical Association 114(526), 713–722.
  * Zhu and Shao (2021) Zhu, C. and X. Shao (2021).  Interpoint distance based two sample tests in high dimension.  Bernoulli 27(2), 1189–1211.