# Hoeffding-type decomposition for $U$-statistics on bipartite networks

Tâm Le Minh Sophie Donnet François Massol and Stéphane Robin

Hoeffding-type decomposition for $U$-statistics on bipartite networks

Tâm Le Minh 1 & Sophie Donnet 1 & François Massol 2 & Stéphane Robin 3

1 Université Paris-Saclay, AgroParisTech, INRAE, UMR MIA Paris-Saclay, 91120,
Palaiseau, France
2 Univ. Lille, CNRS, Inserm, CHU Lille, Institut Pasteur de Lille, U1019 - UMR
9017 - CIIL - Center for Infection and Immunity of Lille, F-59000 Lille,
France
3 Sorbonne Université, CNRS, Laboratoire de Probabilités, Statistique et
Modélisation, UMR 8001, 75005 Paris, France

#### Abstract.

We consider a broad class of random bipartite networks, the distribution of
which is invariant under permutation within each type of nodes. We are
interested in $U$-statistics defined on the adjacency matrix of such a
network, for which we define a new type of Hoeffding decomposition. This
decomposition enables us to characterize non-degenerate $U$-statistics – which
are then asymptotically normal – and provides us with a natural and easy-to-
implement estimator of their asymptotic variance.
We illustrate the use of this general approach on some typical random graph
models and use it to estimate or test some quantities characterizing the
topology of the associated network. We also assess the accuracy and the power
of the proposed estimates or tests, via a simulation study.

##  1 Introduction

Networks are popular objects to represent a set of interacting entities. The
last decades have witnessed an explosion in the number of networks datasets.
The fields of application range from sociology to ecology, from economics to
computer science. Understanding the organization of the network is a first
step towards a better insight of the system it represents. Several strategies
exist to study or describe the topology of a network. Many of them are based
on the calculation of one or several numeric quantities (statistics) such as
density, clustering coefficients or counts of given motifs to name but a few.
These statistics generally rely on several nodes.

The calculation of these numerical quantities on a given network naturally
leads to comparing them to a reference value, or to the value obtained on
another network. The concept of hypothesis tests naturally meets this
expectation. The challenging step of statistical hypothesis testing is
identifying the statistic distribution under the null assumption. In
particular, one class of statistic considered are the $U$-statistics which, in
the context of network analysis, have complex dependencies.

#### Networks and dissociated RCE matrices

In networks, entities are represented by nodes which are linked by edges when
they interact. In bipartite networks, the nodes are divided into two types and
the interactions only happen between nodes of the two different types. Some
examples of bipartite networks connect users and items in recommender systems
[53], papers and scientists in authorship networks [33], or plants and
pollinators in ecological interaction networks [16]. The networks are
naturally encoded in matrices. In the adjacency matrix $Y$ of a bipartite
network (sometimes also called incidence matrix), the two types of nodes are
represented by rows and columns, so that $Y_{ij}$ encodes the interaction
between entity $i$ of the first type and entity $j$ of the second type. In
binary networks, $Y_{ij}=1$ if $i$ and $j$ interact, else $Y_{ij}=0$. Some
networks are weighted, meaning that $Y_{ij}$ represents the intensity of the
interaction.

We consider the asymptotic framework where $Y$ is an infinite adjacency matrix
and the adjacency matrix of an observed network of size $m\times n$ is the
submatrix extracted from the leading $m$ rows and $n$ columns of $Y$.
Probabilistic models define a joint distribution on the values of the matrix
entries. In random graph models, it is common to assume that the nodes of the
networks are exchangeable, i.e. that the distribution of the network does not
change if its nodes are permuted. For instance, the stochastic blockmodel
[47], the random dot product graph model [51] or the latent space model [24]
are all node-exchangeable. On the corresponding adjacency matrix of a
bipartite network, this assumption implies the row-column exchangeability. $Y$
is said to be row-column exchangeable (RCE) if for any couple
$\Phi=(\sigma_{1},\sigma_{2})$ of finite permutations of $\mathbb{N}$,

$\Phi Y\overset{\mathcal{D}}{=}Y,$

where $\Phi Y:=(Y_{\sigma_{1}(i)\sigma_{2}(j)})_{i\geq 1,j\geq 1}$. Many
exchangeable random graph models also have a dissociatedness property, i.e.
their adjacency matrices are also dissociated [44, 28]. A RCE matrix is said
to be dissociated if for all $m$ and $n$, $(Y_{ij})_{1\leq i\leq m,1\leq j\leq
n}$ is independent from $(Y_{ij})_{i>m,j>n}$. In the present work, we only
consider RCE dissociated matrices.

####  $U$-statistics and Hoeffding decomposition

$U$-statistics are a generalization of the empirical mean to functions of more
than one variable. Many estimators fall under the category of $U$-statistics.
Given a sequence of random variables $Y=(Y_{i})_{i\geq 1}$ numbered with a
unique index a $U$-statistic $U^{h}_{n}(Y)$ of order $n$ and kernel function
$h$ is defined as the following average

$U^{h}_{n}(Y)=\binom{n}{k}^{-1}\sum_{1\leq i_{1}<i_{2}<...<i_{k}\leq
n}h(Y_{i_{1}},Y_{i_{2}},...,Y_{i_{k}}),$

where $h:\mathbb{R}^{k}\rightarrow\mathbb{R}$ is a symmetric function referred
to as the kernel. Denote $\llbracket n\rrbracket:=\\{1,...,n\\}$ and for a set
$A$, $\mathcal{P}_{k}(A)$ the set of all subsets of cardinal $k$ of $A$. Let
$\mathbf{i}=\\{i_{1},...,i_{k}\\}\in\mathcal{P}_{k}(\llbracket n\rrbracket)$,
then by symmetry of $h$, $h(Y_{i_{1}},...,Y_{i_{k}})$ does not depend on the
order of the elements of $\mathbf{i}$. Therefore, we will denote
$h(Y_{\mathbf{i}}):=h(Y_{i_{1}},...,Y_{i_{k}})$. Finally, the $U$-statistic
$U^{h}_{n}(Y)$ can be formulated as follows:

$U^{h}_{n}(Y)=\binom{n}{k}^{-1}\sum_{\mathbf{i}\in\mathcal{P}_{k}(\llbracket
n\rrbracket)}h(Y_{\mathbf{i}}).$

When $Y$ is an exchangeable sequence, $h(Y_{\mathbf{i}})$ has the same
distribution for all $\mathbf{i}\in\mathcal{P}_{k}(\llbracket n\rrbracket)$,
therefore $U^{h}_{n}(Y)$ is an unbiased estimate of $h(Y_{\llbracket
k\rrbracket})$. The case where $Y$ is an i.i.d. sequence is well-studied: the
$U$-statistics are known to be asymptotically normal [22] and can be used for
inference tasks such as estimation and hypothesis testing.

In the i.i.d. case, a useful technique to study the asymptotic behavior of
$U$-statistics is the Hoeffding decomposition, formalized for the first time
in [23]. For $1\leq c\leq k$, define the function
$\psi^{c}h:\mathbb{R}^{c}\rightarrow\mathbb{R}$ as

$\psi^{c}h:(y_{1},...,y_{c})\longmapsto\mathbb{E}[h(Y_{1},...,Y_{k})\mid
Y_{1}=y_{1},...,Y_{c}=y_{c}].$

Again, by symmetry of $h$, for some set
$\mathbf{i}\in\mathcal{P}_{c}(\llbracket n\rrbracket)$, we can denote
$\psi^{c}h(Y_{\mathbf{i}}):=\psi^{c}h(Y_{i_{1}},...,Y_{i_{c}})$ since the
order of the elements of $\mathbf{i}$ does not matter. Set
$p^{0}h=\mathbb{E}[h(Y_{\llbracket k\rrbracket})]$ and define recursively

$p^{c}h(Y_{\mathbf{i}})=\psi^{c}h(Y_{\mathbf{i}})-\sum_{c^{\prime}=0}^{c-1}\sum_{\mathbf{i}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{i})}p^{c^{\prime}}h(Y_{\mathbf{i}^{\prime}}).$

for all subsets $\mathbf{i}\in\mathcal{P}_{c}(\llbracket n\rrbracket)$, for
all $1\leq c\leq k$. Then, for $\mathbf{i}\in\mathcal{P}_{k}(\llbracket
n\rrbracket)$, $h(Y_{\mathbf{i}})$ can be written

$h(Y_{\mathbf{i}})=\sum_{0\leq c\leq
k}\sum_{\mathbf{i}^{\prime}\in\mathcal{P}_{c}(\mathbf{i})}p^{c}h(Y_{\mathbf{i}^{\prime}}).$

The $U$-statistic $U^{h}_{n}$ can be written as

$U^{h}_{n}(Y)=\sum_{c=0}^{k}\binom{k}{c}P^{c}_{n}h(Y),$

where for $1\leq c\leq k$,
$P^{c}_{n}h(Y)=\binom{n}{c}^{-1}\sum_{\mathbf{i}\in\mathcal{P}_{c}(\rrbracket
n\llbracket)}p^{c}h(Y_{\mathbf{i}})$.

This decomposition is interesting as all the quantities
$p^{c}h(Y_{\mathbf{i}})$ are orthogonal. By extension the $U$-statistics
$P^{c}_{n}h$ are also orthogonal. The leading terms of this decomposition have
been used by [22] to prove the asymptotic normality of $U$-statistics. The
decomposition also yields a decomposition of the variance of $U$-statistics.

By extension, a network $U$-statistic averages a function $h$ defined over
sub-matrices of size $p\times q$. Let $Y$ be an infinite adjacency matrix from
which we observe the first $m$ rows and $n$ columns. Let
$h:\mathcal{M}_{p,q}(\mathbb{R})\rightarrow\mathbb{R}$ be a function defined
on $p\times q$ matrices, $1\leq p\leq m$, $1\leq q\leq n$, verifying the
following symmetry property: for all
$(\sigma_{1},\sigma_{2})\in\mathbb{S}_{p}\times\mathbb{S}_{q}$,

$h(Y_{(i_{\sigma_{1}(1)},...,i_{\sigma_{1}(p)};j_{\sigma_{2}(1)},...,j_{\sigma_{2}(q))}})=h(Y_{(i_{1},i_{2},...,i_{p};j_{1},j_{2},...,j_{q})}),$
(1)

where $Y_{(i_{1},...,i_{p};j_{1},...,j_{q})}$ is the $p\times q$ sub-matrix
consisting of the rows and columns of $Y$ indexed by $i_{1},...,i_{p}$ and
$j_{1},...,j_{q}$ respectively. Therefore, since the order of the elements of
$\mathbf{i}=\\{i_{1},...,i_{p}\\}$ and $\mathbf{j}=\\{j_{1},...,j_{q}\\}$ does
not matter, we can denote
$h(Y_{\mathbf{i},\mathbf{j}}):=h(Y_{(i_{1},i_{2},...,i_{p};j_{1},j_{2},...,j_{q})})$.
Then the associated $U$-statistic is

$U^{h}_{m,n}(Y)=\binom{m}{p}^{-1}\binom{n}{q}^{-1}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{p}(\llbracket
m\rrbracket)\\\ \mathbf{j}\in\mathcal{P}_{q}(\llbracket
n\rrbracket)\end{subarray}}h(Y_{\mathbf{i},\mathbf{j}}).$ (2)

Note that the assumption on the symmetry of $h$ can be made without loss of
generality. Indeed, if
$h^{0}:\mathcal{M}_{p,q}(\mathbb{R})\rightarrow\mathbb{R}$ is not symmetric,
then $h:\mathcal{M}_{p,q}(\mathbb{R})\rightarrow\mathbb{R}$ defined by

$h(Y_{(i_{1},i_{2},...,i_{p};j_{1},j_{2},...,j_{q})})=(p!q!)^{-1}\sum_{(\sigma_{1},\sigma_{2})\in\mathbb{S}_{p}\times\mathbb{S}_{q}}h^{0}(Y_{(i_{\sigma_{1}(1)},...,i_{\sigma_{1}(p)};j_{\sigma_{2}(1)},...,j_{\sigma_{2}(q))}})$
(3)

verifies Equation (1) and leads to the same $U$-statistic
($U^{h^{0}}_{m,n}(Y)=U^{h}_{m,n}(Y)$).

In the case where $Y$ is a RCE matrix, [29] used a martingale approach to
obtain a weak convergence result for $U^{h}_{m,n}$ when $m$ and $n$ grow to
infinity at the same rate. Applying this result requires specific development
to get the asymptotic variance. In this paper, we propose an Hoeffding
decomposition-based approach. This strategy has the advantage to provide a
method to estimate the asymptotic variance of the $U$-statistics. Indeed,
obtaining an estimation of the (asymptotic) variance of $U$-statistics is a
required condition to perform practical inference such as hypothesis testing.
However, it remains a complex task and has been tackled with various methods
in the literature.

#### Variance estimation of $U$-statistics (related work)

The standard error of $U$-statistics is most often computed using resampling
techniques such as the jackknife [2] and bootstrap [17, 9] estimators of
variance. [41, 42, 11, 40] suggested various estimators of the asymptotic
variance. However, all these estimators are biased for both the variance and
the asymptotic variances. [11, 40] also discussed unbiased estimators, but
they are computationally more demanding than all the previous estimators and
they find them to have a positive probability of being negative, which is
undesirable.

#### Contribution

We show how $U$-statistics of size $p\times q$ can be used for exchangeable
network inference. In particular, we propose a Hoeffding-type decomposition to
identify the asymptotic distribution of these $U$-statistics and build a
consistent estimator of their variance. For that, we assume that the matrix
dimensions $m_{N}$ and $n_{N}$ grow at the same rate with
$m_{N}/N\rightarrow\rho$ and $n_{N}/N\rightarrow 1-\rho$, where
$\rho\in]0,1[$. For simplification, we denote
$U^{h}_{N}(Y):=U^{h}_{m_{N},n_{N}}(Y)$. First, we show that the distribution
of distribution $\sqrt{N}(U^{h}_{N}(Y)-\mathbb{E}[h])$ converges to a
Gaussian. Next, we build a consistent estimator for the variance of
$U^{h}_{N}(Y)$ and we study its properties.

#### Outline

In Section 2, we first extend the Hoeffding decomposition of a sub-graph
statistic. The decomposition we propose is based on the Aldous-Hoover-
Kallenberg (AHK) representation. Section 3 exploits this Hoeffding
decomposition to demonstrate the asymptotic normality of the $U$-statistics we
consider. In Section 4, we use the Hoeffding decomposition to build an
estimator of the asymptotic variance of the $U$-statistics. Section 5 presents
some models for RCE matrices and kernels that can be interesting for network
analysis and comparisons. Section 6 and 7 are dedicated to simulation studies
and an illustration of legislature dataset.

##  2 Hoeffding decomposition of a submatrix $U$-statistic

#### Aldous-Hoover-Kallenberg (AHK) representation

Corollary 7.23 of [25] states that for any dissociated RCE matrix $Y$, there
exists $(\xi_{i})_{i\geq 1}$, $(\eta_{j})_{j\geq 1}$ and
$(\zeta_{ij})_{i,j\geq 1}$ arrays of i.i.d. random variables with uniform
distribution over $[0,1]$ and a real measurable function $\phi$ such that for
all $1\leq i,j<\infty$,
$Y_{ij}\overset{a.s.}{=}\phi(\xi_{i},\eta_{j},\zeta_{ij})$. With such a
representation, the kernel function taken on a $p\times q$ submatrix indexed
by the rows $\mathbf{i}\in\mathcal{P}_{p}(\mathbb{N})$ and columns
$\mathbf{j}\in\mathcal{P}_{q}(\mathbb{N})$ can be written
$h(Y_{\mathbf{i},\mathbf{j}})\overset{a.s.}{=}h_{\phi}((\xi_{i})_{i\in\mathbf{i}};(\eta_{j})_{j\in\mathbf{j}};(\zeta_{ij})_{\begin{subarray}{c}i\in\mathbf{i}\\\
j\in\mathbf{j}\end{subarray}})$, where $h_{\phi}$ is some function depending
on $h$ and $\phi$.

Note that the AHK decomposition is not unique. In the rest of the paper, we
assume that for each dissociated RCE matrix $Y$, we have picked an AHK
representation, i.e. a suitable function $\phi$, and suitable i.i.d. random
variables $(\xi_{i})_{i\geq 1}$, $(\eta_{j})_{j\geq 1}$ and
$(\zeta_{ij})_{\begin{subarray}{c}i\geq 1\\\ j\geq 1\end{subarray}}$. In the
rest of the paper, we will write, abusively but without ambiguity,

$Y_{ij}=\phi(\xi_{i},\eta_{j},\zeta_{ij}).$

and

$h(Y_{\mathbf{i},\mathbf{j}})=h_{\phi}((\xi_{i})_{i\in\mathbf{i}};(\eta_{j})_{j\in\mathbf{j}};(\zeta_{ij})_{\begin{subarray}{c}i\in\mathbf{i}\\\
j\in\mathbf{j}\end{subarray}}).$

For $\mathbf{i}^{\prime}\in\mathcal{P}(\mathbb{N})$ and
$\mathbf{j}^{\prime}\in\mathcal{P}(\mathbb{N})$, we define the
$\sigma$-algebra

$\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}:=\sigma((\xi_{i})_{i\in\mathbf{i}^{\prime}},(\eta_{j})_{j\in\mathbf{j}^{\prime}},(\zeta_{ij})_{\begin{subarray}{c}i\in\mathbf{i}^{\prime}\\\
j\in\mathbf{j}^{\prime}\end{subarray}}).$

Therefore, it follows from our notations that

$\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}]=\mathbb{E}[h_{\phi}((\xi_{i})_{i\in\mathbf{i}};(\eta_{j})_{j\in\mathbf{j}};(\zeta_{ij})_{\begin{subarray}{c}i\in\mathbf{i}\\\
j\in\mathbf{j}\end{subarray}})\mid(\xi_{i})_{i\in\mathbf{i}^{\prime}};(\eta_{j})_{j\in\mathbf{j}^{\prime}};(\zeta_{ij})_{\begin{subarray}{c}i\in\mathbf{i}^{\prime}\\\
j\in\mathbf{j}^{\prime}\end{subarray}}].$

Now, notice that for fixed $\mathbf{i}^{\prime}$ and $\mathbf{j}^{\prime}$,
the quantity
$\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}]$
only depends on the elements shared by $\mathbf{i}$ and $\mathbf{i}^{\prime}$
and the elements shared by $\mathbf{j}$ and $\mathbf{j}^{\prime}$, and not on
the other elements of $\mathbf{i}$, $\mathbf{i}^{\prime}$, $\mathbf{j}$ and
$\mathbf{j}^{\prime}$. Suppose
$r=\mathrm{Card}(\mathbf{i}\cap\mathbf{i}^{\prime})$ and
$c=\mathrm{Card}(\mathbf{j}\cap\mathbf{j}^{\prime})$. Without loss of
generality, we can assume that
$\mathbf{i}^{\prime}\in\mathcal{P}_{r}(\mathbf{i})$ and
$\mathbf{j}^{\prime}\in\mathcal{P}_{c}(\mathbf{j})$ so
$\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}]$
only depends on the $r$ elements of $\mathbf{i}^{\prime}$ and the $c$ elements
of $\mathbf{j}^{\prime}$. Therefore, we can define the quantities
$\psi^{r,c}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})$ such that

$\psi^{r,c}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}):=\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}],$

where the choice of $\mathbf{i}$ and $\mathbf{j}$ does not matter as long as
$\mathbf{i}^{\prime}\subset\mathbf{i}$ and
$\mathbf{j}^{\prime}\subset\mathbf{j}$. Note that
$\psi^{r,c}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})$ is simply a
notation and not a function of $Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}$.
If $\mathbf{i}^{\prime}=\emptyset$ or $\mathbf{j}^{\prime}=\emptyset$, we will
still use this notation, for example

$\psi^{r,c}h(Y_{\mathbf{i}^{\prime},\emptyset})=\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\mathbf{i}^{\prime},\emptyset}],$

despite $Y_{\mathbf{i}^{\prime},\emptyset}$ being undefined.

#### Hoeffding-type decomposition of the kernel

In the following, for elements of $\mathbb{N}^{2}$,
$(x,y)\leq(x^{\prime},y^{\prime})$ means that both $x\leq x^{\prime}$ and
$y\leq y^{\prime}$; $(x,y)<(x^{\prime},y^{\prime})$ means that, in addition,
$(x,y)\neq(x^{\prime},y^{\prime})$.

For all $\mathbf{i}\in\mathcal{P}_{r}(\mathbb{N})$ and
$\mathbf{j}\in\mathcal{P}_{c}(\mathbb{N})$, we define by recursion the
following quantity

$p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})=\psi^{r,c}h(Y_{\mathbf{i},\mathbf{j}})-\sum_{(0,0)\leq(r^{\prime},c^{\prime})<(r,c)}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i})\\\
\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\end{subarray}}p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}).$
(4)

Since $\psi^{p,q}h(Y_{\mathbf{i},\mathbf{j}})=h(Y_{\mathbf{i},\mathbf{j}})$
for $\mathbf{i}\in\mathcal{P}_{p}(\mathbb{N})$ and
$\mathbf{j}\in\mathcal{P}_{q}(\mathbb{N})$, (4) yields the decomposition of
the kernel function $h$

$h(Y_{\mathbf{i},\mathbf{j}})=\sum_{(0,0)\leq(r,c)\leq(p,q)}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r}(\mathbf{i})\\\
\mathbf{j}^{\prime}\in\mathcal{P}_{c}(\mathbf{j})\end{subarray}}p^{r,c}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}).$
(5)

######  Remark 1.

From this formula, we see that $h(Y_{\mathbf{i},\mathbf{j}})$ is a linear
combination of the
$(p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}))_{\begin{subarray}{c}0\leq
r^{\prime}\leq p,0\leq c^{\prime}\leq q\\\
\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i}),\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\end{subarray}}$,
therefore, this is a linear combination of the
$(\psi^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}))_{\begin{subarray}{c}0\leq
r^{\prime}\leq r,0\leq c^{\prime}\leq c\\\
\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i}),\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\end{subarray}}$.

Now, we show that for $\mathbf{i}^{\prime}\subset\mathbf{i}$ and
$\mathbf{j}^{\prime}\subset\mathbf{j}$,
$p^{r,c}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})$ is actually the
projection of $h(Y_{\mathbf{i},\mathbf{j}})$ on the subspace generated by
$L_{2}$ functions of all the AHK variables of
$\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}$, orthogonally to the
sub-spaces generated by $L_{2}$ functions of all the variables of
$\mathcal{A}_{\mathbf{i}^{\prime\prime},\mathbf{j}^{\prime\prime}}$, for
$\mathbf{i}^{\prime\prime}\subset\mathbf{i}^{\prime}$ and
$\mathbf{j}^{\prime\prime}\subset\mathbf{j}^{\prime}$. This system of
projection is analogous to the Hoeffding decomposition for the kernel
functions of usual $U$-statistics on i.i.d. data. We prove that the following
orthogonality properties hold.

######  Proposition 2.1.

Let $h_{1}$ and $h_{2}$ two kernel functions of respective size $p_{1}\times
q_{1}$ and $p_{2}\times q_{2}$.

  1. 1.

Let $(0,0)\leq(r_{1},c_{1})\leq(p_{1},q_{1})$ and
$(0,0)\leq(r_{2},c_{2})\leq(p_{2},q_{2})$ such that
$(r_{1},c_{1})\neq(r_{2},c_{2})$. Let
$(\mathbf{i}_{1},\mathbf{j}_{1})\in\mathcal{P}_{r_{1}}(\llbracket
m\rrbracket)\times\mathcal{P}_{c_{1}}(\llbracket n\rrbracket)$ and
$(\mathbf{i}_{2},\mathbf{j}_{2})\in\mathcal{P}_{r_{2}}(\llbracket
m\rrbracket)\times\mathcal{P}_{c_{2}}(\llbracket n\rrbracket)$, then

$\mathrm{Cov}(p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}}))=0.$

  2. 2.

Let $(r,c)$ such that $(0,0)\leq(r,c)\leq(p_{1},q_{1})$ and
$(0,0)\leq(r,c)\leq(p_{2},q_{2})$. Let $(\mathbf{i}_{1},\mathbf{j}_{1})$ and
$(\mathbf{i}_{2},\mathbf{j}_{2})$ two elements of $\mathcal{P}_{r}(\llbracket
m\rrbracket)\times\mathcal{P}_{c}(\llbracket n\rrbracket)$. If
$(\mathbf{i}_{1},\mathbf{j}_{1})\neq(\mathbf{i}_{2},\mathbf{j}_{2})$, then

$\mathrm{Cov}(p^{r,c}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r,c}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}}))=0.$

This proposition relies on the fact that the projections are "conditionnally
centered". This property is given by the following lemma, the proof of which
is provided in Appendix B.

######  Lemma 2.2.

Let $h$ be a kernel function of size $p\times q$. Let
$(\mathbf{i},\mathbf{j})\in\mathcal{P}_{r}(\mathbb{N})\times\mathcal{P}_{c}(\mathbb{N})$,
where $(0,0)<(r,c)\leq(p,q)$. For all
$\underline{\mathbf{i}}\subset\mathbf{i}$ and
$\underline{\mathbf{j}}\subset\mathbf{j}$, we have

$\mathbb{E}[p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]=0.$

###### Proof of Proposition 2.1.

The two properties are proven similarly and derive from the fact that
$(\mathbf{i}_{1},\mathbf{j}_{1})\neq(\mathbf{i}_{2},\mathbf{j}_{2})$. This is
true for both properties.

Consider any (possibly equal) $(r_{1},c_{1})$ and $(r_{2},c_{2})$ and
associated
$(\mathbf{i}_{1},\mathbf{j}_{1})\neq(\mathbf{i}_{2},\mathbf{j}_{2})$. Then
$\mathbf{i}_{1}\neq\mathbf{i}_{2}$ or $\mathbf{j}_{1}\neq\mathbf{j}_{2}$.
Without loss of generality, assume that $\mathbf{i}_{1}\neq\mathbf{i}_{2}$ so
there is an element $i_{2}\in\mathbf{i}_{2}$ which is not included in
$\mathbf{i}_{1}$. Then

$\begin{split}\mathbb{E}[p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})]&=\mathbb{E}[\mathbb{E}[p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})\mid\mathcal{A}_{\\{i_{2}\\},\emptyset}]]\\\
&=\mathbb{E}[p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})\mathbb{E}[p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})\mid\mathcal{A}_{\\{i_{2}\\},\emptyset}]]\\\
&=0,\end{split}$

since
$\mathbb{E}[p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})\mid\mathcal{A}_{\\{i_{2}\\},\emptyset}]=0$
from Lemma 2.2.

We have proven both properties, since the case
$(r_{1},c_{1})\neq(r_{2},c_{2})$ corresponds to the first property and the
case $(r_{1},c_{1})=(r_{2},c_{2})$ corresponds to the second property. ∎

#### Decomposition of $U$-statistics

Using the Hoeffding-type decomposition of kernel functions (5) the
$U$-statistic (2) can be reformulated as:

$\begin{split}U^{h}_{m,n}(Y)&=\binom{m}{p}^{-1}\binom{n}{q}^{-1}\sum_{\begin{subarray}{c}1\leq
i_{1}<...<i_{p}\leq m\\\ 1\leq j_{1}<...<j_{q}\leq
n\end{subarray}}\sum_{(0,0)\leq(r,c)\leq(p,q)}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{r}(\\{i_{1},...,i_{p}\\})\\\
\mathbf{j}\in\mathcal{P}_{c}(\\{j_{1},...,j_{q}\\})\end{subarray}}p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\\\
&=\binom{m}{p}^{-1}\binom{n}{q}^{-1}\sum_{(0,0)\leq(r,c)\leq(p,q)}\binom{m-r}{p-r}\binom{n-c}{q-c}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{r}(\llbracket
m\rrbracket)\\\ \mathbf{j}\in\mathcal{P}_{c}(\llbracket
n\rrbracket)\end{subarray}}p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\\\
&=\sum_{(0,0)\leq(r,c)\leq(p,q)}\binom{p}{r}\binom{q}{c}\binom{m}{r}^{-1}\binom{n}{c}^{-1}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{r}(\llbracket
m\rrbracket)\\\ \mathbf{j}\in\mathcal{P}_{c}(\llbracket
n\rrbracket)\end{subarray}}p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\\\
&=\sum_{(0,0)\leq(r,c)\leq(p,q)}\binom{p}{r}\binom{q}{c}P^{r,c}_{m,n}h(Y)\end{split}$

where for all $0\leq r\leq p$ and $0\leq c\leq q$,

$P^{r,c}_{m,n}h(Y)=\binom{m}{r}^{-1}\binom{n}{c}^{-1}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{r}(\llbracket
m\rrbracket)\\\ \mathbf{j}\in\mathcal{P}_{c}(\llbracket
n\rrbracket)\end{subarray}}p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})$

is the $U$-statistic of kernel function $p^{r,c}h$ taken on the first $m\times
n$ rows and columns of the matrix $Y$. A consequence of the orthogonality of
the projections of $h$ is the orthogonality of these $U$-statistics, as stated
in the following corollary which is proven in Appendix B.

######  Corollary 2.3.

Let $h_{1}$ and $h_{2}$ two kernel functions of respective sizes $p_{1}\times
q_{1}$ and $p_{2}\times q_{2}$. Let $(0,0)\leq(r_{1},c_{1})<(p_{1},q_{1})$ and
$(0,0)\leq(r_{2},c_{2})<(p_{2},q_{2})$.

  1. 1.

If $(r_{1},c_{1})\neq(r_{2},c_{2})$, then

$\mathrm{Cov}(P^{r_{1},c_{1}}_{m,n}h_{1}(Y),P^{r_{2},c_{2}}_{m,n}h_{2}(Y))=0.$

  2. 2.

If $(r_{1},c_{1})=(r_{2},c_{2})=(r,c)$, then

$\mathrm{Cov}(P^{r,c}_{m,n}h_{1}(Y),P^{r,c}_{m,n}h_{2}(Y))=\binom{m}{r}^{-1}\binom{n}{c}^{-1}\mathrm{Cov}(p^{r,c}h_{1}(Y_{\llbracket
r\rrbracket,\llbracket c\rrbracket}),p^{r,c}h_{2}(Y_{\llbracket
r\rrbracket,\llbracket c\rrbracket})).$

The orthogonality between the $P^{r_{1},c_{1}}_{m,n}h_{1}(Y)$ and
$P^{r_{2},c_{2}}_{m,n}h_{2}(Y)$ allows to decompose the covariance of two
$U$-statistics into a few covariance terms.

######  Corollary 2.4.

$\begin{split}\mathrm{Cov}(U^{h_{1}}_{m,n}(Y),U^{h_{2}}_{m,n}(Y))&=\sum_{(0,0)<(r,c)\leq(p,q)}\binom{p}{r}^{2}\binom{q}{c}^{2}\mathrm{Cov}(P^{r,c}_{m,n}h_{1}(Y),P^{r,c}_{m,n}h_{2}(Y))\\\
&=\sum_{(0,0)<(r,c)\leq(p,q)}\binom{p}{r}^{2}\binom{q}{c}^{2}\binom{m}{r}^{-1}\binom{n}{c}^{-1}\mathrm{Cov}(p^{r,c}h_{1}(Y_{\llbracket
r\rrbracket,\llbracket c\rrbracket}),p^{r,c}h_{2}(Y_{\llbracket
r\rrbracket,\llbracket c\rrbracket})).\end{split}$

Corollary 2.4 is a direct consequence of Corollary 2.3. This result will be
helpful when considering the asymptotic properties of the $U$-statistics.
Indeed, one can see that each covariance term is associated with a binomial
coefficient depending on different orders of $m$ and $n$.

##  3 Asymptotic normality of $U$-statistics

For the following sections, we will use simplified notations, summarizing the
couple $(m_{N},n_{N})$ into $N$. We recall that
$U^{h}_{N}(Y)=U^{h}_{m_{N},n_{N}}(Y)$. We also denote
$P^{r,c}_{N}h(Y):=P^{r,c}_{m_{N},n_{N}}h(Y)$. When this is unambiguous, we
will omit to mention $Y$, so we will simply write $U^{h}_{N}$, $P^{r,c}_{N}h$,
$p^{r,c}_{(\mathbf{i},\mathbf{j})}h$ and
$\psi^{r,c}_{(\mathbf{i},\mathbf{j})}h$ instead of $U^{h}_{N}(Y)$,
$P^{r,c}_{N}h(Y)$, $p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})$ and
$\psi^{r,c}h(Y_{\mathbf{i},\mathbf{j}})$.

The following theorem is a Central Limit Theorem for $U^{h}_{N}$. Denote
$v^{r,c}_{h}:=\mathbb{V}[\psi^{r,c}_{(\llbracket r\rrbracket,\llbracket
c\rrbracket)}h]$.

######  Theorem 3.1.

Let $Y$ be a dissociated RCE matrix. Let $h$ be a $p\times q$ kernel function
such that $\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}]<\infty$. Let
$(m_{N},n_{N})_{N\geq 1}$ be a sequence of dimensions for the $U$-statistics,
such that $\frac{m_{N}}{N}\xrightarrow[N\rightarrow\infty]{}\rho$ and
$\frac{n_{N}}{N}\xrightarrow[N\rightarrow\infty]{}1-\rho$, where
$\rho\in]0,1[$. Let $(U^{h}_{N})_{N\geq 1}$ be the sequence of $U$-statistics
associated with $h$ defined by $U^{h}_{N}:=U^{h}_{m_{N},n_{N}}$. Set
$U^{h}_{\infty}=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]$ and

$\begin{split}V^{h}=&~{}\frac{p^{2}}{\rho}v^{1,0}_{h}+\frac{q^{2}}{1-\rho}v^{0,1}_{h}.\end{split}$

If $V^{h}>0$, then

$\sqrt{N}(U^{h}_{N}-U^{h}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,V^{h}).$

This theorem comes from the decomposition of
$\sqrt{N}(U^{h}_{N}-U^{h}_{\infty})$ into three different terms, the limits of
which are given by the following lemmas (proofs in Appendix C).

######  Lemma 3.2.

If $v^{1,0}_{h}>0$, then we have

$\frac{1}{\sqrt{m}}\sum_{i=1}^{m}p^{1,0}_{(\\{i\\},\emptyset)}h\xrightarrow[m\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,v^{1,0}_{h})$

and if $v^{0,1}_{h}>0$, then we have

$\frac{1}{\sqrt{n}}\sum_{j=1}^{n}p^{0,1}_{(\\{j\\},\emptyset)}h\xrightarrow[n\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,v^{0,1}_{h})$

######  Lemma 3.3.

Let $A_{N}:=\sqrt{N}\sum_{\begin{subarray}{c}(0,0)<(r,c)\leq(p,q)\\\
(r,c)\neq(1,0)\neq(0,1)\end{subarray}}\binom{p}{r}\binom{q}{c}P^{r,c}_{N}h$.
Then $A_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}0$.

###### Proof of Theorem 3.1.

We have

$\begin{split}U^{h}_{N}&=\sum_{(0,0)\leq(r,c)\leq(p,q)}\binom{p}{r}\binom{q}{c}P^{r,c}_{N}h\\\
&=P^{0,0}_{N}h+pP^{1,0}_{N}h+qP^{0,1}_{N}h+\sum_{\begin{subarray}{c}(0,0)<(r,c)\leq(p,q)\\\
(r,c)\neq(1,0)\neq(0,1)\end{subarray}}\binom{p}{r}\binom{q}{c}P^{r,c}_{N}h.\end{split}$

First, we see that $P^{0,0}_{N}h=U^{h}_{\infty}$. Next, $A_{N}$ being defined
in Lemma 3.3, we have

$\sqrt{N}(U^{h}_{N}-U^{h}_{\infty})=\frac{\sqrt{N}p}{m_{N}}\sum_{i=1}^{m_{N}}p^{1,0}_{(\\{i\\},\emptyset)}h+\frac{\sqrt{N}q}{n_{N}}\sum_{j=1}^{n_{N}}p^{0,1}_{(\emptyset,\\{j\\})}h+A_{N}.$

From Lemma 3.3, $A_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}0$. So by
Slutsky’s theorem, $\sqrt{N}(U^{h}_{N}-U^{h}_{\infty})$ has the same limiting
distribution as the two main terms of this decomposition. From Lemma 3.2, this
is the sum of two centered Gaussians of respective variance
$\frac{p^{2}}{\rho}v^{1,0}_{h}$ and $\frac{q^{2}}{1-\rho}v^{0,1}_{h}$.
Furthermore, $\sum_{i=1}^{m}p^{1,0}_{(\\{i\\},\emptyset)}h$ and
$\sum_{j=1}^{n}p^{0,1}_{(\\{j\\},\emptyset)}h$ are independent, so the two
Gaussians are independent, which concludes the proof. ∎

######  Remark 2.

The expression of $V^{h}$ could have been predicted with Corollary 2.4.
Indeed, this corollary implies that
$\mathbb{V}[U^{h}_{N}]=\sum_{(0,0)<(r,c)\leq(p,q)}\binom{p}{r}^{2}\binom{q}{c}^{2}\binom{m_{N}}{r}^{-1}\binom{n_{N}}{c}^{-1}\mathbb{V}[p^{r,c}_{(\llbracket
r\rrbracket,\llbracket c\rrbracket)}h]$, so
$\lim_{N\rightarrow\infty}N\mathbb{V}[U^{h}_{N}]=\lim_{N\rightarrow\infty}\left(\frac{p^{2}N}{m_{N}}\mathbb{V}[p^{1,0}_{(\\{1\\},\emptyset)}h]+\frac{q^{2}N}{n_{N}}\mathbb{V}[p^{0,1}_{(\emptyset,\\{1\\})}h]\right)=\frac{p^{2}}{c}\mathbb{V}[\psi^{1,0}_{(\\{1\\},\emptyset)}h]+\frac{q^{2}}{1-c}\mathbb{V}[\psi^{0,1}_{(\emptyset,\\{1\\})}h]$.

Now, we show that the limiting distribution of a vector of $U$-statistics is a
multivariate normal distribution under the condition that all the kernel
functions are linearly independent. However, if the kernels functions are of
different sizes, the notion of linear independence is unclear. We need to
define the concept of kernel extension to enunciate the corollary.

######  Definition 3.4.

Let $h$ be a kernel function of size $p\times q$. Let $p^{\prime}\geq p$ and
$q^{\prime}\geq q$. We define the extension of $h$ to the size
$p^{\prime}\times q^{\prime}$ by $\tilde{h}$ such that for all
$\mathbf{i}^{\prime}\in\mathcal{P}_{p^{\prime}}(\mathbb{N})$ and
$\mathbf{j}^{\prime}\in\mathcal{P}_{q^{\prime}}(\mathbb{N})$,

$\tilde{h}(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})=\left[\binom{p^{\prime}}{p}\binom{q^{\prime}}{q}\right]^{-1}\sum_{\begin{subarray}{c}\mathbf{i}\subset\mathcal{P}_{p}(\mathbf{i}^{\prime})\\\
\mathbf{j}\subset\mathcal{P}_{q}(\mathbf{j}^{\prime})\end{subarray}}h(Y_{\mathbf{i},\mathbf{j}}).$

The extension of a kernel actually shares some properties with its kernel, as
shown in Lemma C.1. With this definition, we can define the linear
independence of kernel functions needed for the following corollary as the
linear independence of their kernel extensions. Denote
$c^{r,c}_{h_{k},h_{\ell}}:=\mathrm{Cov}\left(\psi^{r,c}_{(\llbracket
r\rrbracket,\llbracket c\rrbracket)}h_{k},\psi^{r,c}_{(\llbracket
r\rrbracket,\llbracket c\rrbracket)}h_{\ell}\right)$. The proof of this
corollary can be found in Appendix C.

######  Corollary 3.5.

Let $Y$ be a dissociated RCE matrix. Let $(h_{1},h_{2},...,h_{D})$ be a vector
of kernel functions of respective sizes $p_{1}\times q_{1},p_{2}\times
q_{2},...,p_{D}\times q_{D}$ such that

  1. 1.

Theorem 3.1 applies for each kernel function, i.e.
$\mathbb{E}[h_{k}(Y_{(1,2;1,2)})^{2}]<\infty$ and $U^{h_{k}}_{\infty}$ and
$V^{h_{k}}$ are as defined in Theorem 3.1 for each kernel $h_{k}$, $1\leq
k\leq D$,

  2. 2.

for $t\in\mathbb{R}^{D}$, $\sum_{k=1}^{D}t_{k}\tilde{h}_{k}\equiv 0$ if and
only if $t=(0,...,0)$, where for $1\leq k\leq D$, $\tilde{h}_{k}$ is the
extension of $h_{k}$ to size $\max_{k}(p_{k})\times\max_{k}(q_{k})$.

Then

$\sqrt{N}\left(\begin{pmatrix}U^{h_{1}}_{N}\\\ U^{h_{2}}_{N}\\\ ...\\\
U^{h_{D}}_{N}\end{pmatrix}-\begin{pmatrix}U^{h_{1}}_{\infty}\\\
U^{h_{2}}_{\infty}\\\ ...\\\
U^{h_{D}}_{\infty}\end{pmatrix}\right)\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,\Sigma^{h_{1},...,h_{D}}),$

with

$\Sigma^{h_{1},...,h_{D}}=\left(C^{h_{k},h_{\ell}}\right)_{1\leq k,\ell\leq
D},$

where
$C^{h_{k},h_{\ell}}=\frac{p^{2}}{\rho}c^{1,0}_{h_{k},h_{\ell}}+\frac{q^{2}}{1-\rho}c^{0,1}_{h_{k},h_{\ell}}$
for all $1\leq k,\ell\leq D$ (and $C^{h_{k},h_{k}}=V^{h_{k}}$).

Although, through the first condition, the theorem requires the kernel
functions of $(h_{1},h_{2},...,h_{D})$ to be linearly independent, the
corresponding $U$-statistics are not independent random variables, even
asymptotically, because $\Sigma^{h_{1},...,h_{D}}$ is not a diagonal matrix.
One consequence of this corollary is that Theorem 3.1 can be extended to
differentiable functions of $U$-statistics.

######  Corollary 3.6.

Let $h_{1},...,h_{D}$ be $D$ kernel functions such that Corollary 3.5 applies
and $\Sigma^{h_{1},...,h_{D}}=\left(C^{h_{k},h_{\ell}}\right)_{1\leq
k,\ell\leq D}$. Denote $\theta=(U^{h_{1}}_{\infty},...,U^{h_{D}}_{\infty})$.
Let $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a differentiable function at
$\theta$. Denote $\nabla g$ the gradient of g and $V^{\delta}:=\nabla
g(\theta)^{T}\Sigma^{h_{1},...,h_{D}}\nabla g(\theta)$. If $V^{\delta}>0$,
then

$\sqrt{N}(g(U^{h_{1}}_{N},...,U^{h_{D}}_{N})-g(\theta))\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,V^{\delta}).$

###### Proof.

The first-order Taylor expansion of $g$ at $\theta$ is written

$g(U^{h_{1}}_{N},...,U^{h_{D}}_{N})-g(\theta)=\nabla
g(\theta)^{T}\left((U^{h_{1}}_{N},...,U^{h_{D}}_{N})-\theta\right)+o_{P}\left(\lvert\lvert(U^{h_{1}}_{N},...,U^{h_{D}}_{N})-\theta\rvert\rvert\right).$

From Corollary 3.5,
$\sqrt{N}\left((U^{h_{1}}_{N},...,U^{h_{D}}_{N})-\theta\right)$ converges to a
multivariate normal distribution with asymptotic covariance matrix
$\Sigma^{h_{1},...,h_{D}}$, so the delta method (see Theorem 3.1 of [50]) can
be applied to prove this proposition. ∎

##  4 Estimation of the asymptotic variance of a non-degenerated $U$-statistic

Theorem 3.1 shows the asymptotic normality of RCE submatrix $U$-statistics. In
order to perform statistical inference using these $U$-statistics, one needs
to estimate their variances. We see that

$\mathbb{V}[U^{h}_{N}]=\left[\binom{m}{p}\binom{n}{q}\right]^{-2}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{p}(\llbracket
m\rrbracket)\\\ \mathbf{j}\in\mathcal{P}_{q}(\llbracket
n\rrbracket)\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{p}(\llbracket
m\rrbracket)\\\ \mathbf{j}^{\prime}\in\mathcal{P}_{q}(\llbracket
n\rrbracket)\end{subarray}}\mathrm{Cov}\left(h(Y_{\mathbf{i},\mathbf{j}}),h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\right).$

By exchangeability, the covariance between the kernels
$h(Y_{\mathbf{i},\mathbf{j}})$ and
$h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})$ only depends on the number of
row and column indices they share. Denote
$\gamma^{r,c}_{h}:=\mathrm{Cov}\left(h(Y_{\mathbf{i},\mathbf{j}}),h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\right)$
where $Y_{\mathbf{i},\mathbf{j}}$ and
$Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}$ share $r$ row indices and $c$
column indices: $\mathrm{Card}(\mathbf{i}\cap\mathbf{i}^{\prime})=r$ and
$\mathrm{Card}(\mathbf{j}\cap\mathbf{j}^{\prime})=c$. We get

$\mathbb{V}[U^{h}_{N}]=\left[\binom{m}{p}\binom{n}{q}\right]^{-1}\sum_{(0,0)\leq(r,c)\leq(p,q)}\binom{p}{r}\binom{q}{c}\binom{m-p}{p-r}\binom{n-q}{q-c}\gamma^{r,c}_{h}.$

Each $\gamma^{r,c}_{h}$, $1\leq r\leq p,1\leq c\leq q$, can be estimated using
empirical covariance estimators, between kernel terms that share $r$ rows and
$c$ columns and in particular. This leads to an unbiased estimator of
$\mathbb{V}[U^{h}_{N}]$ similar to the one discussed by [40] for
$U$-statistics of one-dimensional i.i.d. arrays. However, the estimation of
these covariances is computationally intensive and the estimators can take
negative values, which can lead to a negative variance estimation.

One approach is to estimate the asymptotic variance $V^{h}$. The asymptotic
variance formula given by Theorem 3.1 is
$V^{h}=\frac{p^{2}}{\rho}v^{1,0}_{h}+\frac{q^{2}}{1-\rho}v^{0,1}_{h}$. Remark
that $v^{1,0}_{h}=\gamma^{1,0}_{h}$ and $v^{0,1}_{h}=\gamma^{0,1}_{h}$. It is
often tedious to analytically calculate $V^{h}$, especially as it depends on
$h$ and the distribution of $Y$, see for example Section 3 of [29].

We present here a kernel-free and model-free estimator of $V^{h}$, taking
advantage of the Hoeffding decomposition. Indeed, we will first define
estimators for the conditional expectations
$(\psi^{1,0}_{(\\{i\\},\emptyset)}h)_{1\leq i\leq m}$ and
$(\psi^{0,1}_{(\emptyset,\\{j\\})}h)_{1\leq j\leq n}$. Then, using the fact
that $v^{1,0}_{h}=\mathbb{V}[\psi^{1,0}_{(\\{1\\},\emptyset)}h]$ and
$v^{0,1}_{h}=\mathbb{V}[\psi^{0,1}_{(\emptyset,\\{1\\})}h]$, we can derive a
positive estimator for $V^{h}$.

$V^{h}$ is the asymptotic variance of a single $U$-statistic. Afterwards, we
explain how to generalize this estimator to estimate the asymptotic variance
of a function of $U$-statistics to apply Corollary 3.6.

###  4.1 Some useful notations and results

First, we introduce further notations and a helpful lemma for this section.

For some $N>0$, the size of the overall RCE matrix is $m_{N}\times n_{N}$. $i$
being a row index means that $1\leq i\leq m_{N}$ and $j$ being a column index
means that $1\leq j\leq n_{N}$.

We denote

$X_{\mathbf{i},\mathbf{j}}:=h(Y_{\mathbf{i},\mathbf{j}}).$

For $N$ such that $m_{N}\geq p$ and $n_{N}\geq q$, we further denote

$\mathcal{S}^{p,q}_{N}:=\left\\{(\mathbf{i},\mathbf{j}):\mathbf{i}\in\mathcal{P}_{p}(\llbracket
m_{N}\rrbracket),\mathbf{j}\in\mathcal{P}_{q}(\llbracket
n_{N}\rrbracket)\right\\},$

so that the set of the kernels taken on all the $p\times q$ submatrices can be
written as
$(X_{\mathbf{i},\mathbf{j}})_{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N}}$.

Let $\mathbf{\underline{i}}$ be a set of row indices of size $\underline{p}$
such that $0\leq\underline{p}\leq p$ and $\mathbf{\underline{j}}$ a set of
column indices of size $\underline{q}$ such that $0\leq\underline{q}\leq q$.
The subset of $\mathcal{S}^{p,q}_{N}$ where $\mathbf{\underline{i}}$ is
included in the row indices and $\mathbf{\underline{j}}$ is included in the
column indices is denoted

$\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}},\mathbf{\underline{j}})}:=\left\\{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N}:\mathbf{\underline{i}}\subset\mathbf{i},\mathbf{\underline{j}}\subset\mathbf{j}\right\\}.$

For example, the set of the kernels taken on all the $p\times q$ submatrices
containing the columns $1$ and $2$ can be written
$(X_{\mathbf{i},\mathbf{j}})_{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N,(\emptyset,\\{1,2\\})}}$.

It is obvious that

$\mathrm{Card}(\mathcal{S}^{p,q}_{N})=\binom{m_{N}}{p}\binom{n_{N}}{q}\quad\quad\mbox{and}\quad\quad\mathrm{Card}(\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}},\mathbf{\underline{j}})})=\binom{m_{N}-\underline{p}}{p-\underline{p}}\binom{n_{N}-\underline{q}}{q-\underline{q}}.$

Let $I_{K}=(\mathbf{\underline{i}}_{1},...,\mathbf{\underline{i}}_{K})$ be a
$K$-uplet of sets of row indices and
$J_{K}=(\mathbf{\underline{j}}_{1},...,\mathbf{\underline{j}}_{K})$ a
$K$-uplet of subsets of column indices. Denote

$\mathcal{T}^{p,q}_{N,(I_{K},J_{K})}:=\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}}_{1},\mathbf{\underline{j}}_{1})}\times\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}}_{2},\mathbf{\underline{j}}_{2})}\times...\times\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}}_{K},\mathbf{\underline{j}}_{K})}.$

In the rest of the paper, we will often see averages of the type

$T^{p,q}_{N}(I_{K},J_{K}):=\frac{1}{\prod_{k=1}^{K}\mathrm{Card}(\mathcal{T}^{p,q}_{N,(I_{K},J_{K})})}\sum_{\mathcal{T}^{p,q}_{N,(I_{K},J_{K})}}\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}X_{\mathbf{i}_{2},\mathbf{j}_{2}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}].$
(6)

As a remark, $T^{p,q}_{N}((\emptyset),(\emptyset))=U^{h}_{N}$.

By exchangeability, the quantities
$\mathbb{E}[X_{(\mathbf{i}_{1},\mathbf{j}_{1})}X_{(\mathbf{i}_{2},\mathbf{j}_{2})}...X_{(\mathbf{i}_{K},\mathbf{j}_{K})}]$
do not depend on the row indices that do not belong to any pairwise
intersection of the $(\mathbf{i}_{k})_{1\leq k\leq K}$. The same holds for
column indices that do not belong to any pairwise intersection of the
$(\mathbf{j}_{k})_{1\leq k\leq K}$. Therefore, assuming
$m_{N}\geq\mathrm{Card}(\cup_{k=1}^{K}\mathbf{\underline{i}}_{k})$ and
$n_{N}\geq\mathrm{Card}(\cup_{k=1}^{K}\mathbf{\underline{j}}_{k})$, we can
define

$\alpha(I_{K},J_{K}):=\mathbb{E}[X_{\mathbf{\bar{\underline{i}}}_{1},\mathbf{\bar{\underline{j}}}_{1}}X_{\mathbf{\bar{\underline{i}}}_{2},\mathbf{\bar{\underline{j}}}_{2}}...X_{\mathbf{\bar{\underline{i}}}_{K},\mathbf{\bar{\underline{j}}}_{K}}]$
(7)

where for $1\leq k\leq K$, the $p$-uplet $\mathbf{\bar{\underline{i}}}_{k}$
only consists of elements of $\mathbf{\underline{i}}_{k}$ and elements that
are not in any of the other $\mathbf{\underline{i}}_{k^{\prime}}$, i.e. the
$\mathbf{\bar{\underline{i}}}_{k}$ are of the form
$\mathbf{\bar{\underline{i}}}_{k}=\mathbf{\underline{i}}_{k}\cup\tilde{\mathbf{i}}_{k}$
where $\cap_{k=1}^{K}\tilde{\mathbf{i}}_{k}=\emptyset$ and
$(\cup_{k=1}^{K}\tilde{\mathbf{i}}_{k})\cap(\cup_{k=1}^{K}\underline{\mathbf{i}}_{k})=\emptyset$.

The following lemma will be helpful in later proofs as it provides the
asymptotic behaviour of $T^{p,q}_{N}(I_{K},J_{K})$. It shows that these
averages can be reduced to one dominant expectation term given by
$\alpha(I_{K},J_{K})$ and a remainder vanishing as $N$ grows.

######  Lemma 4.1.

Let $I_{K}=(\mathbf{\underline{i}}_{1},...,\mathbf{\underline{i}}_{K})$ and
$J_{K}=(\mathbf{\underline{j}}_{1},...,\mathbf{\underline{j}}_{K})$ be
$K$-uplets of respectively row and column indices. Let $\alpha(I_{K},J_{K})$
defined by (7). Let :

  * •

$\underline{p}_{k}:=\mathrm{Card}(\mathbf{\underline{i}}_{k})$ and
$\underline{q}_{k}:=\mathrm{Card}(\mathbf{\underline{j}}_{k})$, for $1\leq
k\leq K$,

  * •

$\underline{P}:=\sum_{k=1}^{K}\underline{p}_{k}$ and
$\underline{Q}:=\sum_{k=1}^{K}\underline{q}_{k}$,

  * •

$\bar{p}:=\mathrm{Card}(\cup_{k=1}^{K}\mathbf{\underline{i}}_{k})$ and
$\bar{q}:=\mathrm{Card}(\cup_{k=1}^{K}\mathbf{\underline{j}}_{k})$.

We have

$T^{p,q}_{N}(I_{K},J_{K})=\alpha(I_{K},J_{K})+O\left(m_{N}^{-1}+n_{N}^{-1}\right).$

###  4.2 Estimation of the conditional expectations

In this paragraph, for all $i\in\llbracket m_{N}\rrbracket$ and
$j\in\llbracket n_{N}\rrbracket$, we define estimators for
$\psi^{1,0}_{(\\{i\\},\emptyset)}h=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})\mid\xi_{i}]$
and
$\psi^{0,1}_{(\emptyset,\\{j\\})}h=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})\mid\eta_{j}]$,
where $\xi_{i}$ and $\eta_{j}$ have been defined in Section 2.

Let $\widehat{\mu}^{h,(i)}_{N}$ be the average of the kernel function applied
on the $p\times q$ submatrices containing the row $i$. Symmetrically, let
$\widehat{\nu}^{h,(j)}_{N}$ be the average of the kernel function applied on
the $p\times q$ submatrices containing the column $j$. This means

$\widehat{\mu}^{h,(i)}_{N}:=\binom{m_{N}-1}{p-1}^{-1}\binom{n_{N}}{q}^{-1}\sum_{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N,(\\{i\\},\emptyset)}}h(Y_{\mathbf{i},\mathbf{j}}),$
(8)

and

$\widehat{\nu}^{h,(j)}_{N}:=\binom{m_{N}}{p}^{-1}\binom{n_{N}-1}{q-1}^{-1}\sum_{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N,(\emptyset,\\{j\\})}}h(Y_{\mathbf{i},\mathbf{j}}).$
(9)

Now, we establish some properties for these estimators.

######  Proposition 4.2.

If $Y$ is a RCE matrix, then $\widehat{\mu}^{h,(i)}_{N}$ and
$\widehat{\nu}^{h,(i)}_{N}$ are both conditionally unbiased
$\psi^{1,0}_{(\\{i\\},\emptyset)}h$ and $\psi^{0,1}_{(\emptyset,\\{j\\})}h$,
i.e. we have for all $N\in\mathbb{N}$ :

  * •

$\mathbb{E}[\widehat{\mu}^{h,(i)}_{N}\mid\xi_{i}]=\psi^{1,0}_{(\\{i\\},\emptyset)}h$,

  * •

$\mathbb{E}[\widehat{\nu}^{h,(j)}_{N}\mid\eta_{j}]=\psi^{0,1}_{(\emptyset,\\{j\\})}h$.

###### Proof.

The first result can be found directly by the definition of
$\psi^{1,0}_{(\\{i\\},\emptyset)}h$, since

$\begin{split}\mathbb{E}[\widehat{\mu}^{h,(i)}_{N}\mid\xi_{i}]&=\binom{m_{N}-1}{p-1}^{-1}\binom{n_{N}}{q}^{-1}\sum_{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N,(\\{i\\},\emptyset)}}\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\xi_{i}]\\\
&=\binom{m_{N}-1}{p-1}^{-1}\binom{n_{N}}{q}^{-1}\sum_{(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N,(\\{i\\},\emptyset)}}\psi^{1,0}_{(\\{i\\},\emptyset)}h\\\
&=\psi^{1,0}_{(\\{i\\},\emptyset)}h.\end{split}$

The second result can be obtained analogously. ∎

######  Proposition 4.3.

If $Y$ is a RCE matrix, then :

  * •

$\widehat{\mu}^{h,(i)}_{N}\xrightarrow[N\rightarrow\infty]{a.s.,L_{1}}\psi^{1,0}_{(\\{i\\},\emptyset)}h$,

  * •

$\widehat{\nu}^{h,(j)}_{N}\xrightarrow[N\rightarrow\infty]{a.s.,L_{1}}\psi^{0,1}_{(\emptyset,\\{j\\})}h$.

As a consequence, $\widehat{\mu}^{h,(i)}_{N}$ and $\widehat{\nu}^{h,(j)}_{N}$
are consistent estimators for $\psi^{1,0}_{(\\{i\\},\emptyset)}h$ and
$\psi^{0,1}_{(\emptyset,\\{j\\})}h$.

###### Proof.

Let $N\in\mathbb{N}$ and
$\mathcal{F}_{N}(Y)=\sigma((\widehat{\mu}^{h,(i)}_{K}(Y))_{K\geq N})$. Let
$\Phi_{N}\in\mathbb{S}^{(i)}_{m_{N}}\times\mathbb{S}_{n_{N}}$ where
$\mathbb{S}^{(i)}_{m_{N}}$ is the group of permutations $\sigma_{i}$ of
$\llbracket m_{N}\rrbracket$ such that $\sigma_{i}(i)=i$. If
$\Phi_{N}=(\sigma_{i},\tau)$, denote
$\Phi_{N}Y=(Y_{\sigma_{i}(k)\tau(j)})_{\begin{subarray}{c}1\leq k\leq m_{N}\\\
1\leq j\leq n_{N}\end{subarray}}$.

First, we observe that
$\widehat{\mu}^{h,(i)}_{N}(Y)=\widehat{\mu}^{h,(i)}_{N}(\Phi_{N}Y)$, so
$\mathcal{F}_{N}(\Phi_{N}Y)=\mathcal{F}_{N}(Y)$. Therefore, by the
exchangeability of $Y$, we have

$Y\mid\mathcal{F}_{N}(Y)\overset{\mathcal{D}}{=}\Phi_{N}Y\mid\mathcal{F}_{N}(Y).$

This assertion is true for all $N\in\mathbb{N}$ and
$\Phi_{N}\in\mathbb{S}^{(i)}_{m_{N}}\times\mathbb{S}_{n_{N}}$. Now, note that
for all $(\mathbf{i},\mathbf{j})$ and
$(\mathbf{i}^{\prime},\mathbf{j}^{\prime})$ elements of
$\mathcal{S}^{p,q}_{N,(\\{i\\},\emptyset)}$, we can always find a permutation
$\Phi_{N}\in\mathbb{S}^{(i)}_{m_{N}}\times\mathbb{S}_{n_{N}}$ such that
$h(\Phi_{N}Y_{\mathbf{i},\mathbf{j}})=h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})$.
Thus, we have
$\mathbb{E}[h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{F}_{N}(Y)]=\mathbb{E}[h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{F}_{N}(Y)]$.
Hence, for any
$(\mathbf{i},\mathbf{j})\in\mathcal{S}^{p,q}_{N+1,(\\{i\\},\emptyset)}$, we
deduce that

$\begin{split}\mathbb{E}[\widehat{\mu}^{h,(i)}_{N}(Y)\mid\mathcal{F}_{N+1}(Y)]&=\mathbb{E}[h(Y_{(\mathbf{i},\mathbf{j})})\mid\mathcal{F}_{N+1}(Y)]\\\
&=\mathbb{E}[\widehat{\mu}^{h,(i)}_{N+1}(Y)\mid\mathcal{F}_{N+1}(Y)]\\\
&=\widehat{\mu}^{h,(i)}_{N+1}(Y).\end{split}$

Therefore, $\widehat{\mu}^{h,(i)}_{N}(Y)$ is a backward martingale with
respect to $\mathcal{F}_{N}(Y)$ (see Appendix A). By Theorem A.3, we have that
$\widehat{\mu}^{h,(i)}_{N}(Y)\xrightarrow[N\rightarrow\infty]{a.s.,L_{1}}\mathbb{E}[\widehat{\mu}^{h,(i)}_{1}(Y)\mid\mathcal{F}_{\infty}(Y)]$,
where $\mathcal{F}_{\infty}(Y)=\bigcap_{N=1}^{\infty}\mathcal{F}_{N}(Y)$.

Finally, $\mathcal{F}_{\infty}(Y)=\sigma(\xi_{i})$ so Proposition 4.2 implies
that
$\mathbb{E}[\widehat{\mu}^{h,(i)}_{1}(Y)\mid\mathcal{F}_{\infty}(Y)]=\psi^{1,0}_{(\\{i\\},\emptyset)}h(Y)$
and thus,
$\widehat{\mu}^{h,(i)}_{N}(Y)\xrightarrow[N\rightarrow\infty]{a.s.,L_{1}}\psi^{1,0}_{(\\{i\\},\emptyset)}h(Y)$.
∎

###  4.3 Estimation of $V^{h}$

Finally, since we have defined estimators for
$\psi^{1,0}_{(\\{i\\},\emptyset)}h$ and $\psi^{0,1}_{(\emptyset,\\{j\\})}h$,
natural estimators for
$v^{1,0}_{h}=\mathbb{V}[\psi^{1,0}_{(\\{1\\},\emptyset)}h]$ and
$v^{0,1}_{h}=\mathbb{V}[\psi^{0,1}_{(\emptyset,\\{j\\})}h]$ can be given by :

$\widehat{v}^{h;1,0}_{N}=\binom{m_{N}}{2}^{-1}\sum_{1\leq i_{1}<i_{2}\leq
m_{N}}\frac{(\widehat{\mu}^{h,(i_{1})}_{N}-\widehat{\mu}^{h,(i_{2})}_{N})^{2}}{2}$

$\widehat{v}^{h;0,1}_{N}=\binom{n_{N}}{2}^{-1}\sum_{1\leq j_{1}<j_{2}\leq
n_{N}}\frac{(\widehat{\nu}^{h,(j_{1})}_{N}-\widehat{\nu}^{h,(j_{2})}_{N})^{2}}{2}$

Then, an estimator for $V^{h}$ is

$\widehat{V}^{h}_{N}:=\frac{p^{2}}{\rho}\widehat{v}^{h;1,0}_{N}+\frac{q^{2}}{1-\rho}\widehat{v}^{h;0,1}_{N}.$

The following theorem shows that $\widehat{V}^{h}_{N}$ is a consistent
estimator for $V^{h}$.

######  Theorem 4.4.

We have
$\widehat{v}^{h;1,0}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}v^{1,0}_{h}$
and
$\widehat{v}^{h;0,1}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}v^{0,1}_{h}$.
As a consequence,
$\widehat{V}^{h}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}V^{h}$.

Before proving this theorem, we first need to highlight some properties of
$\widehat{v}^{h;1,0}_{N}$ and $\widehat{v}^{h;0,1}_{N}$, which are proven in
Appendix D.

######  Proposition 4.5.

We have $\mathbb{E}[\widehat{v}^{h;1,0}_{N}]=v^{1,0}_{h}+O\left(N^{-1}\right)$
and $\mathbb{E}[\widehat{v}^{h;0,1}_{N}]=v^{0,1}_{h}+O\left(N^{-1}\right)$. As
a consequence, $\widehat{v}^{h;1,0}_{N}$ and $\widehat{v}^{h;0,1}_{N}$ are
asymptotically unbiased estimators for $v^{1,0}_{h}$ and $v^{0,1}_{h}$.

######  Proposition 4.6.

We have $\mathbb{V}[\widehat{v}^{h;1,0}_{N}]=O\left(N^{-1}\right)$ and
$\mathbb{V}[\widehat{v}^{h;0,1}_{N}]=O\left(N^{-1}\right)$.

###### Proof of Theorem 4.4.

For some $\epsilon>0$, it follows from Proposition 4.5 that for large enough
values of $N$,
$\left\lvert\mathbb{E}[\widehat{v}^{h;1,0}_{N}]-v^{1,0}_{h}\right\rvert<\epsilon$.
The triangular inequality and Chebyshev’s inequality states that

$\mathbb{P}\left(\left\lvert\widehat{v}^{h;1,0}_{N}-v^{1,0}_{h}\right\rvert>\epsilon\right)\leq\mathbb{P}\left(\left\lvert\widehat{v}^{h;1,0}_{N}-\mathbb{E}[\widehat{v}^{h;1,0}_{N}]\right\rvert\geq\epsilon-\left\lvert\mathbb{E}[\widehat{v}^{h;1,0}_{N}]-v^{1,0}_{h}\right\rvert\right)\leq\frac{\mathbb{V}[\widehat{v}^{h;1,0}_{N}]}{\left(\epsilon-\left\lvert\mathbb{E}[\widehat{v}^{h;1,0}_{N}]-v^{1,0}_{h}\right\rvert\right)^{2}}.$

Applying Propositions 4.5 and 4.6 to the right-hand side of the inequality
ensures that
$\mathbb{P}\left(\left\lvert\widehat{v}^{h;1,0}_{N}-v^{1,0}_{h}\right\rvert>\epsilon\right)\xrightarrow[N\rightarrow\infty]{}0$
which concludes the proof. ∎

With Theorem 4.4, it is possible to use $\widehat{V}^{h}_{N}$ for statistical
inference tasks when plugged-in in place of $V^{h}$, a asymptotic normality
result similar to Theorem 3.1 holds.

######  Corollary 4.7.

If $V^{h}>0$, then

$\sqrt{\frac{N}{\widehat{V}^{h}_{N}}}\left(U^{h}_{N}-U^{h}_{\infty}\right)\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1).$

###### Proof.

First, Theorem 3.1 ensures that
$\sqrt{\frac{N}{\widehat{V}^{h}_{N}}}(U^{h}_{N}-U^{h}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1)$.
Second, it can be derived from Theorem 4.4 that
$\sqrt{\frac{V^{h}}{\widehat{V}^{h}_{N}}}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}1$.

Then, applying Slutsky’s theorem, we get
$\sqrt{\frac{N}{\widehat{V}^{h}_{N}}}(U^{h}_{N}-U^{h}_{\infty})=\sqrt{\frac{V^{h}}{\widehat{V}^{h}_{N}}}\times\sqrt{\frac{N}{V^{h}}}(U^{h}_{N}-U^{h}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1)$.
∎

###  4.4 Calculation of the estimator

In practice, the computation of the estimators
$(\widehat{\mu}^{h,(i)}_{N})_{1\leq i\leq m_{N}}$ and
$(\widehat{\nu}^{h,(j)}_{N})_{1\leq j\leq n_{N}}$ using their definition is
pretty tedious, as they are sums of $O(N^{p+q-1})$ terms. The computation cost
of $\widehat{V}^{h}_{N}$ is then $O(N^{p+q})$ which is of same order as the
computation of $U^{h}_{N}$ when naively applying the kernel function to all
the $p\times q$ submatrices and averaging.

However, one would actually try to avoid to compute $U^{h}_{N}$ in that way if
possible, as for simple kernels, it is possible to write $U^{h}_{N}$ in the
form of operation on matrices, which are more optimized, say $O(N^{a+b})$ with
$a\leq p$ and $b\leq q$.

In these cases, it might be time-efficient to compute
$\widehat{V}^{h}_{N}=\frac{p^{2}}{\rho}\widehat{v}^{h;1,0}_{N}+\frac{q^{2}}{1-\rho}\widehat{v}^{h;0,1}_{N}$
with the following alternative form.

######  Proposition 4.8.

An alternative form for $\widehat{v}^{h;1,0}_{N}$ and
$\widehat{v}^{h;0,1}_{N}$ is given by

$\widehat{v}^{h;1,0}_{N}=\frac{(m_{N}-p)^{2}}{p^{2}(m_{N}-1)}\sum_{i=1}^{m_{N}}\left(U^{h}_{N}-U_{N}^{h,(-i,\emptyset)}\right)^{2}$

and

$\widehat{v}^{h;0,1}_{N}=\frac{(n_{N}-q)^{2}}{q^{2}(n_{N}-1)}\sum_{j=1}^{n_{N}}\left(U^{h}_{N}-U_{N}^{h,(\emptyset,-j)}\right)^{2},$

where

$U_{N}^{h,(-i,\emptyset)}:=\left[\binom{m_{N}-1}{p}\binom{n_{N}}{q}\right]^{-1}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{p}(\llbracket
m_{N}\rrbracket\backslash\\{i\\})\\\ \mathbf{j}\in\mathcal{P}_{q}(\llbracket
n_{N}\rrbracket)\end{subarray}}h(Y_{\mathbf{i},\mathbf{j}})$

and

$U_{N}^{h,(\emptyset,-j)}:=\left[\binom{m_{N}}{p}\binom{n_{N}-1}{q}\right]^{-1}\sum_{\begin{subarray}{c}\mathbf{i}\in\mathcal{P}_{p}(\llbracket
m_{N}\rrbracket)\\\ \mathbf{j}\in\mathcal{P}_{q}(\llbracket
n_{N}\rrbracket\backslash\\{j\\})\end{subarray}}h(Y_{\mathbf{i},\mathbf{j}}).$

###### Proof.

The relations are obtained from the definition of $\widehat{v}^{h;1,0}_{N}$
and $\widehat{v}^{h;0,1}_{N}$ and noticing that

$\binom{m_{N}-1}{p-1}\binom{n_{N}}{q}\widehat{\mu}^{h,(i)}_{N}=\binom{m_{N}}{p}\binom{n_{N}}{q}U^{h}_{N}-\binom{m_{N}-1}{p}\binom{n_{N}}{q}U_{N}^{h,(-i,\emptyset)}$

and

$\binom{m_{N}}{p}\binom{n_{N}-1}{q-1}\widehat{\nu}^{h,(j)}_{N}=\binom{m_{N}}{p}\binom{n_{N}}{q}U^{h}_{N}-\binom{m_{N}}{p}\binom{n_{N}-1}{q}U_{N}^{h,(\emptyset,-j)}.$

∎

######  Remark.

This is one alternative method to compute $\widehat{V}^{h}_{N}$, but not
necessarily the optimal one. With this form, the computational cost of
$\widehat{V}^{h}_{N}$ is $O(N^{a+b+1})$. This can be outperformed by the naive
method in some specific cases, e.g. when $a=p$ and $b=q$. This method is
particularly fit when $a+b+1<p+q$.

######  Remark.

The alternative form of $\widehat{v}^{h;1,0}_{N}$ and
$\widehat{v}^{h;0,1}_{N}$ is reminiscent of the jackknife estimator for the
variance of $U$-statistics of one-dimensional arrays [2], but the two are well
distinct. In the case where $Y$ is a one-dimensional array, the $U$-statistic
associated to the kernel $h:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is

$U^{h}_{N}=\binom{N}{p}^{-1}\sum_{\mathbf{i}\in\mathcal{P}_{p}(\llbracket
N\rrbracket)}h(X_{\mathbf{i}}).$

The jackknife estimator of the asymptotic variance of this $U$-statistic is

$\widehat{V}^{h,J}_{N}=(N-1)\sum_{i=1}^{N}(U^{h}_{N}-U_{N}^{h,(-i)})^{2}$

where

$U_{N}^{h,(-i)}=\binom{N}{p}^{-1}\sum_{\mathbf{i}\in\mathcal{P}_{p}(\llbracket
N\rrbracket\backslash\\{i\\})}h(X_{\mathbf{i}}).$

In fact, our estimator is closer to Sen’s estimator of the asymptotic variance
of [41, 42] which depends on the kernel size and is related to the jackknife
estimator $\widehat{V}^{h,S}_{N}=(N-k)^{2}/(N-1)^{2}\widehat{V}^{h,J}_{N}$.
However, it is unclear how these estimators could be translated in a two-
dimensional setup where $Y$ is a matrix instead of a vector, especially how to
define the analog of $U_{N}^{h,(-i)}$.

###  4.5 Extension to functions of $U$-statistics

In the case of a function of $U$-statistics
$g(U^{h_{1}}_{N},...,U^{h_{D}}_{N})$, Corollary 3.6 applies and the asymptotic
variance to be estimated is $V^{\delta}=\nabla
g(\theta)^{T}\Sigma^{h_{1},...,h_{D}}\nabla g(\theta)\neq 0$, where
$\theta=(U^{h_{1}}_{\infty},...,U^{h_{D}}_{\infty})$. Similar to the estimator
for the asymptotic variance of $V^{h}$, we suggest an estimator for the
covariance matrix
$\Sigma^{h_{1},...,h_{D}}=\left(C^{h_{k},h_{\ell}}\right)_{1\leq k,\ell\leq
D}$.

For each kernel $h_{k}$, $1\leq k\leq D$, let $\widehat{\mu}^{h_{k},(i)}_{N}$
and $\widehat{\nu}^{h_{k},(j)}_{N}$ be the respective estimators of the
conditional expectations $\psi^{1,0}_{(\\{i\\},\emptyset)}h_{k}$ and
$\psi^{0,1}_{(\emptyset,\\{j\\})}h_{k}$, as defined in equations (8) and (9).

Now define

$\widehat{c}^{h_{k},h_{\ell};1,0}_{N}:=\binom{m_{N}}{2}^{-1}\sum_{1\leq
i_{1}<i_{2}\leq
m_{N}}\frac{(\widehat{\mu}^{h_{k},(i_{1})}_{N}-\widehat{\mu}^{h_{k},(i_{2})}_{N})(\widehat{\mu}^{h_{\ell},(i_{1})}_{N}-\widehat{\mu}^{h_{\ell},(i_{2})}_{N})}{2}$

and

$\widehat{c}^{h_{k},h_{\ell};0,1}_{N}:=\binom{n_{N}}{2}^{-1}\sum_{1\leq
j_{1}<j_{2}\leq
n_{N}}\frac{(\widehat{\nu}^{h_{k},(j_{1})}_{N}-\widehat{\nu}^{h_{k},(j_{2})}_{N})(\widehat{\nu}^{h_{\ell},(j_{1})}_{N}-\widehat{\nu}^{h_{\ell},(j_{2})}_{N})}{2}.$

Then, for two kernels $h_{k}$ and $h_{\ell}$,

$\widehat{C}^{h_{k},h_{\ell}}_{N}:=\frac{p^{2}}{\rho}\widehat{c}^{h_{k},h_{\ell};1,0}_{N}+\frac{q^{2}}{1-\rho}\widehat{c}^{h_{k},h_{\ell};0,1}_{N},$

is an estimator of the asymptotic covariance term $C^{h_{k},h_{\ell}}$.

With a similar proof as for Theorem 4.4, the following theorem ensures the
consistency of this estimator.

######  Theorem 4.9.

For two kernel functions $h_{k}$ and $h_{\ell}$, we have
$\widehat{c}^{h_{k},h_{\ell};1,0}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}c^{1,0}_{h_{k},h_{\ell}}$
and
$\widehat{c}^{h_{k},h_{\ell};0,1}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}c^{0,1}_{h_{k},h_{\ell}}$.
As a consequence,
$\widehat{C}^{h_{k},h_{\ell}}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}C^{h_{k},h_{\ell}}$.

Thus, for linearly independent kernel functions $(h_{1},...,h_{D})$, the
entries of the matrix
$\widehat{\Sigma}^{h_{1},...,h_{D}}_{N}:=\left(\widehat{C}^{h_{k},h_{\ell}}_{N}\right)_{1\leq
k,\ell\leq D}$ converge to the entries of $\Sigma^{h_{1},...,h_{D}}$. Set
$\widehat{V}^{\delta}_{N}:=\nabla
g(U^{h_{1}}_{N},...,U^{h_{D}}_{N})^{T}\widehat{\Sigma}^{h_{1},...,h_{D}}_{N}\nabla
g(U^{h_{1}}_{N},...,U^{h_{D}}_{N})$, then we have the following
straightforward result.

######  Corollary 4.10.

We have
$\widehat{V}^{\delta}_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}V^{\delta}$.

Finally, the plug-in corollary also stands.

######  Corollary 4.11.

If $V^{\delta}>0$, then

$\sqrt{\frac{N}{\widehat{V}^{\delta}_{N}}}\left(g(U^{h_{1}}_{N},...,U^{h_{D}}_{N})-g(\theta)\right)\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1).$

##  5 RCE models, kernel functions and network comparison

###  5.1 Examples of RCE models

#### Bipartite Expected Degree Distribution model

The Bipartite Expected Degree Distribution (BEDD) model, suggested by [36], is
a binary graph model characterised by two distributions from which the row and
column nodes draw a weight. The probability of a connection between two nodes
is fully determined by the corresponding row and column weight distributions.
The distribution of a graph following a BEDD model can be written using latent
variables $(\xi_{i})_{i\geq 1}$ and $(\eta_{j})_{j\geq 1}$ corresponding to
the row and column nodes of the graph :

$\begin{split}\xi_{i},\eta_{j}&\overset{iid}{\sim}\mathcal{U}[0,1]\\\
Y_{ij}~{}\mid~{}\xi_{i},\eta_{j}&\sim\mathcal{B}(\lambda
f(\xi_{i})g(\eta_{j})).\end{split}$ (10)

where

  * •

$f$ and $g$ are positive, càdlàg, nondecreasing, bounded and normalized ($\int
f=\int g=1$) real functions $[0,1]\rightarrow\mathbb{R}_{+}$,

  * •

$\lambda$ is a positive real number such that $\lambda\leq\lVert
f\rVert^{-1}_{\infty}\lVert g\rVert^{-1}_{\infty}$,

  * •

$\mathcal{B}$ is the Bernoulli distribution.

The BEDD model is a RCE version of the Expected Degree Sequence model of [12]
but in the BEDD, the row weights $f(\xi_{i})$ and column weights $g(\eta_{j})$
are exchangeable. Indeed, $f$ and $g$ characterize the weight distributions of
the row and column nodes whereas the weights are fixed in the Expected Degree
Sequence model. The triplet $(\lambda,f,g)$ is called the BEDD parameters.

#### Latent Block model

The Latent Block model (LBM) [21] is a binary graph model characterised by a
partition of row and column nodes in several groups. It can be considered as a
bipartite extension of the Stochastic Block model [34]. The probability of
interaction between two nodes is fully determined by the groups to which they
belong. All the nodes have the same probability to belong to each group. The
distribution of a LBM is most commonly written using independent latent
variables for the node attribution in a group $(Z_{i})_{i\geq 1}$ and
$(W_{j})_{j\geq 1}$ corresponding to the row and column nodes of the graph :

$\begin{split}Z_{i}&\overset{iid}{\sim}\mathcal{M}(1;\boldsymbol{\alpha})\\\
W_{j}&\overset{iid}{\sim}\mathcal{M}(1;\boldsymbol{\beta})\\\
Y_{ij}~{}\mid~{}Z_{i}=k,W_{j}=\ell&\sim\mathcal{B}(\pi_{k\ell}),\end{split}$
(11)

where $\boldsymbol{\alpha}=(\alpha_{1},...,\alpha_{K})$ and
$\boldsymbol{\beta}=(\beta_{1},...,\beta_{L})$ are the probability vectors of
the rows and the columns and $\boldsymbol{\pi}=(\pi_{k\ell})_{1\leq k\leq
K,1\leq\ell\leq L}\in[0,1]^{KL}$ is a matrix of probabilities.

The LBM is a RCE model, since the group attribution variables of the nodes are
exchangeable.

#### W-graph model

Let $w$ be a function of $[0,1]^{2}\rightarrow[0,1]$. The W-graph model
associated to $w$ is defined by

$\begin{split}\xi_{i},\eta_{j}&\overset{i.i.d.}{\sim}\mathcal{U}[0,1]\\\
Y_{ij}~{}\mid~{}\xi_{i},\eta_{j}&\sim\mathcal{B}\left(w(\xi_{i},\eta_{j})\right),\end{split}$
(12)

$w$ is sometimes referred to as a graphon. For identification reasons, we
assume $\int w(\cdot,\eta)d\eta$ and $\int w(\xi,\cdot)d\xi$ to be càdlàg,
nondecreasing and bounded.

Any RCE model can be written as a W-graph model [14] so can the BEDD model and
the LBM be expressed with a graphon. For the BEDD model,
$w(\xi_{i},\eta_{j})=\lambda f(\xi_{i})g(\eta_{j})$ where $f$ and $g$ are
those of Equation (10). For the LBM,
$w(\xi_{i},\eta_{j})=\sum_{k,\ell}\pi_{k\ell}\mathds{1}\\{s(\xi_{i})=k\\}\mathds{1}\\{t(\eta_{j})=\ell\\}$
where
$s(\xi_{i})=1+\sum_{k=1}^{K}\mathds{1}\\{\xi_{i}>\sum_{k^{\prime}=1}^{k}\alpha_{k^{\prime}}\\}$,
$t(\eta_{j})=1+\sum_{\ell=1}^{L}\mathds{1}\\{\eta_{j}>\sum_{\ell^{\prime}=1}^{\ell}\beta_{\ell^{\prime}}\\}$
and $\boldsymbol{\alpha}$, $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$ are
those of Equation (11).

#### Extension to weighted graphs

All the models presented above are defined for binary graphs. However, one can
extend it to weighted graphs by switching the Bernoulli distributions with
another. In particular, the Poisson distribution is particularly fitted for
count values in $\mathbb{N}$ :

$\begin{split}\xi_{i},\eta_{j}&\overset{i.i.d.}{\sim}\mathcal{U}[0,1]\\\
Y_{ij}~{}\mid~{}\xi_{i},\eta_{j}&\sim\mathcal{P}\left(w(\xi_{i},\eta_{j})\right),\end{split}$
(13)

where $w:[0,1]^{2}\rightarrow\mathbb{R}_{+}$ is the graphon of this model.
Similarly, one can define a Poisson-BEDD model or Poisson-LBM.

###  5.2 Examples of kernel functions

We present three examples of statistics which are $U$-statistics or functions
of $U$-statistics. Sometimes, a kernel function $h$ has a long expression,
especially when it is the symmetric version of some simpler function $h^{0}$,
as in (3). In this case, we introduce the kernels of interest with $h^{0}$
instead of $h$, but the $U$-statistics and asymptotic results always apply to
the symmetric version $h$.

#### Motif counts

Motifs are the name given to small subgraphs. Their occurrences in the
complete network can be counted. Motif counts are useful statistics for random
graphs as they provide information on the network local structure [43, 26].
Many random network models hinge on motif frequencies. In the Exponential
Random Graph Model [19], motif frequencies are sufficient statistics. The
$dk$-random graph model [35] also largely relies on motif frequencies.

Their asymptotic properties are widely studied and a large numbers of studies
use motif counts to perform statistical tests [39, 8, 7, 13, 20, 30, 32, 36].
Our framework is particularly well adapted to the use of motif counts for
statistical tests as frequencies are in fact $U$-statistics with kernel
functions of the same size as the motifs.

For many applications, motifs are considered as elementary building blocks
which can be specifically interpreted. This is the case in molecular biology
[43, 37, 1], neurology [52], sociology [6], evolution [38] and ecology [48,
4]. In particular, Figure 3 of [45] lists all the bipartite motifs consisting
of from 2 to 6 nodes. For example, their motif 6 represents a $2\times 2$
clique, where every node is connected to all others (Figure 1). Their motif 14
represents a path between two column nodes, passing through another column
node and two row nodes (Figure 1). The latter indicates an indirect
interaction between the row nodes, hinging on the middle row node. [27] found
that motif 6 is over-represented in plant-pollinator interaction networks,
while motif 14 is under-represented, compared to Erdös-Rényi graphs of the
same density.

Figure 1: Motif 6 and 14, counted by $U^{h_{6}}_{N}$ and $U^{h_{14}}_{N}$

These motifs can be counted with $U^{h_{6}}_{N}$ and $U^{h_{14}}_{N}$ the
$U$-statistics associated to the kernels $h_{6}$ and $h_{14}$, that we
introduce as the symmetric version of the following functions :

$\displaystyle h^{0}_{6}(Y_{(i_{1},i_{2};j_{1},j_{2})})$
$\displaystyle=Y_{i_{1}j_{1}}Y_{i_{1}j_{2}}Y_{i_{2}j_{1}}Y_{i_{2}j_{2}},$
$\displaystyle h^{0}_{14}(Y_{(i_{1},i_{2};j_{1},j_{2},j_{3})})$
$\displaystyle=Y_{i_{1}j_{1}}Y_{i_{1}j_{2}}Y_{i_{2}j_{2}}Y_{i_{2}j_{3}}(1-Y_{i_{2}j_{1}})(1-Y_{i_{1}j_{3}}).$
(14)

######  Remark.

With this example, it is very apparent that introducing simple, non-symmetric
kernel functions and then symmetrizing them is way simpler than introducing
symmetric kernel functions. $h^{0}_{6}$ is already symmetric, so
$h_{6}=h^{0}_{6}$. However, $h^{0}_{14}$ is not symmetric and one can sum over
all the permutations of the indices to make it symmetric using formula (3).
Because of the automorphisms of motif 14, it involves only $6$ permutations
instead of $2!3!=12$.

$\begin{split}h_{14}(Y_{(i_{1},i_{2};j_{1},j_{2},j_{3})})&=\frac{1}{6}Y_{i_{1}j_{1}}Y_{i_{1}j_{2}}Y_{i_{2}j_{2}}Y_{i_{2}j_{3}}(1-Y_{i_{2}j_{1}})(1-Y_{i_{1}j_{3}})+\frac{1}{6}Y_{i_{1}j_{2}}Y_{i_{1}j_{3}}Y_{i_{2}j_{3}}Y_{i_{2}j_{1}}(1-Y_{i_{1}j_{1}})(1-Y_{i_{2}j_{2}})\\\
&+\frac{1}{6}Y_{i_{1}j_{3}}Y_{i_{1}j_{1}}Y_{i_{2}j_{1}}Y_{i_{2}j_{2}}(1-Y_{i_{1}j_{2}})(1-Y_{i_{2}j_{3}})+\frac{1}{6}Y_{i_{2}j_{1}}Y_{i_{2}j_{2}}Y_{i_{1}j_{2}}Y_{i_{1}j_{3}}(1-Y_{i_{1}j_{1}})(1-Y_{i_{2}j_{3}})\\\
&+\frac{1}{6}Y_{i_{2}j_{2}}Y_{i_{2}j_{3}}Y_{i_{1}j_{3}}Y_{i_{1}j_{1}}(1-Y_{i_{2}j_{1}})(1-Y_{i_{1}j_{2}})+\frac{1}{6}Y_{i_{2}j_{3}}Y_{i_{2}j_{1}}Y_{i_{1}j_{1}}Y_{i_{1}j_{2}}(1-Y_{i_{2}j_{2}})(1-Y_{i_{1}j_{3}}).\end{split}$

Using Corollary 4.7, the following studentized statistics converge to a
standard normal distribution

$Z^{6}_{N}:=\sqrt{\frac{N}{\widehat{V}^{h_{6}}_{N}}}\left(U^{h_{6}}_{N}-U^{h_{6}}_{\infty}\right)\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1),$
$Z^{14}_{N}:=\sqrt{\frac{N}{\widehat{V}^{h_{14}}_{N}}}\left(U^{h_{14}}_{N}-U^{h_{14}}_{\infty}\right)\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1).$

#### Product graphon

The $W$-graph model encompasses all the dissociated RCE models. One can make
use of this model to compare the form of the graphon that has generated a
graph to a known form. One interesting form is the product form, which
corresponds to an absence of specific interaction between row nodes and column
nodes. For example, one can normalize the graphon in the Poisson $W$-graph
model described by (13) such that

$Y_{ij}~{}\mid~{}\xi_{i},\eta_{j}\sim\mathcal{P}(\lambda\bar{w}(\xi_{i},\eta_{j})),$

where $\lambda>0$ and $\iint\bar{w}=1$. Call $\bar{w}$ a normalized graphon.
Define $f$ and $g$ as the marginals of $\bar{w}$, i.e.

$\displaystyle f$ $\displaystyle=\int\bar{w}(\cdot,\eta)d\eta$ $\displaystyle
g$ $\displaystyle=\int\bar{w}(\xi,\cdot)d\xi.$ (15)

A density-free dissimilarity measure between a graphon and its corresponding
product form could be the quantity $d(\bar{w})$ defined as

$d(\bar{w}):=\lvert\lvert\bar{w}-fg\rvert\rvert^{2}_{2}=\iint\left(\bar{w}(\xi,\eta)-f(\xi)g(\eta)\right)^{2}d\xi
d\eta.$ (16)

Observe that $d(\bar{w})=0$ if and only if $\bar{w}$ is of product form, e.g.
$\bar{w}(\xi,\eta)=f(\xi)g(\eta)$ almost everywhere in $[0,1]^{2}$. In this
case, the $W$-graph model (5.2) is a BEDD model and $f$ and $g$ have the same
role as in (10).

We suggest an estimator for $d(\bar{w})$. Let $h_{A},h_{B},h_{C}$ and $h_{D}$
be kernel functions of respective size $2\times 2$, $1\times 2$, $2\times 1$
and $1\times 1$ defined as in Table 1. These kernel functions are linearly
independent, so Corollary 3.5 applies to the associated $U$-statistics
$(U^{h_{A}}_{N},U^{h_{B}}_{N},U^{h_{C}}_{N},U^{h_{D}}_{N})$, and they are
jointly asymptotically normal with asymptotic covariance matrix
$\Sigma^{h_{A},h_{B},h_{C},h_{D}}$.

$h$ |  $h^{0}(Y_{(i_{1},...,i_{r};j_{1},...,j_{c})})$ |  $\mathbb{E}[h(Y_{(i_{1},...,i_{r};j_{1},...,j_{c})})]$
---|---|---
$h_{A}$ |  $h^{0}_{A,1}(Y_{(i_{1},i_{2};j_{1},j_{2})})-2h^{0}_{A,2}(Y_{(i_{1},i_{2};j_{1},j_{2})})$ |  $\lambda^{3}\iint\bar{w}(\xi,\eta)(\bar{w}(\xi,\eta)-2f(\xi)g(\eta))d\xi d\eta$
$h_{A,1}$ |  $Y_{i_{1}j_{1}}(Y_{i_{1}j_{1}}-1)Y_{i_{2}j_{2}}$ |  $\lambda^{3}\iint\bar{w}(\xi,\eta)^{2}d\xi d\eta$
$h_{A,2}$ |  $Y_{i_{1}j_{1}}Y_{i_{1}j_{2}}Y_{i_{2}j_{2}}$ |  $\lambda^{3}\iint\bar{w}(\xi,\eta)f(\xi)g(\eta)d\xi d\eta$
$h_{B}$ |  $Y_{i_{1}j_{1}}Y_{i_{1}j_{2}}$ |  $\lambda^{2}\int f(\xi)^{2}d\xi$
$h_{C}$ |  $Y_{i_{1}j_{1}}Y_{i_{2}j_{1}}$ |  $\lambda^{2}\int g(\eta)^{2}d\eta$
$h_{D}$ |  $\displaystyle{Y_{i_{1}j_{1}}}$ |  $\lambda$
Table 1: Kernel functions used for the estimation of $d(\bar{w})$. The
expectations are given by Lemma E.1 in Appendix.

We define the estimator of $d(\bar{w})$ as
$\widehat{d}_{N}=t(U^{h_{A}}_{N},U^{h_{B}}_{N},U^{h_{C}}_{N},U^{h_{D}}_{N}):=U^{h_{A}}_{N}/(U^{h_{D}}_{N})^{3}+U^{h_{B}}_{N}U^{h_{C}}_{N}/(U^{h_{D}}_{N})^{4}$.
We have
$t(U^{h_{A}}_{\infty},U^{h_{B}}_{\infty},U^{h_{C}}_{\infty},U^{h_{D}}_{\infty})=\iint\bar{w}(\xi,\eta)(\bar{w}(\xi,\eta)-2f(\xi)g(\eta))d\xi
d\eta+\int f(\xi)^{2}d\xi\times\int g(\eta)^{2}d\eta=d(\bar{w})$. Then,
Corollary 4.11 ensures that the studentized statistic converge to a standard
normal distribution

$\begin{split}Z^{d}_{N}:=\sqrt{\frac{N}{\widehat{V}^{d}_{N}}}\left(\widehat{d}_{N}-d(\bar{w})\right)\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1).\end{split}$

where $\widehat{V}^{d}_{N}=\nabla
t(U^{h_{A}}_{N},U^{h_{B}}_{N},U^{h_{C}}_{N},U^{h_{D}}_{N})^{T}\widehat{\Sigma}^{h_{A},h_{B},h_{C},h_{D}}_{N}\nabla
t(U^{h_{A}}_{N},U^{h_{B}}_{N},U^{h_{C}}_{N},U^{h_{D}}_{N})$ and
$\widehat{\Sigma}^{h_{A},h_{B},h_{C},h_{D}}_{N}$ is defined as in Corollary
4.11.

#### Heterogeneity of the rows of a network

The degree distributions of a network hold significant information about its
topology. For a binary network, denote $D_{i}:=\sum_{j}Y_{ij}$ the degree of
the $i$-th row, then
$\mathbb{E}[D_{i}\mid\xi_{i}]=f(\xi_{i})\mathbb{E}[D_{i}]$, where $f$ is given
by (15). Therefore, the marginals $f$ and $g$ accounts for the expected
relative degree distributions of the binary network, which is known to
characterize networks. In the BEDD model, these distributions (and the density
of the network) fully characterize the model.

Although the sum of the edge weights and the number of edges stemming from a
nodes are equivalent for binary networks, for weighted networks, they are two
different quantities. The sum of the edge weights stemming from a node is
sometimes called its strength [5], which do not depend entirely on the number
of edges. In a network with weighted edges, $f(\xi_{i})$ corresponds to
expected relative strength of the $i$-th row node, instead of its expected
relative degree.

One way to characterize the degree/strength distributions of a network is to
calculate their variance. In particular, the variance of the degree
distribution of a network quantifies its heterogeneity (i.e. the unbalance
between strongly interacting nodes and the others) and can be used as an index
to characterize networks [46]. In our framework with random graph models,
rather than directly study the empirical distribution of row
degrees/strengths, one would like to retrieve information on $f$ and $g$, the
distributions of the expected degrees/strengths specified by the model. For
the BEDD model, $f$ and $g$ are directly given by the model. The variance of
the row expected relative degree/strength distribution is $F_{2}-1$, where
$F_{2}:=\int f^{2}(\xi)d\xi$.

$F_{2}$ can be estimated using the estimator
$\widehat{F}_{2,N}:=\kappa(U^{h_{1}}_{N},U^{h_{2}}_{N})=U^{h_{1}}_{N}/U^{h_{2}}_{N}$,
using the $U$-statistics associated to the submatrix kernels functions $h_{1}$
and $h_{2}$ defined as

$h_{1}(Y_{(i;j_{1},j_{2})})=Y_{ij_{1}}Y_{ij_{2}}$

and

$h_{2}(Y_{(i_{1},i_{2};j_{1},j_{2})})=\frac{1}{2}(Y_{i_{1}j_{1}}Y_{i_{2}j_{2}}+Y_{i_{1}j_{2}}Y_{i_{2}j_{1}}).$

By Lemma E.2 in Appendix, we have $U^{h_{1}}_{\infty}=\lambda^{2}F_{2}$ and
$U^{h_{2}}_{\infty}=\lambda^{2}$. Given that
$\nabla\kappa(U^{h_{1}}_{N},U^{h_{2}}_{N})=(1/U^{h_{2}}_{N},-U^{h_{1}}_{N}/(U^{h_{2}}_{N})^{2})$,
the application of Corollary 4.11 yields

$Z^{F_{2}}_{N}:=\sqrt{\frac{N}{\widehat{V}^{F_{2}}_{N}}}\bigg{(}\widehat{F}_{2,N}-F_{2}\bigg{)}\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1),$
(17)

where

$\widehat{V}^{F_{2}}_{N}=\frac{1}{(U^{h_{2}}_{N})^{2}}\widehat{V}^{h_{1}}_{N}-\frac{2\widehat{F}_{2,N}}{(U^{h_{2}}_{N})^{2}}\widehat{C}^{h_{1},h_{2}}_{N}+\frac{(\widehat{F}_{2,N})^{2}}{(U^{h_{2}}_{N})^{2}}\widehat{V}^{h_{2}}_{N}.$

In a similar way, one can define the analogous estimator $\widehat{G}_{2,N}$
for $G_{2}:=\int g^{2}(\eta)d\eta$ and its associated variance estimator to
quantify the heterogeneity of the columns of networks.

###  5.3 Network comparison

Network comparison has a long history in network analysis literature [18, 49].
Few comparison methods use random network models [3, 30, 29]. However, there
are several advantages to using model-based approaches. Indeed, random models
define a probability distribution of networks which can be used to derive
statistical guarantees. Also, models can be used to control the sources of
heterogeneity in networks. The comparison can be made with respect to a
quantity of interest, which makes it easier to interpret. Now, we show how to
define a test statistic for model-based network comparison with our framework.
Naturally, $U$-statistics define network statistics characterizing a network.
They can be used to analyze a single network, but their use is easily extended
to compare different networks.

Let $Y^{A}$ and $Y^{B}$ be two independent networks of respective sizes
$m^{A}_{N}\times n^{A}_{N}$ and $m^{B}_{N}\times n^{B}_{N}$. Define a network
quantity of interest $\theta$. These two networks are generated by two models,
leading to different values $\theta^{A}$ and $\theta^{B}$ of this quantity for
the two models. Let $\widehat{\theta}_{N}$ be an estimator for this quantity
of interest. If $\widehat{\theta}_{N}$ is a $U$-statistic or a function of
$U$-statistics, then as previously, we have both

$\sqrt{N}(\widehat{\theta}_{N}(Y^{A})-\theta^{A})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,V^{A})$

and

$\sqrt{N}(\widehat{\theta}_{N}(Y^{B})-\theta^{B})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,V^{B}),$

where $V^{A}$ and $V^{B}$ are the asymptotic variances.

To compare $\theta^{A}$ and $\theta^{B}$, we can confront the test hypotheses
$\mathcal{H}_{0}:\theta^{A}=\theta^{B}$ and
$\mathcal{H}_{1}:\theta^{A}\neq\theta^{B}$. Because $Y^{A}$ and $Y^{B}$ are
independent, the previous convergence results give

$\sqrt{\frac{N}{V^{A}+V^{B}}}\bigg{(}\widehat{\delta}_{N}(Y^{A},Y^{B})-(\theta^{A}-\theta^{B})\bigg{)}\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1),$

where
$\widehat{\delta}_{N}(Y^{A},Y^{B}):=\widehat{\theta}_{N}(Y^{A})-\widehat{\theta}_{N}(Y^{B})$.
Therefore, using the consistent estimators $\widehat{V}^{A}_{N}$ and
$\widehat{V}^{B}_{N}$ of the asymptotic variances, the test statistic

$Z_{N}(Y^{A},Y^{B}):=\sqrt{\frac{N}{\widehat{V}^{A}_{N}+\widehat{V}^{B}_{N}}}\widehat{\delta}_{N}(Y^{A},Y^{B})$

admits the convergence result

$Z_{N}(Y^{A},Y^{B})-\sqrt{\frac{N}{\widehat{V}^{A}_{N}+\widehat{V}^{B}_{N}}}(\theta^{A}-\theta^{B})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,1).$

Under $\mathcal{H}_{0}$, we have $\theta^{A}-\theta^{B}=0$, so
$Z_{N}(Y^{A},Y^{B})$ converges in distribution to a standard Gaussian
variable.

##  6 Simulations

In this section, we illustrate our theoretical results, in particular
Corollaries 4.7 and 4.11, using the previous examples of graph models and
kernel functions. For each example of kernel function, we check the asymptotic
normality of the studentized statistic for networks simulated under different
configurations (different models, values of $N$ and $\rho$). We give the
resulting Q-Q plots. We also examine the coverage probabilities of the
confidence intervals built with the estimates or the standard deviations of
the statistics, depending on whether they are fitted for estimation or
statistical testing.

###  6.1 Motif counts

#### Model I

This graph model is a LBM with 2 row groups and 2 columns groups of equal
proportion. Using the notations of (11), this means
$\boldsymbol{\alpha}=\boldsymbol{\beta}=(0.5,0.5)$. The probability matrix
$\boldsymbol{\pi}$ has size $2\times 2$. We set $\pi_{k\ell}=0.5$ for all
$1\leq k,\ell\leq 2$ except $\pi_{11}=0.95$.

Under Model I, for each value
$N\in\\{32,45,64,90,128,181,256,362,512,724,1024,1448,2048\\}$ and
$\rho\in\\{1/8,1/2\\}$, we have simulated $K=500$ networks of size
$m_{N}\times n_{N}$ where $m_{N}=\lfloor\rho N\rfloor$ and $n_{N}=N-m_{N}$.
The relative frequencies of motifs 6 and 10 of [45] are respectively given by
$U^{h_{6}}_{N}$ and $U^{h_{14}}_{N}$ where $h_{6}$ and $h_{14}$ were defined
by (14).

The Q-Q plots of the studentized statistic associated to $h_{6}$ and $h_{14}$
are given in Figures 2 and 3. For $\rho=1/2$, we observe that the empirical
distribution of both statistics converge and become close to a normal
distribution as $N\gtrsim 128$. Figure 4 gives the respective coverage
probabilities for $U^{h_{6}}_{\infty}$ and $U^{h_{1}4}_{\infty}$ of the
intervals
$\left[U^{h_{6}}_{N}-\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{h_{6}}_{N}}},U^{h_{6}}_{N}+\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{h_{6}}_{N}}}\right]$
and
$\left[U^{h_{14}}_{N}-\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{h_{14}}_{N}}},U^{h_{14}}_{N}+\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{h_{14}}_{N}}}\right]$
respectively, where $\Phi$ is the quantile function of the standard normal
distribution. The coverage probabilities converge to $\alpha$ but with
different speeds depending on the motif. We also observe that a larger number
of nodes is needed to reach the target coverage probabilities with $\rho=1/8$
(rectangular matrix) than $\rho=1/2$ (square matrix).

Figure 2: Q-Q plots for $Z^{6}_{N}$ the studentized statistic associated with
$U^{h_{6}}_{N}$, with $\rho=0.5$. Figure 3: Q-Q plots for $Z^{14}_{N}$ the
studentized statistic associated with $U^{h_{14}}_{N}$, with $\rho=0.5$.
Figure 4: Empirical coverage probabilities for the asymptotic confidence
intervals at level $\alpha=0.95$ of $U^{h_{6}}_{N}$ (left) and
$U^{h_{14}}_{N}$ (right) for different values of $N$ (x-axis),
$\rho\in\\{1/8,1/2\\}$. Grey dashed lines represent the confidence interval at
level $0.95$ of the frequency $Z=X/K$, if $X$ follows the binomial
distribution with parameters $K$ and $\alpha=0.95$.

###  6.2 Graphon product distance

#### Model II($\epsilon$)

We consider a Poisson-LBM. As this is a weighted graph model,
$\boldsymbol{\pi}$ is a matrix of weights rather than probabilities, i.e.
$(\pi_{k\ell})_{1\leq k\leq K,1\leq\ell\leq L}$ are real non-negative numbers.
We consider 2 row groups and 2 column groups of equal proportion. Let
$\boldsymbol{\pi}^{0}$ be a weight matrix. We set $\pi^{0}_{11}=4$,
$\pi^{0}_{12}=\pi^{0}_{21}=2$ and $\pi^{0}_{22}=1$. The LBM with parameters
$\boldsymbol{\alpha}$, $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$ is also a
BEDD model. Indeed, its graphon has a product form and can be written
$w_{0}(\xi,\eta)=\lambda\bar{w}_{0}(\xi,\eta)=\lambda f_{0}(\xi)g_{0}(\eta)$
where $\lambda=9/4$, $f_{0}=g_{0}$ both take values $4/3$ on $[0,0.5]$, $2/3$
on $]0.5,1]$.

Now, let $\boldsymbol{\tau}$ be the $2\times 2$ matrix where
$\tau_{11}=\tau_{22}=2$ and $\tau_{12}=\tau_{21}=0$. For $\epsilon\geq 0$, we
define Model II($\epsilon$) as a LBM with group probabilities
$\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ and with weight matrix
$\boldsymbol{\pi}^{\epsilon}=\frac{\lambda}{\lambda+\epsilon}(\boldsymbol{\pi}^{0}+\epsilon\boldsymbol{\tau})$.
Thus, Model II(0) is a LBM with product form (BEDD model) and for
$\epsilon>0$, Model II($\epsilon$) with graphon
$w_{\epsilon}=\lambda\bar{w}_{\epsilon}$ is a perturbed version and rescaled
to preserve the same density $\lambda$. As $\epsilon$ grows, the graphon of
Model II($\epsilon$) strays further from a BEDD model. Indeed, one can show
that
$d(\bar{w}_{\epsilon})=64\epsilon^{2}(5+2\epsilon)^{2}/(9+4\epsilon)^{4}$,
which is an increasing function for $\epsilon\geq 0$.

Under Model II($\epsilon$), for each value
$N\in\\{32,45,64,90,128,181,256,362,512,724,1024,1448,2048\\}$,
$\rho\in\\{1/8,1/2\\}$, $\epsilon\in\\{0.5,1,1.5,2,2.5,3\\}$, we have
simulated $K=500$ networks of size $m_{N}\times n_{N}$ where
$m_{N}=\lfloor\rho N\rfloor$ and $n_{N}=N-m_{N}$. For each network, we have
computed $\widehat{d}_{N}$ and $\widehat{V}^{d}_{N}$ as estimates of
$d(\bar{w}_{\epsilon})$ and $V^{d}$ respectively. The Q-Q plots of the
studentized statistic $Z^{d}_{N}$ are given in Figures 5. It is apparent that
$Z^{d}_{N}$ is not centered, especially for $N\lesssim 512$. This is due to
the fact that $\widehat{d}_{N}$ is obtained via the delta method, so it is a
biased estimator of $d(\bar{w}_{\epsilon})$ for finite values of $N$. However,
the bias converges to $0$ when $N$ grows, so we find that the statistic
achieves normality for $N\gtrsim 2048$. Figure 6 gives for different values of
$\epsilon$, the coverage probability for $d(\bar{w}_{\epsilon})$ of the
interval
$\left[\widehat{d}_{N}-\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{d}_{N}}},\widehat{d}_{N}+\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{d}_{N}}}\right]$
where $\Phi$ is the quantile function of the standard normal distribution. We
observe that convergence is fastest when $\rho=1/2$ but also when $\epsilon$
is larger, which seems to indicate that estimation of $d(\bar{w})$ is more
precise for square matrices and when $w$ is further away from a product model.

Figure 5: Q-Q plot for $Z^{d}_{N}$ the studentized statistic associated with
$\widehat{d}_{N}$, $\epsilon=3$, $\rho=0.5$ Figure 6: Empirical coverage
probabilities for the asymptotic confidence intervals at level $\alpha=0.95$
of $\widehat{d}_{N}$ for different values of $N$ (x-axis),
$\rho\in\\{1/8,1/2\\}$, $\epsilon\in\\{1,2,3\\}$. Grey dashed lines represent
the confidence interval at level $0.95$ of the frequency $Z=X/K$, if $X$
follows the binomial distribution with parameters $K$ and $\alpha=0.95$.

###  6.3 Heterogeneity in the row weights of a network

#### Model III

In this example, we consider a weighted BEDD model with power-law strength
distributions, i.e. the marginals $f$ and $g$ have the form
$f(\xi)=(\alpha_{f}+1)\xi^{\alpha_{f}}$ and
$g(\eta)=(\alpha_{g}+1)\eta^{\alpha_{g}}$, where $\alpha_{f}$ and $\alpha_{g}$
are real non-negative numbers. $\alpha_{f}$ is directly related to $F_{2}=\int
f(\xi)^{2}d\xi=(\alpha_{f}+1)^{2}/(2\alpha_{f}+1)$.

Under Model III, for each value
$N\in\\{32,45,64,90,128,181,256,362,512,724,1024,1448,2048\\}$ and
$\rho\in\\{1/8,1/2,7/8\\}$, we have simulated $K=500$ networks of size
$m_{N}\times n_{N}$ where $m_{N}=\lfloor\rho N\rfloor$ and $n_{N}=N-m_{N}$. We
estimate $F_{2}$ and $G_{2}$ using $\widehat{F}_{2,N}$ and its symmetric
counterpart $\widehat{G}_{2,N}$. We also compute $\widehat{V}^{F_{2}}_{N}$ and
$\widehat{V}^{G_{2}}_{N}$ to obtain the studentized statistics $Z^{F_{2}}_{N}$
and $Z^{G_{2}}_{N}$. Q-Q plots of $Z^{F_{2}}_{N}$ are given in Figures 7.
Figure 8 gives the respective coverage probabilities for $F_{2}$ and $G_{2}$
of the intervals
$\left[\widehat{F}_{2,N}-\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{F_{2}}_{N}}},\widehat{F}_{2,N}+\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{F_{2}}_{N}}}\right]$
and
$\left[\widehat{G}_{2,N}-\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{G_{2}}_{N}}},\widehat{G}_{2,N}+\Phi(1-\frac{\alpha}{2})\sqrt{\frac{N}{\widehat{V}^{G_{2}}_{N}}}\right]$
respectively, where $\Phi$ is the quantile function of the standard normal
distribution. Despite the bias due to the delta method, apparent for $N=64$ in
Figure 7, we find that the coverage probabilities fall in the confidence
intervals even for small values of $N$. We would expect that for $F_{2}$, the
worst case correspond to $\rho=1/8$ (there are less rows) and that for
$G_{2}$, the worst case correspond to $\rho=7/8$ (there are less columns), but
Figure 8 does not show a clear difference between the three values of $\rho$.

Figure 7: Q-Q plot for $Z^{F_{2}}_{N}$ the studentized statistic associated
with $\widehat{F}_{2,N}$, $\rho=0.5$ Figure 8: Empirical coverage
probabilities for the asymptotic confidence intervals at level $\alpha=0.95$
of $Z^{F_{2}}_{N}$ (left) and $Z^{G_{2}}_{N}$ (right) for different values of
$N$ (x-axis), $\rho\in\\{1/8,1/2,7/8\\}$. Grey dashed lines represent the
confidence interval at level $0.95$ of the frequency $Z=X/K$, if $X$ follows
the binomial distribution with parameters $K$ and $\alpha=0.95$.

##  7 Illustrations

To illustrate the use and interpretation of some of the $U$-statistics
introduced in this paper, we considered the set of law-makers networks
compiled by [31]. The database contains networks arising from different fields
(ecology, social sciences, life sciences). We focused on the subset of so-
called ’legislature’ networks both because of their sizes and because network
comparison is of interest for this dataset.

#### Data description.

Four law-maker assemblies were considered: the European Parliament (’EP’), the
General Assembly of the United Nations (’UN’), the US House of Representatives
(’USH’) and the US Senate (’USS’). One network has been recorded each year for
each parliament; we considered the 26 years from 1979 to 2004, for which the
data are available for all the four assemblies. The network recorded for a
given assembly in a given year consists of the votes (yes or no) of the
different members (rows of the adjacency matrix) for the different proposed
laws (columns of the adjacency matrix).
Figure 9 gives the dimensions and densities of the 26 networks collected in
each assembly: the main difference is that the European Parliament is both
larger (both in terms of members and laws) and sparser than the three others.

Figure 9: Distribution of the number of members (left), number of laws
(center) and density (right) of the fours lawmakers networks across the 26
years (in log-scale).

#### Degree imbalance.

We then focused on the degree of imbalance among the rows (resp. columns),
which, under the weighted BEDD model defined in Equation (10), can be measured
by the $U$-statistic $F_{2}$ (resp. $G_{2}$). Figure 10 gives the evolution of
each of the two indicators along the years for each parliament. We observe
that, for each of them, both $F_{2}$ and $G_{2}$ remain above 1, all along the
period: as expected no uniformity exist, neither among the members
($F_{2}>1$), nor among the laws ($G_{2}>1$). Regarding the US networks (USH
and USS), the imbalanced if more marked among the laws than between the
members. As expected also, the confidence intervals are narrower for the
largest networks (EP). No systematic pattern is observed, except the shift in
the imbalance among resolutions voted at the General Assembly of the United
Nations (UN) that is observed in 1992 (and which happens to coincide with the
UN membership of former soviet republics).

|  |  |
---|---|---|---
|  |  |
Figure 10: Evolution of degree imbalance in each assembly along the years.
Top: $F_{2}$ $U$-statistics, bottom: $G_{2}$. From left to right: European
Parliament, General Assembly of the United Nations, US House of
Representatives and US Senate. Solid line: $U$-statistic as an estimate of
$F_{2}$ (resp. $G_{2}$). Dotted line: 95%-confidence interval for $F_{2}$
(resp. $G_{2}$).

#### Network comparison.

For each available year, we then compared the networks of the four assemblies
in terms of degree imbalance ($F_{2}$ and $G_{2}$) and frequency of
topological motifs 6 (as given in Figure 1). We chose this motif as it
constitutes a clique, characterizing a group behavior, in which close members
are in favor of the same laws. For each of these parameters, we use the
comparison test procedure described in Section 5.3.
Figure 11 display the results. We observe no significant difference between
the two US assemblies, which are also the smallest ones: the absence of
significant differences can therefore result from a weak power of the tests
when considering small networks. We also observe a higher heterogeneity among
the members of the European Parliament with respect to all other assemblies,
as well as a higher heterogeneity among the member of the United Nations
assembly, with respect to the two US chambers. A different picture is obtained
for the heterogeneity among the laws, which is significantly higher in the UN
assembly and significantly lower in the EP assembly.
The frequency of motif 6 is interesting, as it may reveal a specific socio-
political behavior. To this respect, the group structure turns out to be much
stronger in the UN than in the EP, the members of which represent both
different political orientations and different nations.

Figure 11: Network comparison for the $U$-statistics, $F_{2}$, $G_{2}$ and for
the count the motif 6 (M6). EP-UN= —$\medsquare$—, EP-USS= —$\meddiamond$—,
EP-USH= —$\medtriangleup$—, UN-USS= \- -$\meddiamond$\- -, UN-USH= \-
-$\medtriangleup$\- -, USS-USH= $\cdots\medtriangleup\cdots$. Horizontal lines
= standard normal quantiles with level .025 and .975.

#### Acknowledgements.

This work was partially supported by the grant ANR-18-CE02-0010-01 of the
French National Research Agency ANR (project EcoNet).

## References

  * [1] Waqar Ali, Tiago Rito, Gesine Reinert, Fengzhu Sun, and Charlotte M Deane.  Alignment-free protein interaction network comparison.  Bioinformatics, 30(17):i430–i437, 2014.
  * [2] James N Arvesen.  Jackknifing u-statistics.  The Annals of Mathematical Statistics, 40(6):2076–2100, 1969.
  * [3] Dena Asta and Cosma Rohilla Shalizi.  Geometric network comparison.  arXiv preprint arXiv:1411.1350, 2014.
  * [4] Nick J Baker, Riikka Kaartinen, Tomas Roslin, and Daniel B Stouffer.  Species’ roles in food webs show fidelity across a highly variable oak forest.  Ecography, 38(2):130–139, 2015.
  * [5] Alain Barrat, Marc Barthelemy, Romualdo Pastor-Satorras, and Alessandro Vespignani.  The architecture of complex weighted networks.  Proceedings of the National Academy of Sciences, 101(11):3747–3752, 2004.
  * [6] Peter S Bearman, James Moody, and Katherine Stovel.  Chains of affection: The structure of adolescent romantic and sexual networks.  American journal of sociology, 110(1):44–91, 2004.
  * [7] Sharmodeep Bhattacharyya and Peter J Bickel.  Subsampling bootstrap of count features of networks.  The Annals of Statistics, 43(6):2384–2411, 2015.
  * [8] Peter J Bickel, Aiyou Chen, and Elizaveta Levina.  The method of moments and degree distributions for network models.  The Annals of Statistics, 39(5):2280–2301, 2011.
  * [9] Peter J Bickel and David A Freedman.  Some asymptotic theory for the bootstrap.  The annals of statistics, 9(6):1196–1217, 1981.
  * [10] Patrick Billingsley.  Probability and measure.  John Wiley & Sons, 3rd edition, 1995.
  * [11] Herman Callaert and Noel Veraverbeke.  The order of the normal approximation for a studentized u-statistic.  The Annals of Statistics, pages 194–200, 1981.
  * [12] Fan Chung and Linyuan Lu.  The average distances in random graphs with given expected degrees.  Proceedings of the National Academy of Sciences, 99(25):15879–15882, 2002.
  * [13] Matthew Coulson, Robert E Gaunt, and Gesine Reinert.  Poisson approximation of subgraph counts in stochastic block models and a graphon model.  ESAIM: Probability and Statistics, 20:131–142, 2016.
  * [14] Persi Diaconis and Svante Janson.  Graph limits and exchangeable random graphs.  Rendiconti di Matematica e delle sue Applicazioni. Serie VII, 28(1):33–61, 2008.
  * [15] Joseph L Doob.  Stochastic processes, volume 7.  Wiley New York, 1953.
  * [16] Carsten F Dormann, Jochen Fründ, Nico Blüthgen, and Bernd Gruber.  Indices, graphs and null models: analyzing bipartite ecological networks.  The Open Ecology Journal, 2(1), 2009.
  * [17] B Efron.  Bootstrap methods: Another look at the jackknife.  The Annals of Statistics, 7(1):1–26, 1979.
  * [18] Frank Emmert-Streib, Matthias Dehmer, and Yongtang Shi.  Fifty years of graph matching, network alignment and network comparison.  Information sciences, 346:180–197, 2016.
  * [19] Ove Frank and David Strauss.  Markov graphs.  Journal of the American Statistical association, 81(395):832–842, 1986.
  * [20] Chao Gao and John Lafferty.  Testing network structure using relations between small subgraph probabilities.  arXiv preprint arXiv:1704.06742, 2017.
  * [21] Gérard Govaert and Mohamed Nadif.  Clustering with block mixture models.  Pattern Recognition, 36(2):463–473, 2003.
  * [22] Wassily Hoeffding.  A class of statistics with asymptotically normal distribution.  The Annals of Mathematical Statistics, 19(3):293–325, 1948.
  * [23] Wassily Hoeffding.  The strong law of large numbers for u-statistics.  Technical report, North Carolina State University. Dept. of Statistics, 1961.
  * [24] Peter D Hoff, Adrian E Raftery, and Mark S Handcock.  Latent space approaches to social network analysis.  Journal of the american Statistical association, 97(460):1090–1098, 2002.
  * [25] Olav Kallenberg.  Probabilistic symmetries and invariance principles, volume 9.  Springer, 2005.
  * [26] Nadav Kashtan, Shalev Itzkovitz, Ron Milo, and Uri Alon.  Topological generalizations of network motifs.  Physical Review E, 70(3):031909, 2004.
  * [27] Jose B Lanuza, Alfonso Allen-Perkins, and Ignasi Bartomeus.  The non-random assembly of network motifs in plant-pollinator networks.  Journal of Animal Ecology, 2023.
  * [28] Steffen Lauritzen, Alessandro Rinaldo, and Kayvan Sadeghi.  Random networks, graphical models and exchangeability.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 80(3):481–508, 2018.
  * [29] Tâm Le Minh.  U-statistics on bipartite exchangeable networks.  ESAIM: Probability and Statistics, 27:576–620, 2023.
  * [30] P-AG Maugis, SC Olhede, CE Priebe, and PJ Wolfe.  Testing for equivalence of network distribution using subgraph counts.  Journal of Computational and Graphical Statistics, 29(3):455–465, 2020.
  * [31] M. J Michalska-Smith and S. Allesina.  Telling ecological networks apart by their structure: A computational challenge.  PLoS Computational Biology, 15(6):e1007076, 2019.  Data available at github.com/mjsmith037/Classifying_Bipartite_Networks.
  * [32] Zacharie Naulet, Daniel M Roy, Ekansh Sharma, and Victor Veitch.  Bootstrap estimators for the tail-index and for the count statistics of graphex processes.  Electronic Journal of Statistics, 15:282–325, 2021.
  * [33] Mark EJ Newman.  The structure of scientific collaboration networks.  Proceedings of the national academy of sciences, 98(2):404–409, 2001.
  * [34] Krzysztof Nowicki and Tom A B Snijders.  Estimation and prediction for stochastic blockstructures.  Journal of the American statistical association, 96(455):1077–1087, 2001.
  * [35] Chiara Orsini, Marija M Dankulov, Pol Colomer-de Simón, Almerima Jamakovic, Priya Mahadevan, Amin Vahdat, Kevin E Bassler, Zoltán Toroczkai, Marián Boguná, Guido Caldarelli, et al.  Quantifying randomness in real networks.  Nature communications, 6(1):1–10, 2015.
  * [36] Sarah Ouadah, Pierre Latouche, and Stéphane Robin.  Motif-based tests for bipartite networks.  Electronic Journal of Statistics, 16(1):293–330, 2022.
  * [37] Natasa Pržulj, Derek G Corneil, and Igor Jurisica.  Modeling interactome: scale-free or geometric?  Bioinformatics, 20(18):3508–3515, 2004.
  * [38] Teresa M Przytycka.  An important connection between network motifs and parsimony models.  In Annual International Conference on Research in Computational Molecular Biology, pages 321–335. Springer, 2006.
  * [39] Gesine Reinert and Adrian Röllin.  Random subgraph counts and u-statistics: multivariate normal approximation via exchangeable pairs and embedding.  Journal of applied probability, 47(2):378–393, 2010.
  * [40] William R Schucany and Daniel M Bankson.  Small sample variance estimators for u-statistics.  Australian Journal of Statistics, 31(3):417–426, 1989.
  * [41] Pranab Kumar Sen.  On some convergence properties of ustatistics.  Calcutta Statistical Association Bulletin, 10(1-2):1–18, 1960.
  * [42] Pranab Kumar Sen.  Some invariance principles relating to jackknifing and their role in sequential analysis.  The Annals of Statistics, pages 316–329, 1977.
  * [43] Shai S Shen-Orr, Ron Milo, Shmoolik Mangan, and Uri Alon.  Network motifs in the transcriptional regulation network of escherichia coli.  Nature genetics, 31(1):64–68, 2002.
  * [44] Bernard W Silverman.  Limit theorems for dissociated random variables.  Advances in Applied Probability, 8(4):806–819, 1976.
  * [45] Benno I Simmons, Alyssa R Cirtwill, Nick J Baker, Hannah S Wauchope, Lynn V Dicks, Daniel B Stouffer, and William J Sutherland.  Motifs in bipartite ecological networks: uncovering indirect interactions.  Oikos, 128(2):154–170, 2019.
  * [46] Tom AB Snijders.  The degree variance: an index of graph heterogeneity.  Social networks, 3(3):163–174, 1981.
  * [47] Tom AB Snijders and Krzysztof Nowicki.  Estimation and prediction for stochastic blockmodels for graphs with latent block structure.  Journal of classification, 14(1):75–100, 1997.
  * [48] Daniel B Stouffer, Juan Camacho, Wenxin Jiang, and Luís A Nunes Amaral.  Evidence for the existence of a robust pattern of prey selection in food webs.  Proceedings of the Royal Society B: Biological Sciences, 274(1621):1931–1940, 2007.
  * [49] Mattia Tantardini, Francesca Ieva, Lucia Tajoli, and Carlo Piccardi.  Comparing methods for comparing networks.  Scientific reports, 9(1):1–19, 2019.
  * [50] Aad W Van der Vaart.  Asymptotic statistics, volume 3.  Cambridge university press, 2000.
  * [51] Stephen J Young and Edward R Scheinerman.  Random dot product graph models for social networks.  In Algorithms and Models for the Web-Graph: 5th International Workshop, WAW 2007, San Diego, CA, USA, December 11-12, 2007. Proceedings 5, pages 138–149. Springer, 2007.
  * [52] Liqiong Zhao, Bryce Beverlin, Theoden Netoff, and Duane Q Nykamp.  Synchronization from second order network connectivity statistics.  Frontiers in computational neuroscience, 5:28, 2011.
  * [53] Tao Zhou, Jie Ren, Matúš Medo, and Yi-Cheng Zhang.  Bipartite network projection and personal recommendation.  Physical review E, 76(4):046115, 2007.

##  Appendix A Backward martingales

Here, we present the backward martingales and their convergence theorem, which
is used to prove the convergence of some estimators. The proof of Theorem A.3
can be found in [15], Section 7, Theorem 4.2. We recall beforehand the
definition of a decreasing filtration.

######  Definition A.1.

A decreasing filtration is a decreasing sequence of $\sigma$-fields
$\mathcal{F}=(\mathcal{F}_{n})_{n\geq 1}$, i.e. such that for all $n\geq 1$,
$\mathcal{F}_{n+1}\subset\mathcal{F}_{n}$.

######  Definition A.2.

Let $\mathcal{F}=(\mathcal{F}_{n})_{n\geq 1}$ be a decreasing filtration and
$M=(M_{n})_{n\geq 1}$ a sequence of integrable random variables adapted to
$\mathcal{F}$. $(M_{n},\mathcal{F}_{n})_{n\geq 1}$ is a backward martingale if
and only if for all $n\geq 1$,
$\mathbb{E}[M_{n}\mid\mathcal{F}_{n+1}]=M_{n+1}$.

######  Theorem A.3.

Let $(M_{n},\mathcal{F}_{n})_{n\geq 1}$ be a backward martingale. Then,
$(M_{n})_{n\geq 1}$ is uniformly integrable, and, denoting
$M_{\infty}=\mathbb{E}[M_{1}\mid\mathcal{F}_{\infty}]$ where
$\mathcal{F}_{\infty}=\bigcap_{n=1}^{\infty}\mathcal{F}_{n}$, we have

$M_{n}\xrightarrow[n\rightarrow\infty]{a.s.,L_{1}}M_{\infty}.$

##  Appendix B Proofs of the results presented in Section 2

###### Proof of Lemma 2.2.

We prove this lemma by induction on $(r,c)$ the sizes of $\mathbf{i}$ and
$\mathbf{j}$. For $(r,c)=(1,0)$ and $(r,c)=(0,1)$, we have

$\mathbb{E}[p^{1,0}h(Y_{\mathbf{i},\emptyset})\mid\mathcal{A}_{\emptyset,\emptyset}]=\mathbb{E}[\psi^{1,0}h(Y_{\mathbf{i},\emptyset})]-\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]=0$

and

$\mathbb{E}[p^{0,1}h(Y_{\emptyset,\mathbf{j}})\mid\mathcal{A}_{\emptyset,\emptyset}]=\mathbb{E}[\psi^{0,1}h(Y_{\emptyset,\mathbf{i}})]-\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]=0.$

Suppose that the lemma is true for all $(0,0)<(r^{\prime},c^{\prime})<(r,c)$.
Let
$(\mathbf{i},\mathbf{j})\in\mathcal{P}_{r}(\mathbb{N})\times\mathcal{P}_{c}(\mathbb{N})$,
$\underline{\mathbf{i}}\subset\mathbf{i}$ and
$\underline{\mathbf{j}}\subset\mathbf{j}$. Denote
$\underline{r}=\mathrm{Card}(\underline{\mathbf{i}})$ and
$\underline{c}=\mathrm{Card}(\underline{\mathbf{j}})$. We can write

$\begin{split}\mathbb{E}[p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]&=\mathbb{E}[\psi^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]-\mathbb{E}[p^{\underline{r},\underline{c}}h(Y_{\underline{\mathbf{i}},\underline{\mathbf{j}}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]\\\
&\quad-\sum_{(0,0)<(r^{\prime},c^{\prime})<(r,c)}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i}),\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\\\
(\mathbf{i}^{\prime},\mathbf{j}^{\prime})\neq(\underline{\mathbf{i}},\underline{\mathbf{j}})\end{subarray}}\mathbb{E}[p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]\\\
&=\sum_{(0,0)\leq(r^{\prime},c^{\prime})<(\underline{r},\underline{c})}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\underline{\mathbf{i}}),\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\underline{\mathbf{j}})\end{subarray}}\mathbb{E}[p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]\\\
&\quad-\sum_{(0,0)<(r^{\prime},c^{\prime})<(r,c)}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i}),\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\\\
(\mathbf{i}^{\prime},\mathbf{j}^{\prime})\neq(\underline{\mathbf{i}},\underline{\mathbf{j}})\end{subarray}}\mathbb{E}[p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]\\\
&=-\sum_{(0,0)<(r^{\prime},c^{\prime})<(r,c)}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i}),\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\\\
\mathbf{i}^{\prime}\not\subset\underline{\mathbf{i}},\mathbf{j}^{\prime}\not\subset\underline{\mathbf{j}}\end{subarray}}\mathbb{E}[p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{A}_{\mathbf{i}^{\prime}\cap\underline{\mathbf{i}},\mathbf{j}^{\prime}\cap\underline{\mathbf{j}}}].\end{split}$

where we have used the fact that the
$p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})$ are
measurable by their respective
$\mathcal{A}_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}}$ so that
$\mathbb{E}[p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]=\mathbb{E}[p^{r^{\prime},c^{\prime}}h(Y_{\mathbf{i}^{\prime},\mathbf{j}^{\prime}})\mid\mathcal{A}_{\mathbf{i}^{\prime}\cap\underline{\mathbf{i}},\mathbf{j}^{\prime}\cap\underline{\mathbf{j}}}]$.

Since the last sum excludes the case
$\mathbf{i}^{\prime}=\underline{\mathbf{i}}$ and
$\mathbf{j}^{\prime}=\underline{\mathbf{j}}$, then the induction hypothesis
ensures that all the terms are equal to $0$, so
$\mathbb{E}[p^{r,c}h(Y_{\mathbf{i},\mathbf{j}})\mid\mathcal{A}_{\underline{\mathbf{i}},\underline{\mathbf{j}}}]=0$,
which concludes the proof by induction. ∎

###### Proof of Corollary 2.3.

First, we see that for some $(r_{1},c_{1})$ and $(r_{2},c_{2})$,

$\mathrm{Cov}(P^{r_{1},c_{1}}_{m,n}h_{1}(Y),P^{r_{2},c_{2}}_{m,n}h_{2}(Y))=\left[\binom{m}{r_{1}}\binom{m}{r_{2}}\binom{n}{c_{1}}\binom{n}{c_{2}}\right]^{-1}\sum_{\begin{subarray}{c}\mathbf{i}_{1}\in\mathcal{P}_{r_{1}}(\llbracket
m\rrbracket)\\\ \mathbf{j}_{1}\in\mathcal{P}_{c_{1}}(\llbracket
n\rrbracket)\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{i}_{2}\in\mathcal{P}_{r_{2}}(\llbracket
m\rrbracket)\\\ \mathbf{j}_{2}\in\mathcal{P}_{c_{2}}(\llbracket
n\rrbracket)\end{subarray}}\mathrm{Cov}(p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})).$

If $(r_{1},c_{1})\neq(r_{2},c_{2})$, then from Proposition 2.1, all the terms
$\mathrm{Cov}(p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}}))=0$,
so
$\mathrm{Cov}(P^{r_{1},c_{1}}_{m,n}h_{1}(Y),P^{r_{2},c_{2}}_{m,n}h_{2}(Y))=0$
and that is the first part of the corollary.

If $(r_{1},c_{1})=(r_{2},c_{2})=(r,c)$, then from Proposition 2.1, the terms
$\mathrm{Cov}(p^{r,c}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r,c}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}}))=0$
if $(\mathbf{i}_{1},\mathbf{j}_{1})\neq(\mathbf{i}_{2},\mathbf{j}_{2})$. Using
this fact and the exchangeability of $Y$, we have

$\begin{split}\mathrm{Cov}(P^{r,c}_{m,n}h_{1}(Y),P^{r,c}_{m,n}h_{2}(Y))&=\left[\binom{m}{r}\binom{n}{c}\right]^{-2}\sum_{\begin{subarray}{c}\mathbf{i}_{1}\in\mathcal{P}_{r}(\llbracket
m\rrbracket)\\\ \mathbf{j}_{1}\in\mathcal{P}_{c}(\llbracket
n\rrbracket)\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{i}_{2}\in\mathcal{P}_{r}(\llbracket
m\rrbracket)\\\ \mathbf{j}_{2}\in\mathcal{P}_{c}(\llbracket
n\rrbracket)\end{subarray}}\mathrm{Cov}(p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{2},\mathbf{j}_{2}}))\\\
&=\left[\binom{m}{r}\binom{n}{c}\right]^{-2}\sum_{\begin{subarray}{c}\mathbf{i}_{1}\in\mathcal{P}_{r}(\llbracket
m\rrbracket)\\\ \mathbf{j}_{1}\in\mathcal{P}_{c}(\llbracket
n\rrbracket)\end{subarray}}\mathrm{Cov}(p^{r_{1},c_{1}}h_{1}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}),p^{r_{2},c_{2}}h_{2}(Y_{\mathbf{i}_{1},\mathbf{j}_{1}}))\\\
&=\left[\binom{m}{r}\binom{n}{c}\right]^{-1}\mathrm{Cov}(p^{r,c}h_{1}(Y_{\llbracket
r\rrbracket,\llbracket c\rrbracket}),p^{r,c}h_{2}(Y_{\llbracket
r\rrbracket,\llbracket c\rrbracket})),\end{split}$

which proves the second part of the corollary. ∎

##  Appendix C Proofs of the results presented in Section 3

###### Proof of Lemma 3.2.

We only prove the first convergence result, since the second can be deduced by
analogy.

By definition,

$p^{1,0}_{(\\{i\\},\emptyset)}h=\psi^{1,0}_{(\\{i\\},\emptyset)}h-\mathbb{E}[h(Y_{(1,...,p;1,...,q)})].$

First, we see that the $p^{1,0}_{(\\{i\\},\emptyset)}h(Y)$ only depends on
$\xi_{i}$, so all the $p^{1,0}_{(\\{i\\},\emptyset)}h(Y)$ are independent.
Because $h$ is symmetric and $Y$ is RCE, they are also identically
distributed.

By the tower rule, we also have
$\mathbb{E}[\psi^{1,0}_{(\\{i\\},\emptyset)}h]=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]$
so $\mathbb{E}[p^{1,0}_{(\\{i\\},\emptyset)}h]=0$.

Finally, we have, from the Aldous-Hoover representation theorem, then Cauchy-
Schwarz’s inequality, and the fact that
$\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}]<\infty$,

$\begin{split}\mathbb{V}[p^{1,0}_{(\\{i\\},\emptyset)}h]&=v^{1,0}_{h}\\\
&=\mathbb{V}[\mathbb{E}[h(Y_{(1,...,p;1,...,q)})\mid\xi_{i}]]\\\
&=\mathrm{Cov}(h(Y_{(1,...,p;1,...,q)}),h(Y_{(1,p+1,...,2p-1;1,...,q)}))\\\
&<\mathbb{V}[h(Y_{(1,...,p;1,...,q)})]\\\ &<\infty.\end{split}$

The classical CLT gives the desired result. ∎

###### Proof of Lemma 3.3.

For all $\mathbf{i}\in\mathcal{P}_{r}(\llbracket m\rrbracket)$ and
$\mathbf{j}\in\mathcal{P}_{c}(\llbracket n\rrbracket)$ and
$(r,c)\in\mathbb{N}^{2}$, we have
$\mathbb{E}[\psi^{r,c}_{(\mathbf{i},\mathbf{j})}h]=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]$.

By recursion, we have $\mathbb{E}[A_{N}]=\mathbb{E}[P^{r,c}_{N}h]=0$ for all
$(r,c)>(0,0)$ since

$\begin{split}\mathbb{E}[P^{r,c}_{N}h]&=\mathbb{E}[p^{r,c}_{(\mathbf{i},\mathbf{j})}h]\\\
&=\mathbb{E}[\psi^{r,c}_{(\mathbf{i},\mathbf{j})}h]-\sum_{(0,0)\leq(r^{\prime},c^{\prime})<(r,c)}\sum_{\begin{subarray}{c}\mathbf{i}^{\prime}\in\mathcal{P}_{r^{\prime}}(\mathbf{i})\\\
\mathbf{j}^{\prime}\in\mathcal{P}_{c^{\prime}}(\mathbf{j})\end{subarray}}\mathbb{E}[p^{r^{\prime},c^{\prime}}_{(\mathbf{i}^{\prime},\mathbf{j}^{\prime})}h].\end{split}$

Now, Corollary 2.4 imply that
$\mathrm{Cov}(P^{r,c}_{N}h,P^{r^{\prime},c^{\prime}}_{N}h)=0$ unless
$(r,c)=(r^{\prime},c^{\prime})$. So

$\begin{split}\mathbb{V}[A_{N}]&=N\sum_{\begin{subarray}{c}(0,0)<(r,c)\leq(p,q)\\\
(r,c)\neq(1,0)\neq(0,1)\end{subarray}}\sum_{\begin{subarray}{c}(0,0)<(r^{\prime},c^{\prime})\leq(p,q)\\\
(r^{\prime},c^{\prime})\neq(1,0)\neq(0,1)\end{subarray}}\binom{p}{r}\binom{p}{r^{\prime}}\binom{q}{c}\binom{q}{c^{\prime}}\mathrm{Cov}(P^{r,c}_{N}h,P^{r^{\prime},c^{\prime}}_{N}h)\\\
&=N\sum_{\begin{subarray}{c}(0,0)<(r,c)\leq(p,q)\\\
(r,c)\neq(1,0)\neq(0,1)\end{subarray}}\binom{p}{r}^{2}\binom{q}{c}^{2}\mathbb{V}[P^{r,c}_{N}h]\\\
&=N\sum_{\begin{subarray}{c}(0,0)<(r,c)\leq(p,q)\\\
(r,c)\neq(1,0)\neq(0,1)\end{subarray}}\binom{p}{r}^{2}\binom{q}{c}^{2}\binom{m_{N}}{r}^{-1}\binom{n_{N}}{c}^{-1}\mathbb{V}[p^{r,c}_{(\llbracket
r\rrbracket,\llbracket c\rrbracket)}h].\end{split}$ (18)

From equation 5, we see that $h(Y_{(1,...p,1,...,q)})$ is a linear combination
of all the $(p^{r,c}_{(\mathbf{i},\mathbf{j})}h)_{\begin{subarray}{c}0\leq
r\leq p,0\leq c\leq q\\\ \mathbf{i}\in\mathcal{P}_{r}(\llbracket
p\rrbracket),\mathbf{j}\in\mathcal{P}_{c}(\llbracket
q\rrbracket)\end{subarray}}$. So
$\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}]<\infty$ ensures that
$\mathbb{V}[p^{r,c}_{(\llbracket r\rrbracket,\llbracket
c\rrbracket)}h]<\infty$. Therefore
$\mathbb{V}[A_{N}]=O(Nm_{N}^{-r}n_{N}^{-c})=O(N^{1-r-c})$.

Thus, we deduce from Markov’s inequality that
$A_{N}\xrightarrow[N\rightarrow\infty]{\mathbb{P}}0$. ∎

######  Lemma C.1.

  1. 1.

$U^{\tilde{h}}_{N}=U^{h}_{N}$,

  2. 2.

$\mathbb{E}[\tilde{h}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})]=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]$,

  3. 3.

$\mathbb{E}[\tilde{h}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]\leq\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}]$.

###### Proof of Lemma C.1.

The two first properties are straightforward. The third property stems from

$\begin{split}\mathbb{E}[\tilde{h}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]&=\left[\binom{p^{\prime}}{p}\binom{q^{\prime}}{q}\right]^{-2}\sum_{\begin{subarray}{c}\mathbf{i}_{1}\subset\\{i_{1},...,i_{p^{\prime}}\\}\\\
\mathbf{j}_{1}\subset\\{j_{1},...,j_{q^{\prime}}\\}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{i}_{2}\subset\\{i_{1},...,i_{p^{\prime}}\\}\\\
\mathbf{j}_{2}\subset\\{j_{1},...,j_{q^{\prime}}\\}\end{subarray}}\mathbb{E}[h(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})h(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})]\\\
&=\left[\binom{p^{\prime}}{p}\binom{q^{\prime}}{q}\right]^{-2}\sum_{\begin{subarray}{c}\mathbf{i}_{1}\subset\\{i_{1},...,i_{p^{\prime}}\\}\\\
\mathbf{j}_{1}\subset\\{j_{1},...,j_{q^{\prime}}\\}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{i}_{2}\subset\\{i_{1},...,i_{p^{\prime}}\\}\\\
\mathbf{j}_{2}\subset\\{j_{1},...,j_{q^{\prime}}\\}\end{subarray}}\left(\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}]-\frac{1}{2}\mathbb{E}[(h(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})-h(Y_{\mathbf{i}_{2},\mathbf{j}_{2}}))^{2}]\right)\\\
&\leq\left[\binom{p^{\prime}}{p}\binom{q^{\prime}}{q}\right]^{-2}\sum_{\begin{subarray}{c}\mathbf{i}_{1}\subset\\{i_{1},...,i_{p^{\prime}}\\}\\\
\mathbf{j}_{1}\subset\\{j_{1},...,j_{q^{\prime}}\\}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{i}_{2}\subset\\{i_{1},...,i_{p^{\prime}}\\}\\\
\mathbf{j}_{2}\subset\\{j_{1},...,j_{q^{\prime}}\\}\end{subarray}}\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}]\\\
&\leq\mathbb{E}[h(Y_{(1,...,p;1,...,q)})^{2}].\end{split}$

∎

###### Proof of Corollary 3.5.

First, we show that
$\lim_{N\rightarrow+\infty}N\mathrm{Cov}(U^{h_{k}}_{N},U^{h_{\ell}}_{N})$

Let $(Z^{h_{k}})_{1\leq k\leq D}$ be a vector of random variables following a
centered multivariate Gaussian distribution with covariance matrix
$\Sigma^{h_{1},...,h_{D}}$ defined in the theorem. Then
$Z^{h_{k}}\sim\mathcal{N}(0,V^{h_{k}})$ for all $1\leq k\leq D$ and
$\mathrm{Cov}(Z^{h_{k}},Z^{h_{\ell}})=C^{h_{k},h_{\ell}}$ for all $1\leq k\leq
D$ and $1\leq\ell\leq D$.

Denote $p^{\prime}:=\max_{k}(p_{k})$ and $q^{\prime}:=\max_{k}(q_{k})$. For
some $t=(t_{1},t_{2},...,t_{n})\in\mathbb{R}^{n}$, we set
$\tilde{h}_{t}:=t_{1}\tilde{h}_{1}+t_{2}\tilde{h}_{2}+...+t_{n}\tilde{h}_{D}$.
$\tilde{h}_{t}$ is a kernel function of size $p^{\prime}\times q^{\prime}$.

First, assume that $t\neq(0,...,0)$. Then by hypothesis,
$\tilde{h}_{t}\not\equiv 0$, therefore Proposition C.1 implies that
$\sum_{k=1}^{D}t_{k}U^{h_{k}}_{N}=\sum_{k=1}^{D}t_{k}U^{\tilde{h}_{k}}_{N}=U^{\tilde{h}_{t}}_{N}$,
the $U$-statistic with quadruplet kernel $\tilde{h}_{t}$ (of size
$p^{\prime}\times q^{\prime}$). Using Cauchy-Schwarz inequality and the fact
that from Proposition C.1,
$\mathbb{E}[\tilde{h}_{k}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]\leq\mathbb{E}[h_{k}(Y_{(1,...,p;1,...,q)})^{2}]<\infty$
for all $1\leq k\leq D$, we have furthermore

$\begin{split}\mathbb{E}[\tilde{h}_{t}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]=&~{}\sum_{k=1}^{n}t_{k}^{2}\mathbb{E}[\tilde{h}_{k}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]+2\sum_{1\leq
k\neq\ell\leq
D}t_{k}t_{\ell}\mathbb{E}[\tilde{h}_{k}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})\tilde{h}_{\ell}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})],\\\
\leq&~{}\sum_{k=1}^{n}t_{k}^{2}\mathbb{E}[\tilde{h}_{k}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]+2\sum_{1\leq
k\neq\ell\leq
D}t_{k}t_{\ell}\sqrt{\mathbb{E}[\tilde{h}_{k}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]\mathbb{E}[\tilde{h}_{\ell}(Y_{(1,...,p^{\prime};1,...,q^{\prime})})^{2}]},\\\
<&~{}\infty.\end{split}$

Therefore, Theorem 3.1 also applies for $U^{\tilde{h}_{t}}_{N}$ and
$\sqrt{N}(U^{\tilde{h}_{t}}_{N}-U^{\tilde{h}_{t}}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\mathcal{N}(0,V^{\tilde{h}_{t}})$,
where :

  * •

$U^{\tilde{h}_{t}}_{\infty}=\sum_{k=1}^{D}t_{k}U^{h_{k}}_{\infty}$,

  * •

$V^{\tilde{h}_{t}}=t^{T}\Sigma^{h_{1},...,h_{D}}t$.

The second point comes from the fact that
$V^{\tilde{h}_{t}}=\lim_{N\rightarrow+\infty}N\sum_{k=1}^{D}\sum_{\ell=1}^{D}t_{k}t_{\ell}\mathrm{Cov}(U^{h_{k}}_{N},U^{h_{\ell}}_{N})$
and by Corollary 2.4,
$\lim_{N\rightarrow+\infty}N\mathrm{Cov}(U^{h_{k}}_{N},U^{h_{\ell}}_{N})=\frac{p^{2}}{c}c^{1,0}_{h_{k},h_{\ell}}+\frac{q^{2}}{1-c}c^{0,1}_{h_{k},h_{\ell}}=\Sigma^{h_{1},...,h_{D}}_{k\ell}$.
Therefore, we can conclude that
$\sqrt{N}\sum_{k=1}^{D}t_{k}(U^{h_{k}}_{N}-U^{h_{k}}_{\infty})=\sqrt{N}(U^{\tilde{h}_{t}}_{N}-U^{\tilde{h}_{t}}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\sum_{k=1}^{D}t_{k}Z^{h_{k}}$.

Now assume that $t=(0,...,0)$. Then $\tilde{h}_{t}\equiv 0$ so
$U^{\tilde{h}_{t}}_{N}=0=\sum_{k=1}^{D}t_{k}Z^{h_{k}}$. Therefore,
$\sqrt{N}\sum_{k=1}^{D}t_{k}(U^{h_{k}}_{N}-U^{h_{k}}_{\infty})=\sqrt{N}(U^{\tilde{h}_{t}}_{N}-U^{\tilde{h}_{t}}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\sum_{k=1}^{D}t_{k}Z^{h_{k}}$
is still true.

We have proven that
$\sqrt{N}\sum_{k=1}^{D}t_{k}(U^{h_{k}}_{N}-U^{h_{k}}_{\infty})\xrightarrow[N\rightarrow\infty]{\mathcal{D}}\sum_{k=1}^{D}t_{k}Z^{h_{k}}$
for all $t\in\mathbb{R}^{n}$, so we can finally apply the Cramér-Wold theorem
(Theorem 29.4 of [10]) which states that
$\sqrt{N}\left(U^{h_{k}}_{N}-U^{h_{k}}_{\infty}\right)_{1\leq k\leq D}$
converges jointy in distribution to $(Z^{h_{k}})_{1\leq k\leq D}$, which is a
centered multivariate Gaussian with covariance matrix
$\Sigma^{h_{1},...,h_{D}}$. ∎

##  Appendix D Proofs of the results presented in Section 4

######  Lemma D.1.

Let $k\in\mathbb{N}$. Then as $n\rightarrow\infty$,

$\binom{n}{k}=\frac{n^{k}}{k!}+O(n^{k-1})$

###### Proof of Lemma 4.1.

The sum $T^{p,q}_{N}(I_{K},J_{K})$ defined by equation (6) is a sum over
$T_{N}$ expectation terms where

$T_{N}:=\mathrm{Card}(\mathcal{T}^{p,q}_{N,(I_{K},J_{K})})=\prod_{k=1}^{K}\mathrm{Card}(\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}}_{k},\mathbf{\underline{j}}_{k})})=\prod_{k=1}^{K}\binom{m_{N}-\underline{p}_{k}}{p-\underline{p}_{k}}\binom{n_{N}-\underline{q}_{k}}{q-\underline{q}_{k}}.$
(19)

These expectation terms
$\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]$
only depend on the number of times elements appear in pairwise intersections
between the $(\mathbf{i}_{k},\mathbf{j}_{k})$, $1\leq k\leq K$. In particular,
$\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]=\alpha(I_{K},J_{K})$
when all the $\mathbf{i}_{k}$ contains the elements of
$\mathbf{\underline{i}}_{k}$ and all the other elements do not appear in any
other $\mathbf{\underline{i}}_{k^{\prime}}$ (see equation 7). Denote $A_{N}$
the number of terms of $T^{p,q}_{N}(I_{K},J_{K})$ where
$\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]=\alpha(I_{K},J_{K})$,
then

$T^{p,q}_{N}(I_{K},J_{K})=\frac{A_{N}}{T_{N}}\alpha(I_{K},J_{K})+\frac{1}{T_{N}}\sum_{\begin{subarray}{c}\mathcal{T}^{p,q}_{N,(I_{K},J_{K})}\\\
\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]\neq\alpha(I_{K},J_{K})\end{subarray}}\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}X_{\mathbf{i}_{2},\mathbf{j}_{2}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}].$
(20)

Using Jensen’s inequality, Hölder’s inequality and the exchangeability of the
submatrices, we have for all the expectation terms,

$0\leq\lvert\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]\rvert\leq\mathbb{E}[\lvert
X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}\rvert]\leq\prod_{k=1}^{K}\mathbb{E}[\lvert
X_{\mathbf{i}_{k},\mathbf{j}_{k}}\rvert^{K}]^{\frac{1}{K}}=\mathbb{E}[\lvert
X_{(1,...,p;1,...,q)}\rvert^{K}].$

In particular, this holds for all $(T_{N}-A_{N})$ terms of the remaining sum
in Equation (20) when
$\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]\neq\alpha(I_{K},J_{K})$,
so

$0\leq\lvert
T^{p,q}_{N}(I_{K},J_{K})-\frac{A_{N}}{T_{N}}\alpha(I_{K},J_{K})\rvert\leq\left(1-\frac{A_{N}}{T_{N}}\right)\mathbb{E}[\lvert
X_{(1,...,p;1,...,q)}\rvert^{K}].$ (21)

We need to calculate $T_{N}$ and $A_{N}$ to use this inequality and conclude
the proof. We organize the rest of the proof in 3 parts.

  1. 1.

We find an expression for $T_{N}$. Since for all $1\leq k\leq K$, from Lemma
D.1,

$\binom{m_{N}-\underline{p}_{k}}{p-\underline{p}_{k}}\binom{n_{N}-\underline{q}_{k}}{q-\underline{q}_{k}}=\frac{m_{N}^{p-\underline{p}_{k}}n_{N}^{q-\underline{q}_{k}}}{(p-\underline{p}_{k})!(q-\underline{q}_{k})!}\left(1+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\right),$

we have from Equation (19)

$T_{N}=\left(\prod_{k=1}^{K}(p-\underline{p}_{k})!(q-\underline{q}_{k})!\right)^{-1}m_{N}^{Kp-\underline{P}}n_{N}^{Kq-\underline{Q}}\left(1+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\right).$
(22)

  2. 2.

We find an expression for $A_{N}$. The summation over
$(\mathbf{i}_{k},\mathbf{j}_{k})\in\mathcal{S}^{p,q}_{N,(\mathbf{\underline{i}}_{k},\mathbf{\underline{j}}_{k})}$
is in fact a sum over the $p-\underline{p}_{k}$ elements of $\mathbf{i}_{k}$
and $q-\underline{q}_{k}$ elements of $\mathbf{j}_{k}$ that are not restricted
by $\mathbf{\underline{i}}_{k}$ and $\mathbf{\underline{j}}_{k}$.

     * •

We can pick the first $\mathbf{i}_{1}$ choosing the $p-\underline{p}_{1}$
unrestricted indices among the $m_{N}-\bar{p}$ values of $\\{1,...,m_{N}\\}$
excluding $\cup_{k=1}^{K}\mathbf{\underline{i}}_{k}$. The same follows for the
pick of $\mathbf{j}_{1}$, so there are
$\binom{m_{N}-\bar{p}}{p-\underline{p}_{1}}\binom{n_{N}-\bar{q}}{q-\underline{q}_{1}}$
possible picks for $(\mathbf{i}_{1},\mathbf{j}_{1})$.

     * •

The pick of $\mathbf{i}_{2}$ consists in choosing the $p-\underline{p}_{2}$
unrestricted indices among the $m_{N}-\bar{p}-(p-\underline{p}_{1})$ values of
$\\{1,...,m_{N}\\}$ excluding $\cup_{k=1}^{K}\mathbf{\underline{i}}_{k}$ and
the elements already taken by $\mathbf{i}_{1}$. We deduce that there are
$\binom{m_{N}-\bar{p}-(p-\underline{p}_{1})}{p-\underline{p}_{2}}\binom{n_{N}-\bar{q}-(q-\underline{q}_{1})}{q-\underline{q}_{2}}$
possible picks for $(\mathbf{i}_{2},\mathbf{j}_{2})$.

     * •

Iteratively, for all $1\leq k\leq K$, we find that there are
$\binom{m_{N}-\bar{p}-\sum_{k^{\prime}=1}^{k-1}(p-\underline{p}_{k^{\prime}})}{p-\underline{p}_{k}}\binom{n_{N}-\bar{q}-\sum_{k^{\prime}=1}^{k-1}(q-\underline{q}_{k^{\prime}})}{q-\underline{q}_{k}}$
possible picks for $(\mathbf{i}_{k},\mathbf{j}_{k})$.

We deduce that the number of possible picks for all the
$(\mathbf{i}_{k},\mathbf{j}_{k})$, $1\leq k\leq K$ so that
$\mathbb{E}[X_{\mathbf{i}_{1},\mathbf{j}_{1}}...X_{\mathbf{i}_{K},\mathbf{j}_{K}}]=\alpha(I_{K},J_{K})$
is

$A_{N}=\prod_{k=1}^{K}\binom{m_{N}-\bar{p}-\sum_{k^{\prime}=1}^{k-1}(p-\underline{p}_{k^{\prime}})}{p-\underline{p}_{k}}\binom{n_{N}-\bar{q}-\sum_{k^{\prime}=1}^{k-1}(q-\underline{q}_{k^{\prime}})}{q-\underline{q}_{k}}.$

But we see that from Lemma D.1, for $1\leq k\leq K$,

$\binom{m_{N}-\bar{p}-\sum_{k^{\prime}=1}^{k-1}(p-\underline{p}_{k^{\prime}})}{p-\underline{p}_{k}}\binom{n_{N}-\bar{q}-\sum_{k^{\prime}=1}^{k-1}(q-\underline{q}_{k^{\prime}})}{q-\underline{q}_{k}}=\frac{m_{N}^{p-\underline{p}_{k}}n_{N}^{q-\underline{q}_{k}}}{(p-\underline{p}_{k})!(q-\underline{q}_{k})!}\left(1+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\right).$

So we have

$A_{N}=\left(\prod_{k=1}^{K}(p-\underline{p}_{k})!(q-\underline{q}_{k})!\right)^{-1}m_{N}^{Kp-\underline{P}}n_{N}^{Kq-\underline{Q}}\left(1+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\right).$

  3. 3.

Now with the expressions of $T_{N}$ and $A_{N}$, we can deduce that

$T_{N}-A_{N}=O\left(m_{N}^{Kp-\underline{P}}n_{N}^{Kq-\underline{Q}}(m_{N}^{-1}+n_{N}^{-1})\right),$

so

$1-\frac{A_{N}}{T_{N}}=\frac{O\left(m_{N}^{Kp-\underline{P}}n_{N}^{Kq-\underline{Q}}(m_{N}^{-1}+n_{N}^{-1})\right)}{\left(\prod_{k=1}^{K}(p-\underline{p}_{k})!(q-\underline{q}_{k})!\right)^{-1}m_{N}^{Kp-\underline{P}}n_{N}^{Kq-\underline{Q}}\left(1+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\right)}=O\left(m_{N}^{-1}+n_{N}^{-1}\right).$

Therefore, we can finally conclude using Equation (21),

$\begin{split}T^{p,q}_{N}(I_{K},J_{K})&=\left(1+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\right)\alpha(I_{K},J_{K})+O\left(m_{N}^{-1}+n_{N}^{-1}\right)\\\
&=\alpha(I_{K},J_{K})+O\left(m_{N}^{-1}+n_{N}^{-1}\right).\end{split}$

∎

###### Proof of Proposition 4.5.

In this proof, for the estimator defined by equation (8), we write
$\widehat{\mu}^{(i)}_{N}$ instead of $\widehat{\mu}^{h,(i)}_{N}$ to simplify
the notation without ambiguity.

Notice that

$\begin{split}\mathbb{E}[\widehat{v}^{h;1,0}_{N}]&=\frac{1}{2}\mathbb{E}[(\widehat{\mu}^{(1)}_{N}-\widehat{\mu}^{(2)}_{N})^{2}]\\\
&=\mathbb{E}[(\widehat{\mu}^{(1)}_{N})^{2}]-\mathbb{E}[\widehat{\mu}^{(1)}_{N}\widehat{\mu}^{(2)}_{N}].\end{split}$
(23)

  * •

First, we calculate $\mathbb{E}[(\widehat{\mu}^{(1)}_{N})^{2}]$ :

$\begin{split}\mathbb{E}[\widehat{\mu}^{(1)}_{N}(Y)^{2}]&=\binom{m_{N}-1}{p-1}^{-2}\binom{n_{N}}{q}^{-2}\sum_{(\mathbf{i}_{1},\mathbf{j}_{1})\in\mathcal{S}^{p,q}_{N,(\\{1\\},\emptyset)}}\sum_{(\mathbf{i}_{2},\mathbf{j}_{2})\in\mathcal{S}^{p,q}_{N,(\\{1\\},\emptyset)}}\mathbb{E}[h(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})h(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})]\\\
&=T^{p,q}_{N}(I_{2},J_{2})\end{split}$

where $I_{2}=(\\{1\\},\\{1\\})$ and $J_{2}=(\emptyset,\emptyset)$. Applying
Lemma 4.1, with
$\underline{P}=\mathrm{Card}(\\{1\\})+\mathrm{Card}(\\{1\\})=2$,
$\underline{Q}=\mathrm{Card}(\emptyset)+\mathrm{Card}(\emptyset)=0$ and

$\begin{split}\alpha(I_{2},J_{2})&=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})h(Y_{(1,p+1,...,2p-1;q+1,...,2q)})]\\\
&=\mathbb{E}[\mathbb{E}[h(Y_{(1,...,p;1,...,q)})h(Y_{(1,p+1,...,2p-1;q+1,...,2q)})\mid\xi_{1}]]\\\
&=\mathbb{E}[\mathbb{E}[h(Y_{(1,...,p;1,...,q)})\mid\xi_{1}]^{2}]\\\
&=\mathbb{E}[\psi^{1,0}_{(\\{1\\},\emptyset)}h^{2}].\end{split}$

we find

$\begin{split}\mathbb{E}[(\widehat{\mu}^{(1)}_{N})^{2}]&=\mathbb{E}[\psi^{1,0}_{(\\{1\\},\emptyset)}h^{2}]+O\left(N^{-1}\right).\end{split}$
(24)

  * •

Next, we calculate
$\mathbb{E}[\widehat{\mu}^{(1)}_{N}\widehat{\mu}^{(2)}_{N}]$ :

$\begin{split}\mathbb{E}[\widehat{\mu}^{(1)}_{N}\widehat{\mu}^{(2)}_{N}]&=\binom{m_{N}-1}{p-1}^{-2}\binom{n_{N}}{q}^{-2}\sum_{(\mathbf{i}_{1},\mathbf{j}_{1})\in\mathcal{S}^{p,q}_{N,(\\{1\\},\emptyset)}}\sum_{(\mathbf{i}_{2},\mathbf{j}_{2})\in\mathcal{S}^{p,q}_{N,(\\{2\\},\emptyset)}}\mathbb{E}[h(Y_{\mathbf{i}_{1},\mathbf{j}_{1}})h(Y_{\mathbf{i}_{2},\mathbf{j}_{2}})]\\\
&=T^{p,q}_{N}(I_{2}^{\prime},J_{2}^{\prime})\\\ \end{split}$

where $I_{2}^{\prime}=(\\{1\\},\\{2\\})$ and
$J_{2}^{\prime}=(\emptyset,\emptyset)$. Applying Lemma 4.1 with

$\begin{split}\alpha(I_{2}^{\prime},J_{2}^{\prime})&=\mathbb{E}[h(Y_{(1,3,...,p+1;1,...,q)})h(Y_{(2,p+2,...,2p;q+1,...,2q)})]\\\
&=\mathbb{E}[h(Y_{(1,...,p;1,...,q)})]^{2}\\\
&=\mathbb{E}[\psi^{1,0}_{(\\{1\\},\emptyset)}h]^{2},\end{split}$

we find

$\begin{split}\mathbb{E}[\widehat{\mu}^{(1)}_{N}\widehat{\mu}^{(2)}_{N}]&=\mathbb{E}[\psi^{1,0}_{(\\{1\\},\emptyset)}h]^{2}+O\left(N^{-1}\right).\end{split}$
(25)

Finally, we can combine Equations (23), (24) and (25) to obtain

$\mathbb{E}[\widehat{v}^{h;1,0}_{N}]=\mathbb{E}[\psi^{1,0}_{(\\{1\\},\emptyset)}h^{2}]-\mathbb{E}[\psi^{1,0}_{(\\{1\\},\emptyset)}h]^{2}+O\left(N^{-1}\right)=v^{1,0}_{h}+O\left(N^{-1}\right),$

which proves the theorem. ∎

###### Proof of Proposition 4.6.

In this proof, we write $\widehat{\mu}^{(i)}_{N}$ instead of