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

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

► To cite this version:

Tâm Le Minh, Sophie Donnet, François Massol, Stéphane Robin. Hoeffding-type decomposition for U-statistics on bipartite networks. Electronic Journal of Statistics, 2025, 19 (1), pp.2829-2875. 10.1214/25-EJS2402 . hal-05215215

HAL Id: hal-05215215

https://hal.science/hal-05215215v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Hoeffding-type decomposition for *U*-statistics on bipartite networks

Tâm Le Minh<fnref n="1" />, Sophie Donnet<fnref n="1" />, François Massol<fnref n="2" /> and Stéphane Robin<fnref n="3" />

**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 based on the Aldous-Hoover-Kallenberg representation of row-column exchangeable matrices. 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.

**MSC2020 subject classifications:** Primary 62G20, 60F05; secondary 62G10.

**Keywords and phrases:** Bipartite networks, *U*-statistics, variance estimation, Hoeffding decomposition, row-column exchangeability, Central Limit Theorem.

Received February 2025.

## 1. Introduction

Networks are popular objects to represent a set of interacting entities. The last decades have witnessed an explosion in the number of network 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 into 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 statistics considered is the *U*-statistics which, in the context of network analysis, have complex dependencies.

<footnote n="1">¹Université Paris-Saclay, AgroParisTech, INRAE, UMR MIA Paris-Saclay, 91120, Palaiseau, France, e-mail: tam.le-minh@inria.fr; sophie.donnet@inrae.fr</footnote>

<footnote n="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, e-mail: francois.massol@univ-lille.fr</footnote>

<footnote n="3">³Sorbonne Université, CNRS, Laboratoire de Probabilités, Statistique et Modélisation, e-mail: stephane.robin@sorbonne-universite.fr</footnote>
**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 [67], papers and scientists in authorship networks [43], or plants and pollinators in ecological interaction networks [21]. 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 [60], the random dot product graph model [65] or the latent space model [29] 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 $N$,

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

where $\Phi Y := (Y_{\sigma_1(i)\sigma_2(j)})_{i \ge 1, j \ge 1}$. Many exchangeable random graph models also have a dissociatedness property, i.e. their adjacency matrices are also dissociated [57, 36]. An RCE matrix is said to be dissociated if for all $m$ and $n$, $(Y_{ij})_{1 \le i \le m, 1 \le j \le n}$ is independent of $(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 \ge 1}$ numbered with a unique index a U-statistic $U_n^h(Y)$ of order $n$ and kernel function $h$ is defined as the following average

$$ U_n^h(Y) = \left(\binom{n}{k}\right)^{-1} \sum_{1 \le i_1 < i_2 < \dots < i_k \le n} h(Y_{i_1}, Y_{i_2}, \dots, Y_{i_k}), $$

where $h: \mathbb{R}^k \to \mathbb{R}$ is a symmetric function referred to as the kernel. Denote $\llbracket n \rrbracket := \{1, \dots, 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, \dots, i_k\} \in \mathcal{P}_k(\llbracket n \rrbracket)$, then by symmetry of $h$, $h(Y_{i_1}, \dots, Y_{i_k})$ does not depend on the order of the elements of $\mathbf{i}$. Therefore, we will denote $h(\mathbf{Y}_i) := h(Y_{i_1}, \dots, Y_{i_k})$. Finally, the U-statistic $U_n^h(Y)$ can be formulated as follows:

$$ U_n^h(Y) = \left(\binom{n}{k}\right)^{-1} \sum_{\mathbf{i} \in \mathcal{P}_k(\llbracket n \rrbracket)} h(\mathbf{Y}_i). $$

When $Y$ is an exchangeable sequence, $h(\mathbf{Y}_i)$ has the same distribution for all $\mathbf{i} \in \mathcal{P}_k(\llbracket n \rrbracket)$, therefore $U_n^h(Y)$ is an unbiased estimate of $h(\mathbf{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 [27] 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 [28]. For $1 \le c \le k$, define the function $\psi^c h: \mathbb{R}^c \to \mathbb{R}$ as

$$ \psi^c h : (y_1, \dots, y_c) \mapsto \mathbb{E}[h(Y_1, \dots, Y_k) | Y_1 = y_1, \dots, 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}, \dots, 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'=0}^{c-1} \sum_{\mathbf{i}' \in \mathcal{P}_{c'}(\mathbf{i})} p^{c'} h(Y_{\mathbf{i}'}). $$

for all subsets $\mathbf{i} \in \mathcal{P}_c(\llbracket n \rrbracket)$, for all $1 \le c \le 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 \le c \le k} \sum_{\mathbf{i}' \in \mathcal{P}_c(\mathbf{i})} p^c h(Y_{\mathbf{i}'}). $$

The U-statistic $U_n^h$ can be written as

$$ U_n^h(Y) = \sum_{c=0}^{k} \binom{k}{c} P_n^c h(Y), $$

where for $1 \le c \le k$, $P_n^c h(Y) = \binom{n}{c}^{-1} \sum_{\mathbf{i} \in \mathcal{P}_c(\llbracket n \rrbracket)} 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_n^c h$ are also orthogonal. The leading terms of this decomposition have been used by [27] to prove the asymptotic normality of U-statistics. The decomposition also yields a decomposition of the variance of U-statistics.

This framework naturally extends to a network U-statistics setting, where the kernel $h$ is defined on submatrices 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}) \to \mathbb{R}$ be a function defined on $p \times q$ matrices, $1 \le p \le m$, $1 \le q \le 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)}, \dots, i_{\sigma_1(p)}; j_{\sigma_2(1)}, \dots, j_{\sigma_2(q)})}) = h(Y_{(i_1, i_2, \dots, i_p; j_1, j_2, \dots, j_q)}), \quad (1) $$

where $Y_{(i_1, \dots, i_p; j_1, \dots, j_q)}$ is the $p \times q$ submatrix consisting of the rows and columns of $Y$ indexed by $i_1, \dots, i_p$ and $j_1, \dots, j_q$ respectively. Therefore, since the order of the elements of $\mathbf{i} = \{i_1, \dots, i_p\}$ and $\mathbf{j} = \{j_1, \dots, j_q\}$ does not matter, we can denote $h(Y_{\mathbf{i},\mathbf{j}}) := h(Y_{(i_1, i_2, \dots, i_p; j_1, j_2, \dots, j_q)})$. Then the associated U-statistic is

$$ U_{m,n}^h(Y) = \binom{m}{p}^{-1} \binom{n}{q}^{-1} \sum_{\substack{\mathbf{i} \in \mathcal{P}_p(\llbracket m \rrbracket) \\ \mathbf{j} \in \mathcal{P}_q(\llbracket n \rrbracket)}} h(Y_{\mathbf{i},\mathbf{j}}). \quad (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}) \to \mathbb{R}$ is not symmetric, then $h: \mathcal{M}_{p,q}(\mathbb{R}) \to \mathbb{R}$ defined by

$$ h(Y_{(i_1, i_2, \dots, i_p; j_1, j_2, \dots, 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)}, \dots, i_{\sigma_1(p)}; j_{\sigma_2(1)}, \dots, j_{\sigma_2(q)})}) \quad (3) $$
verifies Equation (1) and leads to the same U-statistic $(U_{m,n}^{h0}(Y) = U_{m,n}^{h}(Y))$.

Formally, the network U-statistic defined by (2) is a weighted U-statistic [55], and more precisely, an incomplete U-statistic [13], where the weights are binary and determined by a sampling design. Flattening the matrix $Y$ into a vector in $\mathbb{R}^{m \times n}$, and interpreting $h$ as a function of $p \times q$ variables, clarifies this perspective. The asymptotic properties of weighted U-statistics [14, 30, 45] have been studied extensively, although a general distributional characterization has been derived only for the case of linear kernels [8]. Moreover, most existing results rely on the assumption of i.i.d. observations. In contrast, the RCE structure of $Y$ implies a different dependency pattern and suggests a natural sampling design that respects both row and column symmetries.

The asymptotic behavior of network U-statistics fundamentally depends on the dependency structure of the data $Y$. Many studies have considered the case where the matrix $Y$ represents an exchangeable network, more precisely through the investigation of subgraph counts [10, 34, 47, 54, 7]. Some have investigated their asymptotic behavior through the theory of generalized U-statistics [31]. However, most of them consider *unipartite* networks. These results do not apply to *bipartite* networks, since they have two disjoint sets of nodes which can be of different sizes. [64] considered bipartite graphs, but the U-statistics they have defined are averages over only the pairs of present edges not sharing any node. They do not make assumptions on the distribution of networks, so they do not derive the limit distribution of these U-statistics. [18] investigated averages of separately exchangeable arrays, which corresponds to the network U-statistic $U_{m,n}^h$ on $Y$, where $Y$ is an RCE matrix but $h$ is a function of only one variable. Therefore, their results do not directly apply. More closely related, [37] used a martingale approach to obtain a weak convergence result for $U_{m,n}^h$ 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 a Hoeffding decomposition-based approach to analyze $U_{m,n}^h$. In our context, the Hoeffding decomposition yields a simpler variance decomposition than for classic generalized U-statistics. This strategy also has the advantage of providing a method to estimate the asymptotic variance of network U-statistics. Indeed, estimating this variance is a required condition to perform practical inference tasks, such as hypothesis testing. However, it remains a complex problem that has been tackled with various methods in the literature.

**Variance estimation of U-statistics** The standard error of U-statistics of i.i.d. random variables is most often computed using resampling techniques such as the jackknife [2] and bootstrap [22, 11] estimators of variance. [52, 53, 15, 51] suggested various estimators of the asymptotic variance. However, all these estimators are biased for both the variance and the asymptotic variances. [15, 51] 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. For U-statistics on RCE matrices, no generic variance estimator has been proposed. The bootstrap procedures considered in [18, 40] are proven to be consistent for the sample average of an RCE matrix, but have not been extended to U-statistics. Previous studies have estimated the variance by identifying estimable quantities in the problem-specific analytical expressions [37, 38].

**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. Due to the RCE structure of $Y$, this decomposition uses projections on orthogonal spaces generated by latent variables obtained through the Aldous-Hoover-Kallenberg (AHK) representation. In our asymptotic framework, the total number of rows and columns, denoted by $N$, tends to infinity. To treat this setting, we assume that the matrix dimensions $m_N$ and $n_N$ grow at the same rate, satisfying $N = m_N + n_N$, with $m_N/N \to \rho$ and $n_N/N \to 1 - \rho$, where $\rho \in (0, 1)$. The case in which the dimensions grow at different rates ($\rho = 0$ or $\rho = 1$) is not considered in this paper, although our results and proofs can be generalized to accommodate this scenario. For simplification, we denote $U_N^h(Y) := U_{m_N,n_N}^h(Y)$. First, we show that $\sqrt{N}(U_N^h(Y) - E[h])$ converges to a Gaussian distribution, except in degenerate cases. Then, for these U-statistics, we build a computationally efficient estimator for their variances. This variance estimator relies on the prior estimation of conditional expectations given the AHK latent variables, which makes it unique compared to traditional ways of estimating the variance of U-statistics.

**Outline** In Section 2, we recall the Aldous-Hoover-Kallenberg (AHK) representation for RCE matrices. In Section 3, we leverage the AHK representation to derive an orthogonal decomposition for U-statistics on RCE matrices. Section 4 exploits this new Hoeffding-type decomposition to establish the asymptotic normality of the U-statistics under consideration. In Section 5, we use this decomposition to construct an efficient estimator for the asymptotic variance. Section 6 demonstrates how our results can be applied to define a methodological framework for bipartite network analysis, providing examples of models for RCE matrices and network analysis questions that can be addressed. Sections 7 and 8 are dedicated to simulation studies and an illustration using a legislature dataset, showcasing this methodology.

## 2. Aldous-Hoover-Kallenberg representation

Traditionally, orthogonal decompositions for U-statistics on i.i.d. observations [28] or generalized U-statistics [31] are obtained by partitioning the probability space into orthogonal subspaces generated by subsets of observations. However, this approach generally fails for U-statistics on RCE matrices, as the dependencies between matrix elements can prevent orthogonality. To ensure orthogonality, the decomposition proposed in this paper uses the Aldous-Hoover-Kallenberg representation of RCE matrices, which expresses the matrix as a function of i.i.d. random variables. Although these variables are unobserved, they provide a natural framework for constructing suitable orthogonal subspaces. Before specifying these subspaces and deriving the decomposition, we first introduce this representation.

**Aldous-Hoover-Kallenberg representation** Corollary 7.23 of [32] states that for any dissociated RCE matrix $Y$, there exists $(\xi_i)_{i \ge 1}$, $(\eta_j)_{j \ge 1}$ and $(\zeta_{ij})_{i,j \ge 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 \le i, j < \infty$,

$$Y_{ij} \stackrel{\text{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 $i \in \mathcal{P}_p(\mathbb{N})$ and columns $j \in \mathcal{P}_q(\mathbb{N})$ can be written $h(Y_{i,j}) \stackrel{\text{a.s.}}{=} h_{\phi}((\xi_i)_{i \in i}; (\eta_j)_{j \in j}; (\zeta_{ij})_{i,j \in i})$, where $h_{\phi}$ is some function depending on $h$ and $\phi$.
**Projection sets of AHK variables** For $i' \in \mathcal{P}(\mathbb{N})$ and $j' \in \mathcal{P}(\mathbb{N})$, we define the sets of AHK variables defined as

$$A(i', j') := ((\xi_i)_{i \in i'}, (\eta_j)_{j \in j'}, (\zeta_{ij})_{\substack{i \in i' \\ j \in j'}}).$$

These sets play a crucial role, as they form the foundation of our U-statistic decomposition. More specifically, the decomposition will be defined through projections onto probability spaces generated by these sets of AHK variables. Consequently, we consider the $\sigma$-algebras given by

$$\mathcal{A}_{i',j'} := \sigma(A(i', j')) = \sigma((\xi_i)_{i \in i'}, (\eta_j)_{j \in j'}, (\zeta_{ij})_{\substack{i \in i' \\ j \in j'}}).$$

**Notations** 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 \ge 1}$, $(\eta_j)_{j \ge 1}$ and $(\zeta_{ij})_{\substack{i \ge 1 \\ j \ge 1}}$. 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_{i,j}) = h_{\phi}((\xi_i)_{i \in i}; (\eta_j)_{j \in j}; (\zeta_{ij})_{\substack{i \in i' \\ j \in j'}}). \quad (4)$$

Therefore, it follows from our notations that

$$\mathbb{E}[h(Y_{i,j}) | \mathcal{A}_{i',j'}] = \mathbb{E}[h_{\phi}((\xi_i)_{i \in i}; (\eta_j)_{j \in j}; (\zeta_{ij})_{\substack{i \in i' \\ j \in j'}}) | (\xi_i)_{i \in i'}; (\eta_j)_{j \in j'}; (\zeta_{ij})_{\substack{i \in i' \\ j \in j'}}].$$

For fixed sets $i'$ and $j'$, the quantity $\mathbb{E}[h(Y_{i,j}) | \mathcal{A}_{i',j'}]$ only depends on the elements shared by $i$ and $i'$ and the elements shared by $j$ and $j'$, and not on the other elements of $i$, $i'$, $j$ and $j'$. Suppose $r = \text{Card}(i \cap i')$ and $c = \text{Card}(j \cap j')$. Without loss of generality, we can assume that $i' \in \mathcal{P}_r(i)$ and $j' \in \mathcal{P}_c(j)$ so $\mathbb{E}[h(Y_{i,j}) | \mathcal{A}_{i',j'}]$ only depends on the $r$ elements of $i'$ and the $c$ elements of $j'$. Therefore, we can define the quantities $\psi^{r,c}h(Y_{i',j'})$ such that

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

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

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

despite $Y_{i',\emptyset}$ being undefined.

### 3. Decomposition for U-statistics on RCE matrices

In this section, we derive a decomposition for U-statistics on RCE matrices based on the AHK representation, using the subspaces generated by the specific projection sets that we have defined. This decomposition follows from orthogonal projections on these subspaces. In this section and the following sections, for elements of $\mathbb{N}^2$, $(x, y) \le (x', y')$ means that both $x \le x'$ and $y \le y'$; $(x, y) < (x', y')$ means that, in addition, $(x, y) \ne (x', y')$.
**Hoeffding-type decomposition of the kernel** 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_{i,j}) = \psi^{r,c}h(Y_{i,j}) - \sum_{(0,0) \le (r',c') < (r,c)} \sum_{\substack{i' \in \mathcal{P}_{r'}(\mathbf{i}) \\ j' \in \mathcal{P}_{c'}(\mathbf{j})}} p^{r',c'}h(Y_{i',j'}). \quad (5)$$

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

$$h(Y_{i,j}) = \sum_{(0,0) \le (r,c) \le (p,q)} \sum_{\substack{i' \in \mathcal{P}_r(\mathbf{i}) \\ j' \in \mathcal{P}_c(\mathbf{j})}} p^{r,c}h(Y_{i',j}). \quad (6)$$

*Remark 1.* From this formula, we see that $h(Y_{i,j})$ is a linear combination of the projections $p^{r',c'}h(Y_{i',j'})$, for $0 \le r' \le p, 0 \le c' \le q, i' \in \mathcal{P}_{r'}(\mathbf{i}), j' \in \mathcal{P}_{c'}(\mathbf{j})$. Therefore, this is a linear combination of the respective conditional expectations $\psi^{r',c'}h(Y_{i',j'})$.

*Remark 2.* In practice, the AHK variables $(\xi_i)_{i \ge 1}, (\eta_j)_{j \ge 1}$ and $(\xi_{ij})_{i \ge 1, j \ge 1}$ are usually unobserved. However, the existence of such variables is sufficient to define the conditional expectations of the form $\psi^{r,c}h(Y_{i,j})$.

Now, we show that $p^{r,c}h(Y_{i,j})$ is the projection of $Y_{i,j}$ on the probability space generated by the projection set of AHK variables $A(\mathbf{i}, \mathbf{j})$, orthogonally to all the spaces generated by the sets $A(\mathbf{i}', \mathbf{j}')$, for $\mathbf{i}' \subset \mathbf{i}, \mathbf{j}' \subset \mathbf{j}$. This system of projection is analogous to the Hoeffding decomposition for the kernel functions of usual U-statistics on i.i.d. data. The next proposition, proven in Appendix B, states that the following properties hold, ensuring the orthogonality of the projection spaces generated by $A(\mathbf{i}_1, \mathbf{j}_1)$ and $A(\mathbf{i}_2, \mathbf{j}_2)$, if $(\mathbf{i}_1, \mathbf{j}_1) \ne (\mathbf{i}_2, \mathbf{j}_2)$.

**Proposition 3.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. Let $(0,0) \le (r_1, c_1) \le (p_1, q_1)$ and $(0,0) \le (r_2, c_2) \le (p_2, q_2)$ such that $(r_1, c_1) \ne (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_{i_1, j_1}), p^{r_2, c_2}h_2(Y_{i_2, j_2})) = 0. $$

2. Let $(r, c)$ such that $(0,0) \le (r, c) \le (p_1, q_1)$ and $(0,0) \le (r, c) \le (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) \ne (\mathbf{i}_2, \mathbf{j}_2)$, then
$$ \mathrm{Cov}(p^{r, c}h_1(Y_{i_1, j_1}), p^{r, c}h_2(Y_{i_2, j_2})) = 0. $$

**Decomposition of U-statistics** Using the Hoeffding-type decomposition of kernel functions (6) the U-statistic (2) can be reformulated as:

$$ \begin{aligned} U_{m,n}^h(Y) &= \left(\frac{m}{p}\right)^{-1} \left(\frac{n}{q}\right)^{-1} \sum_{\substack{1 \le i_1 < \dots < i_p \le m \\ 1 \le j_1 < \dots < j_q \le n}} \sum_{\substack{i \in \mathcal{P}_r(\{i_1, \dots, i_p\}) \\ j \in \mathcal{P}_c(\{j_1, \dots, j_q\})}} \sum_{\mathbf{p} \in \mathcal{P}_r(\llbracket m \rrbracket)} p^{r,c}h(Y_{i,j}) \\ &= \left(\frac{m}{p}\right)^{-1} \left(\frac{n}{q}\right)^{-1} \sum_{(0,0) \le (r,c) \le (p,q)} \binom{m-r}{p-r} \binom{n-c}{q-c} \sum_{\substack{i \in \mathcal{P}_r(\llbracket m \rrbracket) \\ j \in \mathcal{P}_c(\llbracket n \rrbracket)}} p^{r,c}h(Y_{i,j}) \\ &= \sum_{(0,0) \le (r,c) \le (p,q)} \binom{p}{r} \binom{q}{c} P_{m,n}^{r,c}h(Y), \end{aligned} \quad (7) $$
where for all $0 \le r \le p$ and $0 \le c \le q$,

$$P_{m,n}^{r,c}h(Y) = \binom{m}{r}^{-1} \binom{n}{c}^{-1} \sum_{\substack{i \in \mathcal{P}_r(\llbracket m \rrbracket) \\ j \in \mathcal{P}_c(\llbracket n \rrbracket)}} p^{r,c}h(Y_{i,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 by Corollary B.2 provided in Appendix B. The orthogonality between the $P_{m,n}^{r_1, c_1}h_1(Y)$ and $P_{m,n}^{r_2, c_2}h_2(Y)$ allows to decompose the covariance of two U-statistics into a few covariance terms, as specified by Corollary 3.2.

**Corollary 3.2.**

$$\operatorname{Cov}(U_{m,n}^{h_1}(Y), U_{m,n}^{h_2}(Y)) = \sum_{(0,0)<(r,c) \le (p,q)} \binom{p}{r}^2 \binom{q}{c}^2 \binom{m}{r}^{-1} \binom{n}{c}^{-1} \operatorname{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)).$$

This variance decomposition is particularly noteworthy and useful for analyzing the asymptotic properties of U-statistics. In contrast, the classic decomposition of (unipartite) graph U-statistics is typically more complex [31]. Importantly, with Corollary 3.2, each term in the decomposition (7) has a distinct asymptotic order, providing a general understanding of variance contributions. In this regard, our decomposition more closely resembles the traditional decomposition of U-statistics on i.i.d. observations, where each covariance term is associated with a distinct binomial coefficient that determines its order of magnitude, rather than the decomposition in [31]. In addition, the bipartite setting uniquely allows for separate contributions from the rows and the columns of observations. Later, the simplicity of this decomposition is exploited to introduce a novel variance estimator for these U-statistics.

**4. Asymptotic normality**

This section establishes the following Central Limit Theorem for $U_N^h$. In this section and the following sections, we will use simplified notations, summarizing the couple $(m_N, n_N)$ into $N = m_N + n_N$. We recall that $U_N^h(Y) = U_{m_N, n_N}^h(Y)$. We also denote $P_N^{r,c}h(Y) := P_{m_N, n_N}^{r,c}h(Y)$. When this is unambiguous, we will omit to mention $Y$, so we will simply write $U_N^h$, $P_N^{r,c}h$, $p_{(i,j)}^{r,c}h$ and $\psi_{(i,j)}^{r,c}h$ instead of $U_N^h(Y)$, $P_N^{r,c}h(Y)$, $p^{r,c}h(Y_{i,j})$ and $\psi^{r,c}h(Y_{i,j})$. Denote $v_h^{r,c} := \mathbb{V}[\psi_{(\llbracket r \rrbracket, \llbracket c \rrbracket)}^{r,c}]^h$.

**Theorem 4.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 \ge 1}$ be a sequence of dimensions for the U-statistics, such that $\frac{m_N}{N} \xrightarrow{N \to \infty} \rho$ and $\frac{n_N}{N} \xrightarrow{N \to \infty} 1 - \rho$, where $\rho \in (0, 1)$. Let $(U_N^h)_{N \ge 1}$ be the sequence of U-statistics associated with $h$ defined by $U_N^h := U_{m_N, n_N}^h$. Set $U_\infty^h = \mathbb{E}[h(Y_{(1,...,p;1,...,q)})]$ and

$$V^h = \frac{p^2}{\rho} v_h^{1,0} + \frac{q^2}{1-\rho} v_h^{0,1}.$$

If $V^h > 0$, then

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

**Lemma 4.2.** If $v_h^{1,0} > 0$, then we have

$$ \frac{1}{\sqrt{m}} \sum_{i=1}^{m} p_{(\{i\}, \emptyset)}^{1,0} h \xrightarrow{\mathcal{D}} \mathcal{N}(0, v_h^{1,0}), $$

and if $v_h^{0,1} > 0$, then we have

$$ \frac{1}{\sqrt{n}} \sum_{j=1}^{n} p_{(\emptyset, \{j\})}^{0,1} h \xrightarrow{\mathcal{D}} \mathcal{N}(0, v_h^{0,1}). $$

**Lemma 4.3.** Let $A_N := \sqrt{N} \sum_{(0,0)<(r,c) \le (p,q)} (p_r)(q_c) P_N^{r,c} h$. Then $A_N \xrightarrow{\mathbb{P}} 0$.

*Proof of Theorem 4.1.* We have

$$ \begin{aligned} U_N^h &= \sum_{(0,0) \le (r,c) \le (p,q)} \binom{p}{r} \binom{q}{c} P_N^{r,c} h \\ &= P_N^{0,0} h + p P_N^{1,0} h + q P_N^{0,1} h + \sum_{\substack{(0,0)<(r,c) \le (p,q) \\ (r,c) \ne (1,0) \ne (0,1)}} \binom{p}{r} \binom{q}{c} P_N^{r,c} h. \end{aligned} $$

First, we see that $P_N^{0,0}h = U_\infty^h$. Next, $A_N$ being defined in Lemma 4.3, we have

$$ \sqrt{N}(U_N^h - U_\infty^h) = \frac{\sqrt{N}p}{m_N} \sum_{i=1}^{m_N} p_{(\{i\}, \emptyset)}^{1,0} h + \frac{\sqrt{N}q}{n_N} \sum_{j=1}^{n_N} p_{(\emptyset, \{j\})}^{0,1} h + A_N. $$

From Lemma 4.3, we have $A_N \xrightarrow{\mathbb{P}} 0$. So by Slutsky's theorem, $\sqrt{N}(U_N^h - U_\infty^h)$ has the same limiting distribution as the two main terms of this decomposition. From Lemma 4.2, this is the sum of two centered Gaussians of respective variance $\frac{p^2}{\rho}v_h^{1,0}$ and $\frac{q^2}{1-\rho}v_h^{0,1}$. Furthermore, $\sum_{i=1}^m p_{(\{i\}, \emptyset)}^{1,0}h$ and $\sum_{j=1}^n p_{(\{j\}, \emptyset)}^{0,1}h$ are independent, so the two Gaussians are independent, which concludes the proof. □

**Remark 3.** The expression of $V^h$ could have been predicted with Corollary 3.2. Indeed, this corollary implies that

$$ \mathbb{V}[U_N^h] = \sum_{(0,0)<(r,c) \le (p,q)} \binom{p}{r}^2 \binom{q}{c}^2 \binom{m_N}{r}^{-1} \binom{n_N}{c}^{-1} \mathbb{V}[p_{(\llbracket r \rrbracket, \llbracket c \rrbracket)}^{r,c}] h, $$

so

$$ \begin{aligned} \lim_{N \to \infty} N \mathbb{V}[U_N^h] &= \lim_{N \to \infty} \left( \frac{p^2 N}{m_N} \mathbb{V}[p_{(\{1\}, \emptyset)}^{1,0} h] + \frac{q^2 N}{n_N} \mathbb{V}[p_{(\emptyset, \{1\})}^{0,1} h] \right) \\ &= \frac{p^2}{c} \mathbb{V}[\psi_{(\{1\}, \emptyset)}^{1,0} h] + \frac{q^2}{1-c} \mathbb{V}[\psi_{(\emptyset, \{1\})}^{0,1} h]. \end{aligned} $$
**Remark 4.** It is worth noting that, if $V^h = 0$, this convergence result remains valid, despite the theorem excluding this case. In such a scenario, the limit would be trivial:

$$ \sqrt{N}(U_N^h - U_\infty^h) \xrightarrow[N \to \infty]{\mathbb{P}} 0. $$

This case is referred to as degeneracy, and we say that $U_N^h$ is a degenerate U-statistic. From the orthogonal decomposition of the U-statistic (7), each component $P_N^{r,c}$ contributes to the limit distribution with a magnitude determined by $(r, c)$. Lemma 4.2 establishes that the contributions of the two dominant terms are Gaussian. Degeneracy occurs when these components vanish, implying that the shape of the limit distribution is determined by higher-order terms, which requires a larger scaling than $\sqrt{N}$ to observe. In this paper, we only consider non-degenerate U-statistics.

**Remark 5.** The asymptotic distribution of the U-statistic is determined by the first-order projections of the form $p_{\{i\},\emptyset}^{1,0}h$ and $p_{\emptyset,\{j\}}^{1,0}h$, through their separate contributions given by Lemma 4.2. These projections depend on the AHK latent variables $(\xi_i)$ and $(\eta_j)$, which correspond to the row and column effects in the RCE matrix. Interestingly, the entry-specific AHK variables $(\zeta_{ij})$ do not contribute to the asymptotic distribution. To build intuition for this observation, consider the representation of the U-statistic in terms of AHK variables, derived from (4):

$$ U_N^h = \binom{m_N}{r}^{-1} \binom{n_N}{c}^{-1} \sum_{\substack{i \in \mathcal{P}_p(\lfloor m_N \rfloor) \\ j \in \mathcal{P}_q(\lfloor n_N \rfloor)}} h_{\phi((\xi_i)_{i \in i}; (\eta_j)_{j \in j}; (\zeta_{ij})_{i \in i})}. $$

In this summation over **i** and **j**, the entry-specific variable $\zeta_{ij}$ appears only in terms where both $i \in i$ and $j \in j$. In contrast, the row-specific variable $(\xi_i)$ contributes to all terms where $i \in i$, regardless of the choice of $j \in j$. Since the number of such terms is $O(N)$ larger, the influence of $\xi_i$ on $U_N^h$ dominates that of $\zeta_{ij}$. An analogous reasoning applies to column-specific variables $(\eta_j)$. Thus, in the asymptotic regime, the contribution of entry-specific fluctuations captured by $(\zeta_{ij})$ is of higher order and is negligible compared to the first-order row-column variability captured by $(\xi_i)$ and $(\eta_j)$.

**Remark 6.** Theorem 4.1 can be extended to cases where the numbers of rows and columns grow at different rates, meaning that $\rho = 0$ or $\rho = 1$. For instance, consider the case $\rho = 0$, which implies $m_N = o(n_N)$. In this scenario, we can rescale $U_N^h$ by $\sqrt{m_N}$ instead of $\sqrt{N}$ in the proof of Theorem 4.1, leading to

$$ \sqrt{m_N}(U_N^h - U_\infty^h) = \frac{p}{\sqrt{m_N}} \sum_{i=1}^{m_N} p_{\{i\},\emptyset}^{1,0} h + B_N, $$

where Lemma 4.2 remains applicable and, analogous to Lemma 4.3, it can be shown that the sequence $B_N \xrightarrow[N \to \infty]{\mathbb{P}} 0$. Therefore, we obtain

$$ \sqrt{m_N}(U_N^h - U_\infty^h) \xrightarrow[N \to \infty]{\mathcal{D}} \mathcal{N}(0, V'), $$

where $V' = p^2 v_h^{1,0}$.
In Appendix C, we also demonstrate the asymptotic normality of vectors and functions of U-statistics on RCE matrices, as stated by Corollaries C.3 and C.4. Specifically, the asymptotic normality of functions of U-statistics is very useful to perform statistical inference on bipartite network data. We give more details in Section 6.

## 5. Estimation of the asymptotic variance

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

$$ \mathbb{V}[U_N^h] = \left[ \begin{pmatrix} m \\ p \end{pmatrix} \begin{pmatrix} n \\ q \end{pmatrix} \right]^{-2} \sum_{\substack{i \in \mathcal{P}_p(\llbracket m \rrbracket) \\ j \in \mathcal{P}_q(\llbracket n \rrbracket)}} \sum_{\substack{i' \in \mathcal{P}_p(\llbracket m \rrbracket) \\ j' \in \mathcal{P}_q(\llbracket n \rrbracket)}} \text{Cov}(h(Y_{i,j}), h(Y_{i',j'})). $$

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

$$ \mathbb{V}[U_N^h] = \left[ \begin{pmatrix} m \\ p \end{pmatrix} \begin{pmatrix} n \\ q \end{pmatrix} \right]^{-1} \sum_{(0,0) \le (r,c) \le (p,q)} \begin{pmatrix} p \\ r \end{pmatrix} \begin{pmatrix} q \\ c \end{pmatrix} \begin{pmatrix} m-p \\ p-r \end{pmatrix} \begin{pmatrix} n-q \\ q-c \end{pmatrix} \gamma_h^{r,c}. $$

Each $\gamma_h^{r,c}$, $1 \le r \le p$, $1 \le c \le 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_N^h]$ similar to the one discussed by [51] 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 4.1 is $V^h = \frac{p^2}{\rho} v_h^{1,0} + \frac{q^2}{1-\rho} v_h^{0,1}$. Remark that $v_h^{1,0} = \gamma_h^{1,0}$ and $v_h^{0,1} = \gamma_h^{0,1}$. 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 [37].

We present here a kernel-free and model-free estimator of $V^h$, taking advantage of the Hoffding decomposition. Indeed, we will first define estimators for the conditional expectations $(\psi_{\{\{i\},\emptyset\}}^{1,0} h)_{1 \le i \le m}$ and $(\psi_{(\emptyset,\{j\})}^{0,1} h)_{1 \le j \le n}$. Then, using the fact that $v_h^{1,0} = \mathbb{V}[\psi_{\{\{1\},\emptyset\}}^{1,0} h]$ and $v_h^{0,1} = \mathbb{V}[\psi_{(\emptyset,\{1\})}^{0,1} h]$, we can derive a positive estimator for $V^h$. Afterward, we explain how to generalize it to estimate the variance of functions of U-statistics.

### 5.1. Notations

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 \le i \le m_N$ and $j$ being a column index means that $1 \le j \le n_N$. We denote

$$ X_{i,j} := h(Y_{i,j}). $$
For $N$ such that $m_N \ge p$ and $n_N \ge q$, we further denote

$$S_N^{p,q} := \{ (\mathbf{i}, \mathbf{j}) : \mathbf{i} \in \mathcal{P}_p(\llbracket m_N \rrbracket), \mathbf{j} \in \mathcal{P}_q(\llbracket n_N \rrbracket) \},$$

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 S_N^{p,q}}$.

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

$$S_{N, (\underline{\mathbf{i}}, \underline{\mathbf{j}})}^{p,q} := \{ (\mathbf{i}, \mathbf{j}) \in S_N^{p,q} : \underline{\mathbf{i}} \subset \mathbf{i}, \underline{\mathbf{j}} \subset \mathbf{j} \}.$$

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 S_{N, (\emptyset, \{1,2\})}^{p,q}}$.

## 5.2. Estimation of the conditional expectations

In this section, for all $i \in \llbracket m_N \rrbracket$ and $j \in \llbracket n_N \rrbracket$, we define estimators for $\psi_{(\{i\},\emptyset)}^{1,0} h = \mathbb{E}[X_{\llbracket p \rrbracket, \llbracket q \rrbracket} | \xi_i]$ and $\psi_{(\emptyset,\{j\})}^{0,1} h = \mathbb{E}[X_{\llbracket p \rrbracket, \llbracket q \rrbracket} | \eta_j]$, where $\xi_i$ and $\eta_j$ have been defined in Section 2.

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

$$\widehat{\mu}_N^{h,(i)} := \binom{m_N - 1}{p - 1}^{-1} \binom{n_N}{q}^{-1} \sum_{(\mathbf{i}, \mathbf{j}) \in S_{N, (\{i\}, \emptyset)}^{p,q}} X_{\mathbf{i}, \mathbf{j}}, \quad (8)$$

and

$$\widehat{\nu}_N^{h,(j)} := \binom{m_N}{p}^{-1} \binom{n_N - 1}{q - 1}^{-1} \sum_{(\mathbf{i}, \mathbf{j}) \in S_{N, (\emptyset, \{j\})}^{p,q}} X_{\mathbf{i}, \mathbf{j}}. \quad (9)$$

The following propositions, proven in Appendix D, establish guarantees for these estimators.

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

* $\mathbb{E}[\widehat{\mu}_N^{h,(i)} | \xi_i] = \psi_{(\{i\},\emptyset)}^{1,0} h$,

* $\mathbb{E}[\widehat{\nu}_N^{h,(j)} | \eta_j] = \psi_{(\emptyset,\{j\})}^{0,1} h$.

**Proposition 5.2.** If $Y$ is an RCE matrix, then

* $\widehat{\mu}_N^{h,(i)} \xrightarrow[N \to \infty]{a.s., L_1} \psi_{(\{i\},\emptyset)}^{1,0} h$,

* $\widehat{\nu}_N^{h,(j)} \xrightarrow[N \to \infty]{a.s., L_1} \psi_{(\emptyset,\{j\})}^{0,1} h$.

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

**Remark 7.** Despite the AHK variables $(\xi_i)_{i \ge 1}$ and $(\eta_j)_{j \ge 1}$ being unobserved, Propositions 5.1 and 5.2 ensure that the conditional expectations $\psi_{(\{i\},\emptyset)}^{1,0} h$ and $\psi_{(\emptyset,\{j\})}^{0,1} h$ can be estimated using the estimators $\widehat{\mu}_N^{h,(i)}$ and $\widehat{\nu}_N^{h,(j)}$.
### 5.3. Estimation of the asymptotic variance

Using the estimators for the conditional expectations $\psi_{\{(i),\emptyset\}}^{1,0}h$ and $\psi_{\{(\emptyset,\{j\})}^{0,1}h$ that we have defined in the previous section, we suggest the following estimators for $v_h^{1,0} = \mathbb{V}[\psi_{\{(1),\emptyset\}}^{1,0}h]$ and $v_h^{0,1} = \mathbb{V}[\psi_{\{(\emptyset,\{j\})}^{0,1}h]$:

$$ \hat{v}_N^{h;1,0} = \left(\binom{m_N}{2}\right)^{-1} \sum_{1 \le i_1 < i_2 \le m_N} \frac{(\hat{\mu}_N^{h,(i_1)} - \hat{\mu}_N^{h,(i_2)})^2}{2}, \quad (10) $$

$$ \hat{v}_N^{h;0,1} = \left(\binom{n_N}{2}\right)^{-1} \sum_{1 \le j_1 < j_2 \le n_N} \frac{(\hat{v}_N^{h,(j_1)} - \hat{v}_N^{h,(j_2)})^2}{2}. \quad (11) $$

Then, an estimator for $V^h$ is

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

The following theorem shows that $\hat{V}_N^h$ is a consistent estimator for $V^h$. The proof relies on the fact that $\hat{v}_N^{h;1,0}$ and $\hat{v}_N^{h;0,1}$ are both asymptotically unbiased and have vanishing variances. These properties are established by Propositions D.2 and D.3 in Appendix D.

**Theorem 5.3.** We have $\hat{v}_N^{h;1,0} \xrightarrow[N \to \infty]{} v_h^{1,0}$ and $\hat{v}_N^{h;0,1} \xrightarrow[N \to \infty]{} v_h^{0,1}$. As a consequence, $\hat{v}_N^h \xrightarrow[N \to \infty]{} V^h$.

*Proof of Theorem 5.3.* For some $\epsilon > 0$, it follows from Proposition D.2 that for large enough values of $N$, $|\mathbb{E}[\hat{v}_N^{h;1,0}] - v_h^{1,0}| < \epsilon$. The triangular inequality and Chebyshev's inequality state that

$$ \begin{aligned} \mathbb{P}(|\hat{v}_N^{h;1,0} - v_h^{1,0}| > \epsilon) &\le \mathbb{P}(|\hat{v}_N^{h;1,0} - \mathbb{E}[\hat{v}_N^{h;1,0}]| \ge \epsilon - |\mathbb{E}[\hat{v}_N^{h;1,0}] - v_h^{1,0}|) \\ &\le \frac{\mathbb{V}[\hat{v}_N^{h;1,0}]}{(\epsilon - |\mathbb{E}[\hat{v}_N^{h;1,0}] - v_h^{1,0}|)^2}. \end{aligned} $$

Applying Propositions D.2 and D.3 to the right-hand side of the inequality ensures that $\mathbb{P}(|\hat{v}_N^{h;1,0} - v_h^{1,0}| > \epsilon) \xrightarrow[N \to \infty]{} 0$ which concludes the proof. $\square$

With Theorem 5.3, it is possible to use $\hat{V}_N^h$ for statistical inference tasks when plugged-in in place of $V^h$, an asymptotic normality result similar to Theorem 4.1 holds, by Slutsky's theorem.

**Proposition 5.4.** If $V^h > 0$, then

$$ \sqrt{\frac{N}{\hat{V}_N^h}} (U_N^h - U_\infty^h) \xrightarrow{\mathcal{D}} \mathcal{N}(0, 1). $$
## 5.4. Computation of the estimator

In practice, the computation of the estimators $(\hat{\mu}_N^{h,(i)})_{1\le i\le m_N}$ and $(\hat{\nu}_N^{h,(j)})_{1\le j\le n_N}$ using their definition is pretty tedious, as they are sums of $O(N^{p+q-1})$ terms. The computation cost of $\hat{V}_N^h$ is then $O(N^{p+q})$ which is of the same order as the computation of $U_N^h$ when naively applying the kernel function to all the $p \times q$ submatrices and averaging. However, computing $U_N^h$ in this manner is generally undesirable. In many problem-specific settings, more efficient methods exist for computing U-statistics. Typically, for simple kernels, $U_N^h$ can be expressed in terms of matrix operations, which are optimized in most computing software (see Section 7.4 for an example). In such cases, it is beneficial to rewrite the proposed variance estimators in the form of those U-statistics, to leverage these advantages.

**Proposition 5.5.** An alternative form for $\hat{v}_N^{h;1,0}$ and $\hat{v}_N^{h;0,1}$ is given by

$$ \hat{v}_N^{h;1,0} = \frac{(m_N - p)^2}{p^2(m_N - 1)} \sum_{i=1}^{m_N} (U_N^h - U_N^{h,(-i,0)})^2 \quad (12) $$

and

$$ \hat{v}_N^{h;0,1} = \frac{(n_N - q)^2}{q^2(n_N - 1)} \sum_{j=1}^{n_N} (U_N^h - U_N^{h,(0,-j)})^2, \quad (13) $$

where

$$ U_N^{h,(-i,0)} := \left[ \begin{pmatrix} m_N - 1 \\ p \end{pmatrix} \begin{pmatrix} n_N \\ q \end{pmatrix} \right]^{-1} \sum_{\substack{i \in \mathcal{P}_p(\llbracket m_N \rrbracket \setminus \{i\}) \\ j \in \mathcal{P}_q(\llbracket n_N \rrbracket)}} h(Y_{i,j}) $$

and

$$ U_N^{h,(0,-j)} := \left[ \begin{pmatrix} m_N \\ p \end{pmatrix} \begin{pmatrix} n_N - 1 \\ q \end{pmatrix} \right]^{-1} \sum_{\substack{i \in \mathcal{P}_p(\llbracket m_N \rrbracket) \\ j \in \mathcal{P}_q(\llbracket n_N \rrbracket \setminus \{j\})}} h(Y_{i,j}). $$

*Proof.* The relations are obtained from the definition of $\hat{v}_N^{h;1,0}$ and $\hat{v}_N^{h;0,1}$ and noticing that

$$ \begin{pmatrix} m_N - 1 \\ p - 1 \end{pmatrix} \begin{pmatrix} n_N \\ q \end{pmatrix} \hat{\mu}_N^{h,(i)} = \begin{pmatrix} m_N \\ p \end{pmatrix} \begin{pmatrix} n_N \\ q \end{pmatrix} U_N^h - \begin{pmatrix} m_N - 1 \\ p \end{pmatrix} \begin{pmatrix} n_N \\ q \end{pmatrix} U_N^{h,(-i,0)} $$

and

$$ \begin{pmatrix} m_N \\ p \end{pmatrix} \begin{pmatrix} n_N - 1 \\ q - 1 \end{pmatrix} \hat{v}_N^{h,(j)} = \begin{pmatrix} m_N \\ p \end{pmatrix} \begin{pmatrix} n_N \\ q \end{pmatrix} U_N^h - \begin{pmatrix} m_N \\ p \end{pmatrix} \begin{pmatrix} n_N - 1 \\ q \end{pmatrix} U_N^{h,(0,-j)}. \quad \square $$

*Remark 8.* This is one alternative method to compute $\hat{V}_N^h$, but not necessarily the optimal one. With this form, the computational cost of $\hat{V}_N^h$ is $O(N\Gamma(N))$, assuming the cost of computing a U-statistic is $O(\Gamma(N))$, where $\Gamma$ is an increasing function. This can be outperformed by the naive method in some situations, for example when $N^{p+q} = O(\Gamma(N))$.
Remark 9. The alternative form of $\hat{v}_N^{h;1,0}$ and $\hat{v}_N^{h;0,1}$ 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 \to \mathbb{R}$ is

$$U_N^h = \left(\frac{N}{p}\right)^{-1} \sum_{i \in \mathcal{P}_p(\llbracket N \rrbracket)} h(X_i).$$

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

$$\hat{V}_N^{h,J} = (N-1) \sum_{i=1}^{N} (U_N^h - U_N^{h,(-i)})^2,$$

where

$$U_N^{h,(-i)} = \left(\frac{N}{p}\right)^{-1} \sum_{i \in \mathcal{P}_p(\llbracket N \rrbracket \setminus \{i\})} h(X_i).$$

In fact, our estimator is closer to Sen's estimator of the asymptotic variance of [52, 53] which depends on the kernel size and is related to the jackknife estimator $\hat{V}_N^{h,S} = (N-k)^2/(N-1)^2\hat{V}_N^{h,J}$. 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)}$.

### 5.5. Extension to functions of U-statistics

In the case of a function of U-statistics $g(U_N^{h_1}, \dots, U_N^{h_D})$, the asymptotic normality result of Corollary C.4 applies and the asymptotic variance to be estimated is

$$V^{\delta} = \nabla g(\theta)^T \Sigma^{h_1, \dots, h_D} \nabla g(\theta) \neq 0,$$

where $\theta = (U_{\infty}^{h_1}, \dots, U_{\infty}^{h_D})$. Similar to the estimator for the asymptotic variance of $V^h$, we suggest an estimator for the covariance matrix $\Sigma^{h_1, \dots, h_D} = (C^{h_k, h_l})_{1 \le k, l \le D}$.

For each kernel $h_k$, $1 \le k \le D$, let $\hat{\mu}_N^{h_k,(i)}$ and $\hat{v}_N^{h_k,(j)}$ be the respective estimators of the conditional expectations $\psi_{(\{i\},\emptyset)}^{1,0} h_k$ and $\psi_{(\emptyset,\{j\})}^{0,1} h_k$, as defined in equations (8) and (9). Define

$$\hat{c}_N^{h_k, h_l; 1,0} := \left(\frac{m_N}{2}\right)^{-1} \sum_{1 \le i_1 < i_2 \le m_N} \frac{(\hat{\mu}_N^{h_k, (i_1)} - \hat{\mu}_N^{h_k, (i_2)}) (\hat{\mu}_N^{h_l, (i_1)} - \hat{\mu}_N^{h_l, (i_2)})}{2}$$

and

$$\hat{c}_N^{h_k, h_l; 0,1} := \left(\frac{n_N}{2}\right)^{-1} \sum_{1 \le j_1 < j_2 \le n_N} \frac{(\hat{v}_N^{h_k, (j_1)} - \hat{v}_N^{h_k, (j_2)}) (\hat{v}_N^{h_l, (j_1)} - \hat{v}_N^{h_l, (j_2)})}{2}.$$

Then, for two kernels $h_k$ and $h_l$,

$$\hat{C}_N^{h_k, h_l} := \frac{p^2}{\rho} \hat{c}_N^{h_k, h_l; 1,0} + \frac{q^2}{1-\rho} \hat{c}_N^{h_k, h_l; 0,1},$$

is an estimator of the asymptotic covariance term $C^{h_k, h_l}$. With a similar proof as for Theorem 5.3, the following theorem ensures the consistency of this estimator.
**Theorem 5.6.** For two kernel functions $h_k$ and $h_\ell$, we have $\tilde{c}_N^{h_k,h_\ell;1,0} \xrightarrow[N \to \infty]{P} c_{h_k,h_\ell}^{1,0}$ and $\tilde{c}_N^{h_k,h_\ell;0,1} \xrightarrow[N \to \infty]{P} c_{h_k,h_\ell}^{0,1}$. As a consequence, $\tilde{C}_N^{h_k,h_\ell} \xrightarrow[N \to \infty]{P} C^{h_k,h_\ell}$.

Thus, for linearly independent kernel functions $(h_1, ..., h_D)$, the entries of the matrix $\tilde{\Sigma}_N^{h_1, ..., h_D} := (\tilde{C}_N^{h_1, h_\ell})_{1 \le k, \ell \le D}$ converge to the entries of $\Sigma^{h_1, ..., h_D}$. Set

$$ \tilde{V}_N^\delta := \nabla g(U_N^{h_1}, ..., U_N^{h_D})^T \tilde{\Sigma}_N^{h_1, ..., h_D} \nabla g(U_N^{h_1}, ..., U_N^{h_D}), $$

then we have $\tilde{V}_N^\delta \xrightarrow[N \to \infty]{P} V^\delta$. Finally, the “studentized” statistic obtained by plug-in also stands.

**Proposition 5.7.** If $V^\delta > 0$, then

$$ \sqrt{\frac{N}{\tilde{V}_N^\delta}} (g(U_N^{h_1}, ..., U_N^{h_D}) - g(\theta)) \xrightarrow[N \to \infty]{D} \mathcal{N}(0, 1). $$

# 6. Methodology for bipartite network analysis

The previous section introduced a generic and consistent estimator for the asymptotic variance of U-statistics on RCE matrices. This estimator can be used to develop a methodology for statistical inference on bipartite network data. In particular, the RCE framework, combined with the versatility provided by Proposition 5.7 in selecting network statistics, supports a broad range of network analysis tasks through model parameter estimation and hypothesis testing. In this section, we first examine the network comparison problem as a statistical inference task. We then illustrate the flexibility of this framework by presenting examples of RCE models and U-statistic-based network statistics applicable to various network analysis problems.

## 6.1. Example of inference task: Network comparison

Network comparison has a long history in network analysis literature [23, 62]. Few comparison methods use random network models [3, 39, 37]. However, there are several advantages to using model-based approaches. Indeed, random models define a probability distribution of networks that 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_N^A \times n_N^A$ and $m_N^B \times n_N^B$. 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 $\hat{\theta}_N$ be an estimator for this quantity of interest. If $\hat{\theta}_N$ is a U-statistic or a function of U-statistics, then as previously, we have both

$$ \sqrt{N}(\hat{\theta}_N(Y^A) - \theta^A) \xrightarrow[N \to \infty]{D} \mathcal{N}(0, V^A) $$
and

$$ \sqrt{N}(\hat{\theta}_N(Y^B) - \theta^B) \xrightarrow[N \to \infty]{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 $H_0: \theta^A = \theta^B$ and $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}} \left( \hat{\delta}_N(Y^A, Y^B) - (\theta^A - \theta^B) \right) \xrightarrow[N \to \infty]{D} \mathcal{N}(0, 1), $$

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

$$ Z_N(Y^A, Y^B) := \sqrt{\frac{N}{\hat{V}_N^A + \hat{V}_N^B}} \hat{\delta}_N(Y^A, Y^B) $$

admits the convergence result

$$ Z_N(Y^A, Y^B) - \sqrt{\frac{N}{\hat{V}_N^A + \hat{V}_N^B}} (\theta^A - \theta^B) \xrightarrow[N \to \infty]{D} \mathcal{N}(0, 1). $$

Under $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.2. Examples of RCE models for networks

**Bipartite Expected Degree Distribution model** The Bipartite Expected Degree Distribution (BEDD) model, suggested by [47], is a binary graph model characterized 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 \ge 1}$ and $(\eta_j)_{j \ge 1}$ corresponding to the row and column-nodes of the graph:

$$ \begin{gathered} \xi_i, \eta_j \stackrel{\text{i.i.d.}}{\sim} \mathcal{U}[0, 1] \tag{14} \\ Y_{ij} \mid \xi_i, \eta_j \sim \mathcal{B}(\lambda f(\xi_i) g(\eta_j)). \end{gathered} $$

where

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

* $\lambda$ is a positive real number such that $\lambda \le \|f\|_{\infty}^{-1} \|g\|_{\infty}^{-1}$,

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

The BEDD model is an RCE version of the Expected Degree Sequence model of [16] 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) [26] is a binary graph model characterized by a partition of row and column-nodes in several groups. It can be considered as a bipartite extension of the Stochastic Block model [44]. The probability of interaction between two nodes is fully determined by the groups to which they belong. All the nodes have the same probability of belonging to each group. The distribution of an LBM is most commonly written using independent latent variables for the node attribution in a group ($Z_i$)$_{i \ge 1}$ and ($W_j$)$_{j \ge 1}$ corresponding to the row and column-nodes of the graph:

$$ 
\begin{align} 
Z_i &\stackrel{\text{i.i.d.}}{\sim} \mathcal{M}(1; \alpha) \nonumber \\ 
W_j &\stackrel{\text{i.i.d.}}{\sim} \mathcal{M}(1; \beta) \tag{15} \\ 
Y_{ij} | Z_i = k, W_j = \ell &\sim \mathcal{B}(\pi_{k\ell}), \nonumber 
\end{align} 
$$

where $\alpha = (\alpha_1, ..., \alpha_K)$ and $\beta = (\beta_1, ..., \beta_L)$ are the probability vectors of the rows and the columns and $\pi = (\pi_{k\ell})_{1 \le k \le K, 1 \le \ell \le L} \in [0, 1]^{KL}$ is a matrix of probabilities. The LBM is an RCE model since the group attribution variables of the nodes are exchangeable.

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

$$ 
\begin{align} 
\xi_i, \eta_j &\stackrel{\text{i.i.d.}}{\sim} \mathcal{U}[0, 1] \nonumber \\ 
Y_{ij} | \xi_i, \eta_j &\sim \mathcal{B}(w(\xi_i, \eta_j)), \tag{16} 
\end{align} 
$$

$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 [19] 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 (14). For the LBM, $w(\xi_i, \eta_j) = \sum_{k, \ell} \pi_{k\ell} \mathbb{1}\{s(\xi_i) = k\} \mathbb{1}\{t(\eta_j) = \ell\}$ where $s(\xi_i) = 1 + \sum_{k=1}^K \mathbb{1}\{\xi_i > \sum_{k'=1}^k \alpha_{k'}\}$, $t(\eta_j) = 1 + \sum_{\ell=1}^L \mathbb{1}\{\eta_j > \sum_{\ell'=1}^\ell \beta_\ell\}$ and $\alpha, \beta$ and $\pi$ are those of Equation (15).

**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 suitable for count values in $\mathbb{N}$:

$$ 
\begin{align} 
\xi_i, \eta_j &\stackrel{\text{i.i.d.}}{\sim} \mathcal{U}[0, 1] \nonumber \\ 
Y_{ij} | \xi_i, \eta_j &\sim \mathcal{P}(w(\xi_i, \eta_j)), \tag{17} 
\end{align} 
$$

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

## 6.3. Examples of kernel functions for network analysis

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$.
<figure><img src="image_2.png" /><figcaption>Fig 1. Motif 6 and motif 14, counted by $U_N^{h_6}$ and $U_N^{h_{14}}$.</figcaption></figure>

**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 [56, 33]. Many random network models hinge on motif frequencies. In the Exponential Random Graph Model [24], motif frequencies are sufficient statistics. The dk-random graph model [46] also largely relies on motif frequencies.

Their asymptotic properties are widely studied and a large number of studies use motif counts to perform statistical tests [50, 10, 9, 17, 25, 39, 42, 47]. Our framework is particularly well adapted to the use of motif counts for statistical tests because frequencies are U-statistics with kernel functions of the same size as the motifs.

In many applications, motifs are viewed as elementary building blocks that can be specifically interpreted. This is the case in molecular biology [56, 48, 1], neurology [66], sociology [6], evolution [49] and ecology [61, 4]. In particular, Figure 3 of [58] lists all the bipartite motifs consisting of from 2 to 6 nodes. For example, their motif 6 represents a 2 × 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. [35] 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.

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

$$ 
\begin{aligned}
h_6^0(Y_{(i_1, i_2; j_1, j_2)}) &= Y_{i_1 j_1} Y_{i_1 j_2} Y_{i_2 j_1} Y_{i_2 j_2}, \\
h_{14}^0(Y_{(i_1, i_2; j_1, j_2, j_3)}) &= 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}).
\end{aligned} 
\quad (18) 
$$

**Remark 10.** 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_6^0$ is already symmetric, so $h_6 = h_6^0$. However, $h_{14}^0$ 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{aligned}
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}) \\
&\quad + \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}) \\
&\quad + \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}) \\
&\quad + \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})
\end{aligned} 
$$
$$ + \frac{1}{6}Y_{i_2j_2}Y_{i_2j_3}Y_{i_1j_3}Y_{i_1j_1}(1-Y_{i_2j_1})(1-Y_{i_1j_2}) + \frac{1}{6}Y_{i_2j_3}Y_{i_2j_1}Y_{i_1j_1}Y_{i_1j_2}(1-Y_{i_2j_2})(1-Y_{i_1j_3}). $$

Using Proposition 5.4, the following studentized statistics converge to a standard normal distribution

$$ Z_N^6 := \sqrt{\frac{N}{\hat{V}_N^{h_6}}} (U_N^{h_6} - U_\infty^{h_6}) \xrightarrow{D} \mathcal{N}(0, 1), $$

$$ Z_N^{14} := \sqrt{\frac{N}{\hat{V}_N^{h_{14}}}} (U_N^{h_{14}} - U_\infty^{h_{14}}) \xrightarrow{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 (17) such that

$$ Y_{ij} | \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.

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

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}) := ||\bar{w} - fg||_2^2 = \iint (\bar{w}(\xi, \eta) - f(\xi)g(\eta))^2 d\xi d\eta. \qquad (20) $$

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 (6.3) is a BEDD model, and $f$ and $g$ have the same role as in (14).

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 C.3 applies to the associated U-statistics $(U_N^{h_A}, U_N^{h_B}, U_N^{h_C}, U_N^{h_D})$, and they are jointly asymptotically normal with asymptotic covariance matrix $\Sigma^{h_A,h_B,h_C,h_D}$.

We define the estimator $\hat{d}_N$ of $d(\bar{w})$ as

$$ \hat{d}_N = t(U_N^{h_A}, U_N^{h_B}, U_N^{h_C}, U_N^{h_D}) := U_N^{h_A} / (U_N^{h_D})^3 + U_N^{h_B} U_N^{h_C} / (U_N^{h_D})^4. $$

We have

$$ \begin{aligned} t(U_{\infty}^{h_A}, U_{\infty}^{h_B}, U_{\infty}^{h_C}, U_{\infty}^{h_D}) &= \iint \bar{w}(\xi, \eta)(\bar{w}(\xi, \eta) - 2f(\xi)g(\eta))d\xi d\eta \\ &\quad + \int f(\xi)^2 d\xi \times \int g(\eta)^2 d\eta \\ &= d(\bar{w}). \end{aligned} $$
<table><caption>Table 1
Kernel functions (before symmetrization) used for the estimation of d(\bar{w}). The expectations are given by Lemma E.1 in Appendix.</caption><thead><tr><th></th><th>Kernel</th><th>Expectation</th></tr></thead><tbody><tr><td>h<sub>A</sub></td><td>h<sup>0</sup><sub>A,1</sub> - 2h<sup>0</sup><sub>A,2</sub></td><td>λ<sup>3</sup> ∬ w&#773;(ξ, η)(w&#773;(ξ, η) - 2f(ξ)g(η))dξdη</td></tr><tr><td>h<sub>A,1</sub></td><td>Y<sub>i<sub>1</sub>j<sub>1</sub></sub>(Y<sub>i<sub>1</sub>j<sub>1</sub></sub> - 1)Y<sub>i<sub>2</sub>j<sub>2</sub></sub></td><td>λ<sup>3</sup> ∬ w&#773;(ξ, η)<sup>2</sup>dξdη</td></tr><tr><td>h<sub>A,2</sub></td><td>Y<sub>i<sub>1</sub>j<sub>1</sub></sub>Y<sub>i<sub>1</sub>j<sub>2</sub></sub>Y<sub>i<sub>2</sub>j<sub>2</sub></sub></td><td>λ<sup>3</sup> ∬ w&#773;(ξ, η)f(ξ)g(η)dξdη</td></tr><tr><td>h<sub>B</sub></td><td>Y<sub>i<sub>1</sub>j<sub>1</sub></sub>Y<sub>i<sub>1</sub>j<sub>2</sub></sub></td><td>λ<sup>2</sup> ∫ f(ξ)<sup>2</sup>dξ</td></tr><tr><td>h<sub>C</sub></td><td>Y<sub>i<sub>1</sub>j<sub>1</sub></sub>Y<sub>i<sub>2</sub>j<sub>1</sub></sub></td><td>λ<sup>2</sup> ∫ g(η)<sup>2</sup>dη</td></tr><tr><td>h<sub>D</sub></td><td>Y<sub>i<sub>1</sub>j<sub>1</sub></sub></td><td>λ</td></tr></tbody></table>

Then, Proposition 5.7 ensures that the studentized statistic converges to a standard normal distribution

$$Z_N^d := \sqrt{\frac{N}{\hat{V}_N^d}} \left( \hat{d}_N - d(\bar{w}) \right) \xrightarrow{D} \mathcal{N}(0, 1),$$

where

$$\hat{V}_N^d = \nabla_t (U_N^{h_A}, U_N^{h_B}, U_N^{h_C}, U_N^{h_D})^T \hat{\Sigma}_N^{h_A, h_B, h_C, h_D} \nabla_t (U_N^{h_A}, U_N^{h_B}, U_N^{h_C}, U_N^{h_D}),$$

and $\hat{\Sigma}_N^{h_A, h_B, h_C, h_D}$ is defined as in Proposition 5.7.

**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 | \xi_i] = f(\xi_i)\mathbb{E}[D_i]$, where f is given by (19). Therefore, the marginals f and g account 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 node 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 does not depend entirely on the number of edges. In a network with weighted edges, $f(\xi_i)$ corresponds to the 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 [59]. 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 $\hat{F}_{2,N} := \kappa(U_N^{h_1}, U_N^{h_2}) = U_N^{h_1}/U_N^{h_2}$, 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_{\infty}^{h_1} = \lambda^2 F_2$ and $U_{\infty}^{h_2} = \lambda^2$. Given that $\nabla_{\kappa}(U_N^{h_1}, U_N^{h_2}) = (1/U_N^{h_2}, -U_N^{h_1}/(U_N^{h_2})^2)$, the application of Proposition 5.7 yields

$$Z_N^{F_2} := \sqrt{\frac{N}{\hat{V}_N^{F_2}}} (\hat{F}_{2,N} - F_2) \xrightarrow{D} \mathcal{N}(0, 1), \quad (21)$$

where

$$\hat{V}_N^{F_2} = \frac{1}{(U_N^{h_2})^2} \hat{V}_N^{h_1} - \frac{2\hat{F}_{2,N}}{(U_N^{h_2})^2} \hat{C}_N^{h_1,h_2} + \frac{(\hat{F}_{2,N})^2}{(U_N^{h_2})^2} \hat{V}_N^{h_2}.$$

In a similar way, one can define the analogous estimator $\hat{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.

# 7. Simulations

The methodology proposed in the previous section relies on Propositions 5.4 and 5.7, which are asymptotic results. In this section, we illustrate the behavior of the statistics defined in the previous section on synthetic networks, investigating the effect of network size. For each example, 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. Finally, we illustrate the computational efficiency of our variance estimator, compared to naive approaches.

## 7.1. Motif counts

The following simulation results illustrate the application to network motif counting developed in Section 6.3.

**Model I:** This graph model is an LBM with two row groups and two column groups of equal proportion. Using the notations of (15), this means $\alpha = \beta = (0.5, 0.5)$. The probability matrix $\pi$ has size $2 \times 2$. We set $\pi_{k\ell} = 0.5$ for all $1 \le k, \ell \le 2$ except $\pi_{11} = 0.95$.

Under Model I, for each value $N \in \{2^k : k \in \llbracket 5, 12 \rrbracket\}$ and $\rho \in \{1/8, 1/2\}$, we 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 [58] are respectively given by $U_N^{h_6}$ and $U_N^{h_{14}}$ where $h_6$ and $h_{14}$ were defined by (18).

The Q-Q plots of the studentized statistic associated to $h_6$ and $h_{14}$ are given in Figure 2 and 3. For $\rho = 1/2$, we observe that the empirical distributions of both statistics converge and become close to a normal distribution as $N \ge 128$. Figure 4 gives the respective coverage prob-

$$\text{abilities for } U_{\infty}^{h_6} \text{ and } U_{\infty}^{h_{14}} \text{ of the intervals } \left[ U_N^{h_6} - \Phi(1 - \frac{\alpha}{2}) \sqrt{\frac{N}{\hat{V}_N^{h_6}}}, U_N^{h_6} + \Phi(1 - \frac{\alpha}{2}) \sqrt{\frac{N}{\hat{V}_N^{h_6}}} \right]$$
<figure><img src="image_3.png" /><figcaption>Fig 2. Q-Q plots for $Z_N^6$ the studentized statistic associated with $U_N^{h_6}$, with $\rho = 0.5$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig 3. Q-Q plots for $Z_N^{14}$ the studentized statistic associated with $U_N^{h_{14}}$, with $\rho = 0.5$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig 4. Empirical coverage probabilities for the asymptotic confidence intervals at level $\alpha = 0.95$ of $U_N^{h_6}$ (left) and $U_N^{h_{14}}$ (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$.</figcaption></figure>

and $\left[ U_N^{h_{14}} - \Phi(1 - \frac{\alpha}{2}) \sqrt{\frac{N}{\hat{V}_N^{h_{14}}}} \right], U_N^{h_{14}} + \Phi(1 - \frac{\alpha}{2}) \sqrt{\frac{N}{\hat{V}_N^{h_{14}}}}$ 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><img src="image_6.png" /><figcaption>Fig 5. Q-Q plot for $Z_N^d$, the studentized statistic associated with $\hat{d}_N$, $\epsilon = 3$, $\rho = 0.5$.</figcaption></figure>

## 7.2. Graphon product distance

The following simulation results illustrate the application to graphon analysis developed in Section 6.3.

**Model II(ε):** We consider a Poisson-LBM. As this is a weighted graph model, $\pi$ is a matrix of weights rather than probabilities, i.e. $(\pi_{k\ell})_{1\le k\le K, 1\le\ell\le L}$ are real non-negative numbers. We consider two row groups and two column groups of equal proportion. Let $\pi^0$ be a weight matrix. We set $\pi_{11}^0 = 4$, $\pi_{12}^0 = \pi_{21}^0 = 2$ and $\pi_{22}^0 = 1$. The LBM with parameters $\alpha, \beta$ and $\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 $\tau$ be the 2 × 2 matrix where $\tau_{11} = \tau_{22} = 2$ and $\tau_{12} = \tau_{21} = 0$. For $\epsilon \ge 0$, we define Model II(ε) as a LBM with group probabilities $\alpha$ and $\beta$ and with weight matrix $\pi^\epsilon = \frac{\lambda}{\lambda+\epsilon}(\pi^0 + \epsilon\tau)$. Thus, Model II(0) is an LBM with product form (BEDD model) and for $\epsilon > 0$, Model II(ε) 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(ε) 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 \ge 0$.

Under Model II(ε), for each value $N \in \{2^k : k \in \{\lfloor 5, 12 \rfloor\}\}$, $\rho \in \{1/8, 1/2\}$, $\epsilon \in \{0.5, 1, 1.5, 2, 2.5, 3\}$, we 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 computed $\hat{d}_N$ and $\hat{V}_N^d$ as estimates of $d(\bar{w}_\epsilon)$ and $V^d$ respectively. The Q-Q plots of the studentized statistic $Z_N^d$ are given in Figure 5. It is apparent that $Z_N^d$ is not centered, especially for $N \lesssim 512$. This is due to the fact that $\hat{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 \ge 2048$. Figure 6 gives for different values of $\epsilon$, the coverage probability for $d(\bar{w}_\epsilon)$ of the interval $[\hat{d}_N - \Phi(1 - \frac{\alpha}{2})\sqrt{\frac{N}{\hat{V}_N^d}}, \hat{d}_N + \Phi(1 - \frac{\alpha}{2})\sqrt{\frac{N}{\hat{V}_N^d}}]$ 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><img src="image_7.png" /><figcaption>Fig 6. Empirical coverage probabilities for the asymptotic confidence intervals at level $\alpha = 0.95$ of $\hat{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$.</figcaption></figure>

### 7.3. Heterogeneity in the row weights of a network

The following simulation results illustrate the application to network row heterogeneity analysis developed in Section 6.3.

**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)$. We set $\lambda = 1$, $F_2 = 3$ and $G_2 = 2$.

Under Model III, for each value $N \in \{2^k : k \in \llbracket 5, 12 \rrbracket\}$ and $\rho \in \{1/8, 1/2, 7/8\}$, we 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 estimated $F_2$ and $G_2$ using $\hat{F}_{2,N}$ and its symmetric counterpart $\hat{G}_{2,N}$. We also computed $\hat{V}_N^{F_2}$ and $\hat{V}_N^{G_2}$ to obtain the studentized statistics $Z_N^{F_2}$ and $Z_N^{G_2}$. Q-Q plots of $Z_N^{F_2}$ are given in Figure 7. Figure 8 gives the respective coverage probabilities for $F_2$ and $G_2$ of the intervals

$$ \left[ \hat{F}_{2,N} - \Phi\left(1 - \frac{\alpha}{2}\right) \sqrt{\frac{N}{\hat{V}_N^{F_2}}}, \hat{F}_{2,N} + \Phi\left(1 - \frac{\alpha}{2}\right) \sqrt{\frac{N}{\hat{V}_N^{F_2}}} \right] $$

and

$$ \left[ \hat{G}_{2,N} - \Phi\left(1 - \frac{\alpha}{2}\right) \sqrt{\frac{N}{\hat{V}_N^{G_2}}}, \hat{G}_{2,N} + \Phi\left(1 - \frac{\alpha}{2}\right) \sqrt{\frac{N}{\hat{V}_N^{G_2}}} \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 corresponds to $\rho = 1/8$ (there are fewer rows) and that for $G_2$, the worst case corresponds to $\rho = 7/8$ (there are fewer columns), but Figure 8 does not show a clear difference between the three values of $\rho$.
<figure><img src="image_8.png" /><figcaption>Fig 7. Q-Q plot for $Z_N^{F_2}$ the studentized statistic associated with $\hat{F}_{2,N}$, $\rho = 0.5$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig 8. Empirical coverage probabilities for the asymptotic confidence intervals at level $\alpha = 0.95$ of $Z_N^{F_2}$ (left) and $Z_N^{G_2}$ (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$.</figcaption></figure>

## 7.4. Computational efficiency of variance estimators

In the following simulation study, we assess the computational performance of variance estimation using the techniques developed in this paper. Specifically, we consider the $U$-statistics $U_N^{h_1}$ and $U_N^{h_2}$ defined in Section 6.3. To estimate their asymptotic variances $V^{h_i}$ for $i \in \{1, 2\}$, we compare three algorithms:

* Algorithm A: We estimate $V^{h_i}$ using empirical covariance estimators for $v_{h_i}^{1,0} = \gamma_h^{1,0}$ and $v_{h_i}^{0,1} = \gamma_h^{0,1}$, as described in the beginning of Section 5. This is the direct approach to estimating the asymptotic variance of a $U$-statistic on RCE matrices.

* Algorithm B: We compute $\hat{V}_N^{h_i}$ defined in Section 5.3 by directly applying the kernel function to all subgraphs in the summations of equations (10) and (11).

* Algorithm C: We compute $\hat{V}_N^{h_i}$ defined in Section 5.3, using an alternative approach that evaluates $U$-statistics on smaller networks and then combines them using the expressions of (12) and (13). This method leverages matrix operations for efficient computation:

$$U_N^{h_1} = \frac{1}{n_N m_N (m_N - 1)} \left[ |Y_N^T Y_N|_1 - \text{Tr}(Y_N^T Y_N) \right], \quad (22)$$

$$U_N^{h_2} = \frac{1}{m_N(m_N - 1)n_N(n_N - 1)} \left[ (|Y_N|_1)^2 - |Y_N^T Y_N|_1 + \text{Tr}(Y_N^T Y_N) \right] \quad (23)$$
<figure><img src="image_10.png" /><figcaption>Fig 9. Empirical average of the estimates for $V^{h_1}$ (left) and $V^{h_2}$ (right) given by Algorithms A, B and C for different network sizes N (in log-scale). The dashed lines represent the true values of $V^{h_1}$ and $V^{h_2}$. Algorithm A was not tested for $N \ge 32$ due to its excessively high computation time.</figcaption></figure>

$$-|Y_N Y_N^T|_1 + \text{Tr}(Y_N Y_N^T) - |Y_N^{\odot 2}|_1],$$

where $Y_N$ is the observed network, with $m_N$ rows and $n_N$ columns.

For each value $N \in \{2^{k/2} : k \in \llbracket 6, 16 \rrbracket\}$, we simulated $K = 100$ networks of size $m_N \times n_N$ where $m_N = \lfloor N/2 \rfloor$ and $n_N = N - m_N$, under Model III. We then estimated $V^{h_1}$ and $V^{h_2}$ using Algorithms A, B and C. Due to its excessive computational cost, Algorithm A has been used only for $N < 32$. Appendix F presents the average computation times for each variance estimate. They indicate that, while Algorithm A is by far the slowest among the three, Algorithm C is significantly faster than Algorithm B due to the matrix operation trick, despite both computing the same estimator.

Figure 9 presents the average values of the estimates across network sizes. They show quick convergence for the estimates $\hat{V}_N^{h_1}$ and $\hat{V}_N^{h_2}$ provided by Algorithms B and C. In contrast, Algorithm A not only requires a higher computational cost but also produces substantially less precise estimates over the considered network size range. Furthermore, Algorithm A occasionally yields negative variance estimates, which is highly undesirable. These results highlight the advantages of the proposed approach compared to the classic method.

## 8. Illustrations

To illustrate our methodology and the interpretation of some of the U-statistics introduced in this paper, we considered the set of law-makers networks compiled by [41]. 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 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 10 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><img src="image_11.png" /><figcaption>Fig 10. Distribution of the number of members (left), number of laws (center), and density (right) of the four lawmakers networks across the 26 years (in log-scale).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig 11. Evolution of degree imbalance in each assembly over 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$).</figcaption></figure>

**Degree imbalance** We then focused on the degree of imbalance among the rows (resp. columns), which, under the weighted BEDD model defined in Equation (14), can be measured by the *U*-statistic $F_2$ (resp. $G_2$). Figure 11 gives the evolution of each of the two indicators over 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 exists, neither among the members ($F_2$ > 1), nor among the laws ($G_2$ > 1). Regarding the US networks (USH and USS), the imbalance is 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).

**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 6.1.
<figure><img src="image_13.png" /><figcaption>Fig 12. Network comparison for the U-statistics, $F_2$, $G_2$ and for the count the motif 6 (M6). EP-UN= —□—, EP-USS= —◇—, EP-USH= —△—, UN-USS= —◇—, UN-USH= —△—, USS-USH= —△—. Horizontal lines = standard normal quantiles with levels .025 and .975.</figcaption></figure>

Figure 12 displays 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 members 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. In 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.

## 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 [20], 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 \ge 1}$, i.e. such that for all $n \ge 1$, $\mathcal{F}_{n+1} \subset \mathcal{F}_n$.

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

**Theorem A.3.** Let $(M_n, \mathcal{F}_n)_{n \ge 1}$ be a backward martingale. Then, $(M_n)_{n \ge 1}$ is uniformly integrable, and, denoting $M_\infty = \mathbb{E}[M_1 | \mathcal{F}_\infty]$ where $\mathcal{F}_\infty = \bigcap_{n=1}^\infty \mathcal{F}_n$, we have

$$M_n \xrightarrow[n \to \infty]{a.s., L_1} M_\infty.$$

## Appendix B: Additional results presented in Section 3 and proofs

### B.1. Orthogonality of projection spaces

First, we prove Proposition 3.1. This proposition relies on the fact that the projections are conditionally centered. This property is given by the following lemma, the proof of which is
provided in Appendix B.

**Lemma B.1.** 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) \le (p, q)$. For all $\mathbf{i} \subset \mathbf{i}$ and $\mathbf{j} \subset \mathbf{j}$, we have

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

*Proof.* 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},\mathbf{0}}) | \mathcal{A}_{\mathbf{0},\mathbf{0}}] = \mathbb{E}[\psi^{1,0}h(Y_{\mathbf{i},\mathbf{0}})] - \mathbb{E}[h(Y_{\llap{[}p\llap{]},\llap{[}q\llap{]}})] = 0 $$

and

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

Suppose that the lemma is true for all $(0, 0) < (r', c') < (r, c)$. Let $(\mathbf{i}, \mathbf{j}) \in \mathcal{P}_r(\mathbb{N}) \times \mathcal{P}_c(\mathbb{N})$, $\mathbf{i} \subset \mathbf{i}$ and $\mathbf{j} \subset \mathbf{j}$. Denote $\underline{\mathcal{R}} = \text{Card}(\mathbf{i})$ and $\underline{\mathcal{C}} = \text{Card}(\mathbf{j})$. We can write

$$
\begin{align*}
\mathbb{E}[p^{r,c}h(Y_{\mathbf{i},\mathbf{j}}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] &= \mathbb{E}[\psi^{r,c}h(Y_{\mathbf{i},\mathbf{j}}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] - \mathbb{E}[p^{\underline{r},\underline{c}}h(Y_{\mathbf{i},\mathbf{j}}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] \\
&\quad - \sum_{(0,0)<(r',c')<(r,c)} \sum_{\substack{i' \in \mathcal{P}_{r'}(\mathbf{i}), j' \in \mathcal{P}_{c'}(\mathbf{j}) \\ (i',j') \neq (\mathbf{i},\mathbf{j})}} \mathbb{E}[p^{r',c'}h(Y_{i',j'}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] \\
&= \sum_{(0,0) \le (r',c') < (\underline{r},\underline{c})} \sum_{\substack{i' \in \mathcal{P}_{r'}(\mathbf{i}), j' \in \mathcal{P}_{c'}(\mathbf{j}) \\ (i',j') \neq (\mathbf{i},\mathbf{j})}} \mathbb{E}[p^{r',c'}h(Y_{i',j'}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] \\
&\quad - \sum_{(0,0)<(r',c')<(r,c)} \sum_{\substack{i' \in \mathcal{P}_{r'}(\mathbf{i}), j' \in \mathcal{P}_{c'}(\mathbf{j}) \\ (i',j') \neq (\mathbf{i},\mathbf{j})}} \mathbb{E}[p^{r',c'}h(Y_{i',j'}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] \\
&= - \sum_{(0,0)<(r',c')<(r,c)} \sum_{\substack{i' \in \mathcal{P}_{r'}(\mathbf{i}), j' \in \mathcal{P}_{c'}(\mathbf{j}) \\ i' \notin \mathbf{i}, j' \notin \mathbf{j}}} \mathbb{E}[p^{r',c'}h(Y_{i',j'}) | \mathcal{A}_{\mathbf{i}' \cap \mathbf{i}, \mathbf{j}' \cap \mathbf{j}}].
\end{align*}
$$

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

Since the last sum excludes the case $i' = \underline{i}$ and $j' = \underline{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}}) | \mathcal{A}_{\mathbf{i},\mathbf{j}}] = 0$, which concludes the proof by induction. □

*Proof of Proposition 3.1.* The two properties are proven similarly and derive from the fact that $(\mathbf{i}_1, \mathbf{j}_1) \ne (\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) \ne (\mathbf{i}_2, \mathbf{j}_2)$. Then $\mathbf{i}_1 \ne \mathbf{i}_2$ or $\mathbf{j}_1 \ne \mathbf{j}_2$. Without loss of generality, assume that $\mathbf{i}_1 \ne \mathbf{i}_2$ so there is an element $i_2 \in \mathbf{i}_2$ which is not included in $\mathbf{i}_1$. Then

$$
\begin{align*}
\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}) | \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}) | \mathcal{A}_{\{i_2\}, \emptyset}]] \\
&= 0,
\end{align*}
$$

since $\mathbb{E}[p^{r_2, c_2} h_2(Y_{\mathbf{i}_2, \mathbf{j}_2}) | \mathcal{A}_{\{i_2\}, \emptyset}] = 0$ from Lemma B.1.

We have proven both properties, since the case $(r_1, c_1) \ne (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. □
## B.2. Variance decomposition

Here, we state and prove Corollary B.2, of which Corollary 3.2, providing a variance decomposition for U-statistics, is a direct consequence.

**Corollary B.2.** 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) \le (r_1, c_1) < (p_1, q_1)$ and $(0,0) \le (r_2, c_2) < (p_2, q_2)$.

1. If $(r_1, c_1) \ne (r_2, c_2)$, then

$$ \mathrm{Cov}(P_{m,n}^{r_1, c_1} h_1(Y), P_{m,n}^{r_2, c_2} h_2(Y)) = 0. $$

2. If $(r_1, c_1) = (r_2, c_2) = (r, c)$, then

$$ \begin{aligned} & \mathrm{Cov}(P_{m,n}^{r,c} h_1(Y), P_{m,n}^{r,c} 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)). \end{aligned} $$

*Proof.* First, we see that for some $(r_1, c_1)$ and $(r_2, c_2)$,

$$ \mathrm{Cov}(P_{m,n}^{r_1, c_1} h_1(Y), P_{m,n}^{r_2, c_2} h_2(Y)) = \left[ \binom{m}{r_1} \binom{m}{r_2} \binom{n}{c_1} \binom{n}{c_2} \right]^{-1} \sum_{\mathbf{i}_1 \in \mathcal{P}_{r_1}(\llbracket m \rrbracket)} \sum_{\mathbf{i}_2 \in \mathcal{P}_{r_2}(\llbracket m \rrbracket)} \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) \ne (r_2, c_2)$, then from Proposition 3.1, all the covariance 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_{m,n}^{r_1, c_1} h_1(Y), P_{m,n}^{r_2, c_2} 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 3.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) \ne (\mathbf{i}_2, \mathbf{j}_2)$. Using this fact and the exchangeability of $Y$, we have

$$ \begin{aligned} & \mathrm{Cov}(P_{m,n}^{r,c} h_1(Y), P_{m,n}^{r,c} h_2(Y)) \\ &= \left[ \binom{m}{r} \binom{n}{c} \right]^{-2} \sum_{\mathbf{i}_1 \in \mathcal{P}_r(\llbracket m \rrbracket)} \sum_{\mathbf{i}_2 \in \mathcal{P}_r(\llbracket m \rrbracket)} \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_{\substack{\mathbf{i}_1 \in \mathcal{P}_r(\llbracket m \rrbracket) \\ \mathbf{j}_1 \in \mathcal{P}_c(\llbracket n \rrbracket)}} \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]^{-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{aligned} $$

which proves the second part of the corollary.

□
# Appendix C: Additional results presented in Section 4 and proofs

## C.1. Proofs for the asymptotic normality of U-statistics

*Proof of Lemma 4.2.* We only prove the first convergence result, since the second can be deduced by analogy.

By definition,

$$p_{(\{i\}, \emptyset)}^{1,0} h = \psi_{(\{i\}, \emptyset)}^{1,0} h - \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)]$$

First, we see that the $p_{(\{i\}, \emptyset)}^{1,0}h(Y)$ only depends on $\xi_i$, so all the $p_{(\{i\}, \emptyset)}^{1,0}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_{(\{i\}, \emptyset)}^{1,0}h] = \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)]$ so

$$\mathbb{E}[p_{(\{i\}, \emptyset)}^{1,0}h] = 0.$$

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

$$\begin{aligned}
\mathbb{V}[p_{(\{i\}, \emptyset)}^{1,0}h] &= v_h^{1,0} \\
&= \mathbb{V}[\mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket) | \xi_i]] \\
&= \mathrm{Cov}(h(Y_{(1,...,p;1,...,q)}), h(Y_{(1,p+1,...,2p-1;1,...,q)})) \\
&< \mathbb{V}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)] \\
&< \infty.
\end{aligned}$$

The classical CLT gives the desired result. □

*Proof of Lemma 4.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_{(\mathbf{i},\mathbf{j})}^{r,c}h] = \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)]$.

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

$$\begin{aligned}
\mathbb{E}[P_N^{r,c}h] &= \mathbb{E}[p_{(\mathbf{i},\mathbf{j})}^{r,c}h] \\
&= \mathbb{E}[\psi_{(\mathbf{i},\mathbf{j})}^{r,c}h] - \sum_{(0,0) \le (r',c') < (r,c)} \sum_{\substack{i' \in \mathcal{P}_{r'}(\mathbf{i}) \\ j' \in \mathcal{P}_{c'}(\mathbf{j})}} \mathbb{E}[p_{(\mathbf{i}',\mathbf{j}')}^{r',c'}h].
\end{aligned}$$

Now, Corollary 3.2 imply that $\mathrm{Cov}(P_N^{r,c}h, P_N^{r',c'}h) = 0$ unless $(r, c) = (r', c')$. So

$$\begin{aligned}
\mathbb{V}[A_N] &= N \sum_{\substack{(0,0) < (r,c) \le (p,q) \\ (r,c) \ne (1,0) \ne (0,1)}} \sum_{\substack{(0,0) < (r',c') \le (p,q) \\ (r',c') \ne (1,0) \ne (0,1)}} \binom{p}{r} \binom{p}{r'} \binom{q}{c} \binom{q}{c'} \mathrm{Cov}(P_N^{r,c}h, P_N^{r',c'}h) \\
&= N \sum_{\substack{(0,0) < (r,c) \le (p,q) \\ (r,c) \ne (1,0) \ne (0,1)}} \binom{p}{r}^2 \binom{q}{c}^2 \mathbb{V}[P_N^{r,c}h] \tag{24} \\
&= N \sum_{\substack{(0,0) < (r,c) \le (p,q) \\ (r,c) \ne (1,0) \ne (0,1)}} \binom{p}{r}^2 \binom{q}{c}^2 \binom{m_N}{r}^{-1} \binom{n_N}{c}^{-1} \mathbb{V}[p_{(\llbracket r \rrbracket, \llbracket c \rrbracket)}^{r,c}h].
\end{aligned}$$
From equation 6, we see that $h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)$ is a linear combination of all the $p_{(i,j)}^{r,c}h$, for $0 \le r \le p, 0 \le c \le q, i \in \mathcal{P}_r(\llbracket p \rrbracket), j \in \mathcal{P}_c(\llbracket q \rrbracket)$. So $\mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)^2] < \infty$ ensures that $\mathbb{V}[p_{(r,c)}^{r,c}, \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 \to \infty]{P} 0$. □

### C.2. The asymptotic normality result for functions of U-statistics

In this section, we prove that functions of U-statistics are also asymptotically normal under non-degeneracy conditions. First, 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 kernel 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 C.1.** Let $h$ be a kernel function of size $p \times q$. Let $p' \ge p$ and $q' \ge q$. We define the extension of $h$ to the size $p' \times q'$ by $\tilde{h}$ such that for all $i' \in \mathcal{P}_{p'}(\mathbb{N})$ and $j' \in \mathcal{P}_{q'}(\mathbb{N})$,

$$ \tilde{h}(Y_{i',j'}) = \left[ \binom{p'}{p} \binom{q'}{q} \right]^{-1} \sum_{\substack{i \in \mathcal{P}_p(i') \\ j \in \mathcal{P}_q(j')}} h(Y_{i,j}). $$

The extension of a kernel shares some properties with its kernel, as shown in the following lemma.

**Lemma C.2.**

$$ 2. \mathbb{E}[\tilde{h}(Y_{\llbracket p' \rrbracket}, \llbracket q' \rrbracket)] = \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)], $$

$$ 3. \mathbb{E}[\tilde{h}^2(Y_{\llbracket p' \rrbracket}, \llbracket q' \rrbracket)] \le \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)^2]. $$

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

$$ \begin{align*} \mathbb{E}[\tilde{h}(Y_{\llbracket p' \rrbracket}, \llbracket q' \rrbracket)^2] & = \left[ \binom{p'}{p} \binom{q'}{q} \right]^{-2} \sum_{\substack{i_1 \in \llbracket p' \rrbracket \\ j_1 \in \llbracket q' \rrbracket}} \sum_{\substack{i_2 \in \llbracket p' \rrbracket \\ j_2 \in \llbracket q' \rrbracket}} \mathbb{E}[h(Y_{i_1, j_1}) h(Y_{i_2, j_2})] \\ & = \left[ \binom{p'}{p} \binom{q'}{q} \right]^{-2} \sum_{\substack{i_1 \in \llbracket p' \rrbracket \\ j_1 \in \llbracket q' \rrbracket}} \sum_{\substack{i_2 \in \llbracket p' \rrbracket \\ j_2 \in \llbracket q' \rrbracket}} \left( \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)^2] \right. \\ & \qquad \left. - \frac{1}{2} \mathbb{E}[(h(Y_{i_1, j_1}) - h(Y_{i_2, j_2}))^2] \right) \\ & \le \left[ \binom{p'}{p} \binom{q'}{q} \right]^{-2} \sum_{\substack{i_1 \in \llbracket p' \rrbracket \\ j_1 \in \llbracket q' \rrbracket}} \sum_{\substack{i_2 \in \llbracket p' \rrbracket \\ j_2 \in \llbracket q' \rrbracket}} \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)^2] \\ & \le \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)^2]. \end{align*} $$

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_{h_k, h_l}^{r,c} := \text{Cov}(\psi_{\llbracket r \rrbracket, \llbracket c \rrbracket}^{r,c} h_k, \psi_{\llbracket r \rrbracket, \llbracket c \rrbracket}^{r,c} h_l)$. □
**Corollary C.3.** 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. Theorem 4.1 applies for each kernel function, i.e. $\mathbb{E}[h_k(Y_{[p_k]}, Y_{[q_k]})^2] < \infty$ and $U_{\infty}^{h_k}$ and $V^{h_k}$ are as defined in Theorem 4.1 for each kernel $h_k$, $1 \le k \le D$,

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 \le k \le 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_N^{h_1} \\ U_N^{h_2} \\ \vdots \\ U_N^{h_D} \end{pmatrix} - \begin{pmatrix} U_{\infty}^{h_1} \\ U_{\infty}^{h_2} \\ \vdots \\ U_{\infty}^{h_D} \end{pmatrix} \right) \underset{N \to \infty}{\to} \mathcal{N}(0, \Sigma^{h_1, ..., h_D}), $$

with

$$ \Sigma^{h_1, ..., h_D} = (C^{h_k, h_l})_{1 \le k, l \le D}, $$

where $C^{h_k, h_l} = \frac{p^2}{\rho} c_{h_k, h_l}^{1, 0} + \frac{q^2}{1-\rho} c_{h_k, h_l}^{0, 1}$ for all $1 \le k, l \le D$ (and $C^{h_k, h_k} = V^{h_k}$).

**Proof.** First, we show that $\lim_{N \to +\infty} \text{NCov}(U_N^{h_k}, U_N^{h_l})$

Let $(Z^{h_k})_{1 \le k \le 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 \le k \le D$ and $\text{Cov}(Z^{h_k}, Z^{h_l}) = C^{h_k, h_l}$ for all $1 \le k \le D$ and $1 \le l \le D$.

Denote $p := \max_k(p_k)$ and $q := \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 \times q$.

First, assume that $t \ne (0, ..., 0)$. Then by hypothesis, $\tilde{h}_t \ne 0$, therefore Lemma C.2 implies that $\sum_{k=1}^D t_k U_N^{h_k} = \sum_{k=1}^D t_k U_N^{\tilde{h}_k} = U_N^{\tilde{h}_t}$, the U-statistic with quadruplet kernel $\tilde{h}_t$ (of size $p \times q$). Using Cauchy-Schwarz inequality and the fact that from Lemma C.2, $\mathbb{E}[\tilde{h}_k(Y_{[p_k]}, Y_{[q_k]})^2] \le \mathbb{E}[h_k(Y_{[p_k]}, Y_{[q_k]})^2] < \infty$ for all $1 \le k \le D$, we have furthermore

$$ \begin{aligned} & \mathbb{E}[\tilde{h}_t (Y_{[p_k]}, Y_{[q_k]})^2] \\ &= \sum_{k=1}^n t_k^2 \mathbb{E}[\tilde{h}_k (Y_{[p_k]}, Y_{[q_k]})^2] + 2 \sum_{1 \le k \ne l \le D} t_k t_l \mathbb{E}[\tilde{h}_k (Y_{[p_k]}, Y_{[q_k]}) \tilde{h}_l (Y_{[p_k]}, Y_{[q_k]})], \\ & \le \sum_{k=1}^n t_k^2 \mathbb{E}[\tilde{h}_k (Y_{[p_k]}, Y_{[q_k]})^2] + 2 \sum_{1 \le k \ne l \le D} t_k t_l \sqrt{\mathbb{E}[\tilde{h}_k (Y_{[p_k]}, Y_{[q_k]})^2] \mathbb{E}[\tilde{h}_l (Y_{[p_k]}, Y_{[q_k]})^2]}, \\ & < \infty. \end{aligned} $$

Therefore, Theorem 4.1 also applies for $U_N^{\tilde{h}_t}$ and $\sqrt{N}(U_N^{\tilde{h}_t} - U_{\infty}^{\tilde{h}_t}) \underset{N \to \infty}{\to} \mathcal{N}(0, V^{\tilde{h}_t})$, where

$$ \bullet U_{\infty}^{\tilde{h}_t} = \sum_{k=1}^{D} t_k U_{\infty}^{h_k}, $$

$$ \bullet V^{\tilde{h}_t} = t^T \sum^{h_1, ..., h_D} t. $$

The second point comes from the fact that

$$ V^{\tilde{h}_t} = \lim_{N \to +\infty} N \sum_{k=1}^{D} \sum_{l=1}^{D} t_k t_l \text{Cov}(U_N^{h_k}, U_N^{h_l}) $$
and by Corollary 3.2, $\lim_{N \to +\infty} N\text{Cov}(U_N^{h_k}, U_N^{h_l}) = \frac{p^2}{c}c_{h_k,h_l}^{1,0} + \frac{q^2}{1-c}c_{h_k,h_l}^{0,1} = \sum_{k\ell}^{h_1, \dots, h_D}$. Therefore, we can conclude that $\sqrt{N}\sum_{k=1}^D t_k(U_N^{h_k} - U_\infty^{h_k}) = \sqrt{N}(U_N^{\tilde{h}_t} - U_\infty^{\tilde{h}_t}) \xrightarrow{D} \sum_{k=1}^D t_k Z^{h_k}$.

Now assume that $t = (0, \dots, 0)$. Then $\tilde{h}_t = 0$ so $U_N^{\tilde{h}_t} = 0 = \sum_{k=1}^D t_k Z^{h_k}$. Therefore, $\sqrt{N}\sum_{k=1}^D t_k(U_N^{h_k} - U_\infty^{h_k}) = \sqrt{N}(U_N^{\tilde{h}_t} - U_\infty^{\tilde{h}_t}) \xrightarrow{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_N^{h_k} - U_\infty^{h_k}) \xrightarrow{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 [12]) which states that $\sqrt{N}(U_N^{h_k} - U_\infty^{h_k})_{1 \le k \le D}$ converges jointly in distribution to $(Z^{h_k})_{1 \le k \le D}$, which is a centered multivariate Gaussian with covariance matrix $\Sigma^{h_1, \dots, h_D}$. $\square$

Although, this corollary requires the kernel functions of $(h_1, h_2, \dots, h_D)$ to be linearly independent, the corresponding U-statistics are not independent random variables, even asymptotically, because $\Sigma^{h_1, \dots, h_D}$ is not a diagonal matrix. One consequence of this corollary is that Theorem 4.1 can be extended to differentiable functions of U-statistics.

**Corollary C.4.** Let $h_1, \dots, h_D$ be $D$ kernel functions such that Corollary C.3 applies and $\Sigma^{h_1, \dots, h_D} = (C^{h_k,h_l})_{1 \le k, l \le D}$. Denote $\theta = (U_\infty^{h_1}, \dots, U_\infty^{h_D})$. Let $g: \mathbb{R}^d \to \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, \dots, h_D} \nabla g(\theta)$. If $V^\delta > 0$, then

$$ \sqrt{N}(g(U_N^{h_1}, \dots, U_N^{h_D}) - g(\theta)) \xrightarrow{D} \mathcal{N}(0, V^\delta). $$

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

$$ g(U_N^{h_1}, \dots, U_N^{h_D}) - g(\theta) = \nabla g(\theta)^T \left( (U_N^{h_1}, \dots, U_N^{h_D}) - \theta \right) + o_P \left( \| (U_N^{h_1}, \dots, U_N^{h_D}) - \theta \| \right). $$

From Corollary C.3, $\sqrt{N}((U_N^{h_1}, \dots, U_N^{h_D}) - \theta)$ converges to a multivariate normal distribution with asymptotic covariance matrix $\Sigma^{h_1, \dots, h_D}$, so the delta method (see Theorem 3.1 of [63]) can be applied to prove this proposition. $\square$

# Appendix D: Additional results presented in Section 5 and proofs

## D.1. Some useful notations and results

Consider the sets $S_N^{p,q}$ and $S_{N, (\underline{i}, \underline{j})}^{p,q}$ defined in Section 5.1. It is obvious that

$$ \text{Card}(S_N^{p,q}) = \binom{m_N}{p} \binom{n_N}{q} \quad \text{and} \quad \text{Card}(S_{N, (\underline{i}, \underline{j})}^{p,q}) = \binom{m_N - \underline{p}}{p - \underline{p}} \binom{n_N - \underline{q}}{q - \underline{q}}. $$

Let $I_K = (\underline{i}_1, \dots, \underline{i}_K)$ be a $K$-uplet of sets of row indices and $J_K = (\underline{j}_1, \dots, \underline{j}_K)$ a $K$-uplet of subsets of column indices. Denote

$$ \mathcal{T}_{N, (I_K, J_K)}^{p,q} := S_{N, (\underline{i}_1, \underline{j}_1)}^{p,q} \times S_{N, (\underline{i}_2, \underline{j}_2)}^{p,q} \times \dots \times S_{N, (\underline{i}_K, \underline{j}_K)}^{p,q}. $$

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

$$ T_N^{p,q}(I_K, J_K) := \frac{1}{\prod_{k=1}^K \text{Card}(\mathcal{T}_{N, (I_K, J_K)}^{p,q})} \sum_{\mathcal{T}_{N, (I_K, J_K)}^{p,q}} \mathbb{E}[X_{\underline{i}_1, \underline{j}_1} X_{\underline{i}_2, \underline{j}_2} \dots X_{\underline{i}_K, \underline{j}_K}]. \quad (25) $$
As a remark, $T_N^{p,q}((\emptyset), (\emptyset)) = U_N^h$.

By exchangeability, the quantities $\mathbb{E}[X_{(i_1,j_1)} X_{(i_2,j_2)} \dots X_{(i_K,j_K)}]$ do not depend on the row indices that do not belong to any pairwise intersection of the $(\mathbf{i}_k)_{1 \le k \le K}$. The same holds for column indices that do not belong to any pairwise intersection of the $(\mathbf{j}_k)_{1 \le k \le K}$. Therefore, assuming $m_N \ge \text{Card}(\cup_{k=1}^K \mathbf{i}_k)$ and $n_N \ge \text{Card}(\cup_{k=1}^K \mathbf{j}_k)$, we can define

$$ \alpha(I_K, J_K) := \mathbb{E}[X_{\underline{i}_1 \underline{j}_1} X_{\underline{i}_2 \underline{j}_2} \dots X_{\underline{i}_K \underline{j}_K}] \quad (26) $$

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

The following lemma will be helpful in later proofs as it provides the asymptotic behavior of $T_N^{p,q}(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 D.1.** Let $I_K = (\underline{i}_1, \dots, \underline{i}_K)$ and $J_K = (\underline{j}_1, \dots, \underline{j}_K)$ be K-uplets of respectively row and column indices. Let $\alpha(I_K, J_K)$ defined by (26). We have

$$ T_N^{p,q}(I_K, J_K) = \alpha(I_K, J_K) + O(m_N^{-1} + n_N^{-1}). $$

*Proof.* Let:

* $\underline{p}_k := \text{Card}(\underline{i}_k)$ and $\underline{q}_k := \text{Card}(\underline{j}_k)$, for $1 \le k \le K$,

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

* $\underline{\bar{p}} := \text{Card}(\cup_{k=1}^K \underline{i}_k)$ and $\underline{\bar{q}} := \text{Card}(\cup_{k=1}^K \underline{j}_k)$.

The sum $T_N^{p,q}(I_K, J_K)$ defined by equation (25) is a sum over $T_N$ expectation terms where

$$ T_N := \text{Card}(\mathcal{T}_{N, (I_K, J_K)}^{p,q}) = \prod_{k=1}^{K} \text{Card}(S_{N, (\underline{i}_k, \underline{j}_k)}^{p,q}) = \prod_{k=1}^{K} \binom{m_N - \underline{p}_k}{p - \underline{p}_k} \binom{n_N - \underline{q}_k}{\underline{q} - \underline{q}_k}. \quad (27) $$

These expectation terms $\mathbb{E}[X_{i_1,j_1} \dots X_{i_K,j_K}]$ only depend on the number of times elements appear in pairwise intersections between the $(\mathbf{i}_k, \mathbf{j}_k)$, $1 \le k \le K$. In particular, we have

$$ \mathbb{E}[X_{i_1,j_1} \dots X_{i_K,j_K}] = \alpha(I_K, J_K) $$

when all the $\mathbf{i}_k$ contains the elements of $\underline{i}_k$ and all the other elements do not appear in any other $\underline{i}_{k'}$ (see equation 26). Denote $A_N$ the number of terms of $T_N^{p,q}(I_K, J_K)$ where $\mathbb{E}[X_{i_1,j_1} \dots X_{i_K,j_K}] = \alpha(I_K, J_K)$, then

$$ T_N^{p,q}(I_K, J_K) = \frac{A_N}{T_N} \alpha(I_K, J_K) + \frac{1}{T_N} \sum_{\substack{\mathcal{T}_{N, (I_K, J_K)}^{p,q} \\ \mathbb{E}[X_{i_1,j_1} \dots X_{i_K,j_K}] \neq \alpha(I_K, J_K)}} \mathbb{E}[X_{i_1,j_1} X_{i_2,j_2} \dots X_{i_K,j_K}]. \quad (28) $$

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

$$ 0 \le |\mathbb{E}[X_{i_1,j_1} \dots X_{i_K,j_K}]| \le \prod_{k=1}^{K} \mathbb{E}[|X_{i_k,j_k}|^K]^{\frac{1}{K}} = \mathbb{E}[|X_{[p]},[q]|^K]. $$
In particular, this holds for all $(T_N - A_N)$ terms of the remaining sum in Equation (28) when $\mathbb{E}[X_{i_1,j_1}...X_{i_K,j_K}] \neq \alpha(I_K, J_K)$, so

$$0 \leq |T_N^{p,q}(I_K, J_K) - \frac{A_N}{T_N}\alpha(I_K, J_K)| \leq \left(1 - \frac{A_N}{T_N}\right) \mathbb{E}[|X_{\llbracket p \rrbracket}, |\llbracket q \rrbracket||^K] \quad (29)$$

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. We find an expression for $T_N$. Since for all $1 \leq k \leq K$, we have

$$\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^{\underline{q}-\underline{q}_k}}{(p-\underline{p}_k)!(q-\underline{q}_k)!} (1 + O(m_N^{-1} + n_N^{-1})),$$

then we have from Equation (27)

$$T_N = \left( \prod_{k=1}^{K} (p - \underline{p}_k)!(q - \underline{q}_k)! \right)^{-1} m_N^{Kp-\underline{p}} n_N^{K\underline{q}-\underline{Q}} (1 + O(m_N^{-1} + n_N^{-1})). \quad (30)$$

2. We find an expression for $A_N$. The summation over $(\mathbf{i}_k, \mathbf{j}_k) \in S_{N, (\mathbf{i}_k, \mathbf{j}_k)}^{p,q}$ 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{i}_k$ and $\mathbf{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{i}_k$. The same follows for the pick of $\mathbf{j}_1$, so there are $\binom{m_N - \bar{p}}{p - \bar{p}_1} \binom{n_N - \bar{q}}{q - \bar{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{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'=1}^{k-1} (p - \underline{p}_{k'})}{p - \underline{p}_k} \binom{n_N - \bar{q} - \sum_{k'=1}^{k-1} (q - \underline{q}_{k'})}{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_{i_1,j_1}...X_{i_K,j_K}] = \alpha(I_K, J_K)$ is

$$A_N = \prod_{k=1}^{K} \binom{m_N - \bar{p} - \sum_{k'=1}^{k-1} (p - \underline{p}_{k'})}{p - \underline{p}_k} \binom{n_N - \bar{q} - \sum_{k'=1}^{k-1} (q - \underline{q}_{k'})}{q - \underline{q}_k}.$$

But we see that, for $1 \leq k \leq K$, we have

$$\binom{m_N - \bar{p} - \sum_{k'=1}^{k-1} (p - \underline{p}_{k'})}{p - \underline{p}_k} \binom{n_N - \bar{q} - \sum_{k'=1}^{k-1} (q - \underline{q}_{k'})}{q - \underline{q}_k} = \frac{m_N^{p-\underline{p}_k} n_N^{\underline{q}-\underline{q}_k}}{(p-\underline{p}_k)!(q-\underline{q}_k)!} (1 + O(m_N^{-1} + n_N^{-1})).$$
So we have

$$A_N = \left( \prod_{k=1}^{K} (p - \underline{p}_k)!(q - \underline{q}_k)! \right)^{-1} m_N^{Kp-P} n_N^{Kq-Q} \left( 1 + O(m_N^{-1} + n_N^{-1}) \right).$$

3. Now with the expressions of $T_N$ and $A_N$, we can deduce that

$$T_N - A_N = O \left( m_N^{Kp-P} n_N^{Kq-Q} (m_N^{-1} + n_N^{-1}) \right),$$

so

$$1 - \frac{A_N}{T_N} = \frac{\left( \prod_{k=1}^{K} (p - \underline{p}_k)!(q - \underline{q}_k)! \right) O \left( m_N^{Kp-P} n_N^{Kq-Q} (m_N^{-1} + n_N^{-1}) \right)}{m_N^{Kp-P} n_N^{Kq-Q} (1 + O(m_N^{-1} + n_N^{-1}))} \\ = O \left( m_N^{-1} + n_N^{-1} \right).$$

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

$$T_N^{p,q}(I_K, J_K) = (1 + O(m_N^{-1} + n_N^{-1})) \alpha(I_K, J_K) + O(m_N^{-1} + n_N^{-1}) \\ = \alpha(I_K, J_K) + O(m_N^{-1} + n_N^{-1}). \quad \square$$

## D.2. Proofs for the convergence of the conditional expectation estimators

Here, we provide the proofs for Propositions 5.1 and 5.2, establishing the unbiasedness and asymptotic consistency of the conditional expectation estimators defined by (8) and (9).

*Proof of Proposition 5.1.* The first result can be found directly by the definition of $\psi_{(\{i\}, \emptyset)}^{1,0} h$, since

$$\begin{aligned} \mathbb{E}[\hat{\mu}_N^{h,(i)} | \xi_i] &= \binom{m_N-1}{p-1}^{-1} \binom{n_N}{q}^{-1} \sum_{(i,j) \in S_{N, (\{i\}, \emptyset)}^{p,q}} \mathbb{E}[h(Y_{i,j} | \xi_i)] \\ &= \binom{m_N-1}{p-1}^{-1} \binom{n_N}{q}^{-1} \sum_{(i,j) \in S_{N, (\{i\}, \emptyset)}^{p,q}} \psi_{(\{i\}, \emptyset)}^{1,0} h \\ &= \psi_{(\{i\}, \emptyset)}^{1,0} h. \end{aligned}$$

The second result can be obtained analogously. □

*Proof of Proposition 5.2.* Let $N \in \mathbb{N}$ and $\mathcal{F}_N(Y) = \sigma((\hat{\mu}_K^{h,(i)}(Y))_{K \ge N})$. Let $\Phi_N \in \mathbb{S}_{m_N}^{(i)} \times \mathbb{S}_{n_N}$ where $\mathbb{S}_{m_N}^{(i)}$ is the group of permutations $\sigma_i$ of $[m_N]$ such that $\sigma_i(i) = i$. If $\Phi_N = (\sigma_i, \tau)$, denote $\Phi_N Y = (Y_{\sigma_i(k)\tau(j)})_{\substack{1 \le k \le m_N \\ 1 \le j \le n_N}}$.

First, we observe that $\hat{\mu}_N^{h,(i)}(Y) = \hat{\mu}_N^{h,(i)}(\Phi_N Y)$, so $\mathcal{F}_N(\Phi_N Y) = \mathcal{F}_N(Y)$. Therefore, by the exchangeability of $Y$, we have

$$Y | \mathcal{F}_N(Y) \stackrel{\Delta}{=} \Phi_N Y | \mathcal{F}_N(Y).$$

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

$$ 
\begin{aligned}
\mathbb{E}[\hat{\mu}_N^{h,(i)}(Y) | \mathcal{F}_{N+1}(Y)] &= \mathbb{E}[h(Y_{(i,j)}) | \mathcal{F}_{N+1}(Y)] \\
&= \mathbb{E}[\hat{\mu}_{N+1}^{h,(i)}(Y) | \mathcal{F}_{N+1}(Y)] \\
&= \hat{\mu}_{N+1}^{h,(i)}(Y).
\end{aligned}
$$

Therefore, $\hat{\mu}_N^{h,(i)}(Y)$ is a backward martingale with respect to $\mathcal{F}_N(Y)$ (see Appendix A). By Theorem A.3, we have that $\hat{\mu}_N^{h,(i)}(Y) \xrightarrow[N \to \infty]{a.s.,L_1} \mathbb{E}[\hat{\mu}_1^{h,(i)}(Y) | \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 5.1 implies $\mathbb{E}[\hat{\mu}_1^{h,(i)}(Y) | \mathcal{F}_\infty(Y)] = \psi_{(\{i\},\emptyset)}^{1,0}h(Y)$ and thus, $\hat{\mu}_N^{h,(i)}(Y) \xrightarrow[N \to \infty]{a.s.,L_1} \psi_{(\{i\},\emptyset)}^{1,0}h(Y)$. $\square$

### D.3. Results for the convergence of the variance estimators and proofs

In this section, we state and prove Propositions D.2 and D.3 ensuring that the variance estimators defined by (10) and (11) are asymptotically unbiased and have vanishing variances.

**Proposition D.2.** We have $\mathbb{E}[\hat{v}_N^{h;1,0}] = v_h^{1,0} + O(N^{-1})$ and $\mathbb{E}[\hat{v}_N^{h;0,1}] = v_h^{0,1} + O(N^{-1})$. As a consequence, $\hat{v}_N^{h;1,0}$ and $\hat{v}_N^{h;0,1}$ are asymptotically unbiased estimators for $v_h^{1,0}$ and $v_h^{0,1}$.

*Proof.* In this proof, for the estimator defined by equation (8), we write $\hat{\mu}_N^{(i)}$ instead of $\hat{\mu}_N^{h,(i)}$ to simplify the notation without ambiguity. Notice that

$$ \mathbb{E}[\hat{v}_N^{h;1,0}] = \frac{1}{2} \mathbb{E}[(\hat{\mu}_N^{(1)} - \hat{\mu}_N^{(2)})^2] = \mathbb{E}[(\hat{\mu}_N^{(1)})^2] - \mathbb{E}[\hat{\mu}_N^{(1)} \hat{\mu}_N^{(2)}]. \quad (31) $$

* First, we calculate:

$$ \mathbb{E}[\hat{\mu}_N^{(1)}(Y)^2] = T_N^{p,q}(I_2, J_2) $$

where $I_2 = \{(1), \{1\}\}$ and $J_2 = (\emptyset, \emptyset)$. Applying Lemma D.1, with $\underline{P} = \text{Card}(\{1\}) + \text{Card}(\{1\}) = 2$, $\underline{Q} = \text{Card}(\emptyset) + \text{Card}(\emptyset) = 0$ and

$$ 
\begin{aligned}
\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)}) | \xi_1]] \\
&= \mathbb{E}[\mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket) | \xi_1]^2] \\
&= \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0}h^2],
\end{aligned}
$$

we find

$$ \mathbb{E}[(\hat{\mu}_N^{(1)})^2] = \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0}h^2] + O(N^{-1}). \quad (32) $$

* Next, we calculate:

$$ \mathbb{E}[\hat{\mu}_N^{(1)} \hat{\mu}_N^{(2)}] = T_N^{p,q}(I'_2, J'_2) $$
where $I'_2 = \{(1), \{2\}\}$ and $J'_2 = (\emptyset, \emptyset)$. Applying Lemma D.1 with

$$
\begin{align*}
\alpha(I'_2, J'_2) &= \mathbb{E}[h(Y_{(1,3,...,p+1;1,...,q)})h(Y_{(2,p+2,...,2p;q+1,...,2q)})] \\
&= \mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket)]^2 \\
&= \mathbb{E}[\psi_{\{(1),\emptyset\}}^{1,0} h]^2,
\end{align*}
$$

we find

$$
\mathbb{E}[\hat{\mu}_N^{(1)} \hat{\mu}_N^{(2)}] = \mathbb{E}[\psi_{\{(1),\emptyset\}}^{1,0} h]^2 + O(N^{-1}). \quad (33)
$$

Finally, we can combine Equations (31), (32) and (33) to obtain

$$
\mathbb{E}[\hat{v}_N^{h;1,0}] = \mathbb{E}[\psi_{\{(1),\emptyset\}}^{1,0} h^2] - \mathbb{E}[\psi_{\{(1),\emptyset\}}^{1,0} h]^2 + O(N^{-1}) = v_h^{1,0} + O(N^{-1}),
$$

which proves the result.

**Proposition D.3.** We have $\mathbb{V}[\hat{v}_N^{h;1,0}] = O(N^{-1})$ and $\mathbb{V}[\hat{v}_N^{h;0,1}] = O(N^{-1})$.

*Proof.* In this proof, we write $\hat{\mu}_N^{(i)}$ instead of $\hat{\mu}_N^{h,(i)}$ to simplify the notation without ambiguity.
Notice that

$$
\begin{align}
\mathbb{E}[(\hat{v}_N^{h;1,0})^2] &= \left(m_N/2\right)^{-2} \sum_{1 \le i_1 < i_2 \le m_N} \sum_{1 \le i'_1 < i'_2 \le m_N} \mathbb{E} \left[ \frac{(\hat{\mu}_N^{(i_1)} - \hat{\mu}_N^{(i_2)})^2 (\hat{\mu}_N^{(i'_1)} - \hat{\mu}_N^{(i'_2)})^2}{4} \right] \nonumber \\
&= \frac{1}{4} \mathbb{E}[(\hat{\mu}_N^{(1)} - \hat{\mu}_N^{(2)})^2 (\hat{\mu}_N^{(3)} - \hat{\mu}_N^{(4)})^2] + O(N^{-1}) \nonumber \\
&= \mathbb{E}[(\hat{\mu}_N^{(1)})^2 (\hat{\mu}_N^{(2)})^2] - 2\mathbb{E}[(\hat{\mu}_N^{(1)})^2 \hat{\mu}_N^{(2)} \hat{\mu}_N^{(3)}] + \mathbb{E}[\hat{\mu}_N^{(1)} \hat{\mu}_N^{(2)} \hat{\mu}_N^{(3)} \hat{\mu}_N^{(4)}] \nonumber \\
&\quad + O(N^{-1}). \tag{34}
\end{align}
$$

Now, we calculate each of the three expectation terms in this equation.

• First, we calculate:

$$
\mathbb{E}[(\hat{\mu}_N^{(1)})^2 (\hat{\mu}_N^{(2)})^2] = T_N^{p,q}(I_4, J_4),
$$

where $I_4 = \{(1), \{1\}, \{2\}, \{2\}\}$ and $J_4 = (\emptyset, \emptyset, \emptyset, \emptyset)$. Applying Lemma D.1 with

$$
\begin{align*}
\alpha(I_4, J_4) \\
&= \mathbb{E}[h(Y_{(1,3,...,p+1;1,...,q)})h(Y_{(1,p+2,...,2p;q+1,...,2q)}) \\
&\quad h(Y_{(2,2p+1,...,3p-1;2q+1,...,3q)})h(Y_{(2,3p,...,4p-2;3q+1,...,4q)})] \\
&= \mathbb{E}[h(Y_{(1,...,p;1,...,q)})h(Y_{(1,p+1,...,2p-1;q+1,...,2q)})]^2 \\
&= \mathbb{E}[\mathbb{E}[h(Y_{(1,...,p;1,...,q)})h(Y_{(1,p+1,...,2p-1;q+1,...,2q)} | \xi_1)]^2 \\
&= \mathbb{E}[\mathbb{E}[h(Y_{\llbracket p \rrbracket}, \llbracket q \rrbracket) | \xi_1]^2] \\
&= \mathbb{E}[\psi_{\{(1),\emptyset\}}^{1,0} h^2]^2,
\end{align*}
$$

we find

$$
\mathbb{E}[(\hat{\mu}_N^{(1)})^2 (\hat{\mu}_N^{(2)})^2] = \mathbb{E}[\psi_{\{(1),\emptyset\}}^{1,0} h^2] + O(N^{-1}). \quad (35)
$$
• Next, we calculate:

$$ \mathbb{E}[(\hat{\mu}_N^{(1)})^2 \hat{\mu}_N^{(2)} \hat{\mu}_N^{(3)}] = T_N^{p,q}(I'_4, J'_4), $$

where $I'_4 = (\{1\}, \{1\}, \{2\}, \{3\})$ and $J'_4 = (\emptyset, \emptyset, \emptyset, \emptyset)$. Applying Lemma D.1 with

$$
\begin{align*}
\alpha(I'_4, J'_4) &= \mathbb{E}[h(Y_{(1,4,...,p+2;1,...,q)})h(Y_{(1,p+3,...,2p+1;q+1,...,2q)}) \\
&\quad h(Y_{(2,2p+2,...,3p;2q+1,...,3q)})h(Y_{(3,3p+1,...,4p-1;3q+1,...,4q)})] \\
&= \mathbb{E}[h(Y_{(1,...,p;1,...,q)})h(Y_{(1,p+1,...,2p-1;q+1,...,2q)})]^2 \mathbb{E}[h(Y_{(1,...,p;1,...,q)})]^2 \\
&= \mathbb{E}[\mathbb{E}[h(Y_{(1,...,p;1,...,q)})h(Y_{(1,p+1,...,2p-1;q+1,...,2q)}) | \xi_1]] \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)]^2 \\
&= \mathbb{E}[\mathbb{E}[h(Y_{\llbracket p \rrbracket},\llbracket q \rrbracket) | \xi_1]^2] \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h]^2 \\
&= \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)^2] \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h]^2,
\end{align*}
$$

we find

$$ \mathbb{E}[(\hat{\mu}_N^{(1)}(Y))^2 \hat{\mu}_N^{(2)}(Y) \hat{\mu}_N^{(3)}(Y)] = \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)^2] \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)]^2 + O(N^{-1}). \quad (36) $$

• Now, we calculate:

$$ \mathbb{E}[\hat{\mu}_N^{(1)}(Y) \hat{\mu}_N^{(2)}(Y) \hat{\mu}_N^{(3)}(Y) \hat{\mu}_N^{(4)}(Y)] = T_N^{p,q}(I''_4, J''_4), $$

where $I''_4 = (\{1\}, \{2\}, \{3\}, \{4\})$ and $J''_4 = (\emptyset, \emptyset, \emptyset, \emptyset)$. Applying Lemma D.1 with

$$
\begin{align*}
\alpha(I''_4, J''_4) &= \mathbb{E}[h(Y_{(1,5,...,p+3;1,...,q)})h(Y_{(2,p+4,...,2p+2;q+1,...,2q)}) \\
&\quad h(Y_{(3,2p+3,...,3p+1;2q+1,...,3q)})h(Y_{(4,3p+2,...,4p;3q+1,...,4q)})] \\
&= \mathbb{E}[h(Y_{\llbracket p \rrbracket},\llbracket q \rrbracket)]^4 \\
&= \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)]^4,
\end{align*}
$$

we find

$$ \mathbb{E}[\hat{\mu}_N^{(1)}(Y) \hat{\mu}_N^{(2)}(Y) \hat{\mu}_N^{(3)}(Y) \hat{\mu}_N^{(4)}(Y)] = \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)]^4 + O(N^{-1}). \quad (37) $$

Finally, injecting the calculated expressions (35), (36) and (37) in (34), we obtain

$$
\begin{align*}
\mathbb{E}[(\hat{v}_N^{h;1,0})^2] &= \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)^2] - 2\mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)^2] \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)]^2 \\
&\quad + \mathbb{E}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)]^4 + O(N^{-1}) \\
&= \mathbb{V}[\psi_{(\{1\},\emptyset)}^{1,0} h(Y)^2]^2 + O(N^{-1}) \\
&= (v_h^{1,0})^2 + O(N^{-1}) \\
&= \mathbb{E}[\hat{v}_N^{h;1,0}]^2 + O(N^{-1}),
\end{align*}
$$

where we have applied Proposition D.2 in the last step. This proves that $\mathbb{V}[\hat{v}_N^{h;1,0}] = O(N^{-1})$, concluding the proof. □
## Appendix E: Proofs of the results presented in Section 6

**Lemma E.1.** Let $Y$ be a matrix sampled from a Poisson W-graph model. Let $\bar{w}$, $f$ and $g$ be defined as in Section 6.3. For the kernel functions defined in Table 1, we have

$$ \bullet \mathbb{E}[h_{A,1}(Y_{(i_1,i_2:j_1,j_2)})) = \lambda^3 \iint \bar{w}(\xi, \eta)^2 d\xi d\eta, $$

$$ \bullet \mathbb{E}[h_{A,2}(Y_{(i_1,i_2:j_1,j_2)})) = \lambda^3 \iint \bar{w}(\xi, \eta) f(\xi) g(\eta) d\xi d\eta, $$

$$ \bullet \mathbb{E}[h_B(Y_{(i_1:j_1,j_2)})) = \lambda^2 \int f(\xi)^2 d\xi, $$

$$ \bullet \mathbb{E}[h_C(Y_{(i_1,i_2:j_1)})) = \lambda^2 \int g(\eta)^2 d\eta, $$

$$ \bullet \mathbb{E}[h_D(Y_{(i_1:j_1)})) = \lambda. $$

*Proof.* The result for $h_D$ is straightforward. For the other kernel functions:

$$ \begin{aligned} \mathbb{E}[h_{A,1}(Y_{(i_1,i_2:j_1,j_2)})) &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1}(Y_{i_1j_1} - 1)Y_{i_2j_2} | \xi, \eta]] \\ &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1}(Y_{i_1j_1} - 1) | \xi, \eta] \mathbb{E}[Y_{i_2j_2} | \xi, \eta]] \\ &= \mathbb{E}[\lambda^2 \bar{w}(\xi_{i_1}, \eta_{j_1})^2 \times \lambda] \\ &= \lambda^3 \iint \bar{w}(\xi, \eta)^2 d\xi d\eta. \end{aligned} $$

$$ \begin{aligned} \mathbb{E}[h_{A,2}(Y_{(i_1,i_2:j_1,j_2)})) &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1}Y_{i_1j_2}Y_{i_2j_2} | \xi, \eta]] \\ &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1} | \xi, \eta] \mathbb{E}[Y_{i_1j_2} | \xi, \eta] \mathbb{E}[Y_{i_2j_2} | \xi, \eta]] \\ &= \mathbb{E}[\lambda \bar{w}(\xi_{i_1}, \eta_{j_1}) \times \lambda \bar{w}(\xi_{i_1}, \eta_{j_2}) \times \lambda \bar{w}(\xi_{i_2}, \eta_{j_2})] \\ &= \lambda^3 \iint \left[ \bar{w}(\xi_{i_1}, \eta_{j_2}) \left( \int \bar{w}(\xi_{i_1}, \eta_{j_1}) d\eta_{j_1} \right) \left( \int \bar{w}(\xi_{i_2}, \eta_{j_2}) d\xi_{i_2} \right) \right] d\xi_{i_1} d\eta_{j_2} \\ &= \lambda^3 \iint \bar{w}(\xi, \eta) f(\xi) g(\eta) d\xi d\eta. \end{aligned} $$

$$ \begin{aligned} \mathbb{E}[h_B(Y_{(i_1:j_1,j_2)})) &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1}Y_{i_1j_2} | \xi, \eta]] \\ &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1} | \xi, \eta] \mathbb{E}[Y_{i_1j_2} | \xi, \eta]] \\ &= \mathbb{E}[\lambda \bar{w}(\xi_{i_1}, \eta_{j_1}) \times \lambda \bar{w}(\xi_{i_1}, \eta_{j_2})] \\ &= \lambda^2 \int \left[ \left( \int \bar{w}(\xi_{i_1}, \eta_{j_1}) d\eta_{j_1} \right) \left( \int \bar{w}(\xi_{i_1}, \eta_{j_2}) d\eta_{j_2} \right) \right] d\xi_{i_1} \\ &= \lambda^2 \int f(\xi)^2 d\xi. \end{aligned} $$

$$ \begin{aligned} \mathbb{E}[h_C(Y_{(i_1,i_2:j_1)})) &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1}Y_{i_2j_1} | \xi, \eta]] \\ &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1} | \xi, \eta] \mathbb{E}[Y_{i_2j_1} | \xi, \eta]] \\ &= \mathbb{E}[\lambda \bar{w}(\xi_{i_1}, \eta_{j_1}) \times \lambda \bar{w}(\xi_{i_2}, \eta_{j_1})] \\ &= \lambda^2 \int \left[ \left( \int \bar{w}(\xi_{i_1}, \eta_{j_1}) d\xi_{i_1} \right) \left( \int \bar{w}(\xi_{i_2}, \eta_{j_1}) d\xi_{i_2} \right) \right] d\eta_{j_1} \\ &= \lambda^2 \int g(\eta)^2 d\eta. \end{aligned} \quad \square $$

**Lemma E.2.** Let $Y$ be a matrix sampled from a W-graph model. Let $h_1$ and $h_2$ be the kernel functions defined as in Section 6.3. We have
$$ \bullet \mathbb{E}[h_1(Y_{(i_1,j_1,j_2)})] = \lambda^2 \int f(\xi)^2 d\xi = \lambda^2 F_2, $$

$$ \bullet \mathbb{E}[h_2(Y_{(i_1,i_2,j_1,j_2)})] = \lambda^2. $$

*Proof.* The proof for $h_1$ is identical to that for $h_B$ in the proof of Lemma E.1. For $h_2$,

$$ \begin{aligned} \mathbb{E}[h_2(Y_{(i_1,i_2,j_1,j_2)})] &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1}Y_{i_2j_2} | \xi, \eta]] \\ &= \mathbb{E}[\mathbb{E}[Y_{i_1j_1} | \xi, \eta]\mathbb{E}[Y_{i_2j_2} | \xi, \eta]] \\ &= \mathbb{E}[\lambda \times \lambda] \\ &= \lambda^2. \end{aligned} \quad \square $$

## Appendix F: Computation times for variance estimators

This section contains Tables 2 and 3 reporting the computation times for Algorithms A, B and C of 7.4 for estimating the asymptotic variance of $U_N^{h_1}$ and $U_N^{h_2}$.

As expected, Algorithm A is much slower than the proposed approach of this paper (Algorithm B and C). Notably, Algorithm C outspeeds Algorithm B by efficiently leveraging matrix operations, which are highly optimized in most computing libraries. While Algorithms A and B exhibit similar computational costs between estimating $V^{h_1}$ and $V^{h_2}$, Algorithm C is approximately twice as fast for $V^{h_1}$ compared to $V^{h_2}$. This reflects the fact that the matrix operations required for computing $U_N^{h_1}$ are simpler than for computing $U_N^{h_2}$.

<table><caption>Table 2
Average time in seconds for estimating $V^{h_1}$ with the three algorithms.</caption><thead><tr><th>N</th><th>Algorithm A</th><th>Algorithm B</th><th>Algorithm C</th></tr></thead><tbody><tr><td>8</td><td>1.92 × 10<sup>-2</sup></td><td>1.24 × 10<sup>-3</sup></td><td>3.52 × 10<sup>-4</sup></td></tr><tr><td>11</td><td>5.49 × 10<sup>-1</sup></td><td>2.75 × 10<sup>-3</sup></td><td>3.23 × 10<sup>-4</sup></td></tr><tr><td>16</td><td>7.82</td><td>1.08 × 10<sup>-2</sup></td><td>4.48 × 10<sup>-4</sup></td></tr><tr><td>22</td><td>1.56 × 10<sup>2</sup></td><td>4.54 × 10<sup>-2</sup></td><td>6.33 × 10<sup>-4</sup></td></tr><tr><td>32</td><td>-</td><td>1.90 × 10<sup>-1</sup></td><td>8.19 × 10<sup>-4</sup></td></tr><tr><td>45</td><td>-</td><td>8.19 × 10<sup>-1</sup></td><td>1.27 × 10<sup>-3</sup></td></tr><tr><td>64</td><td>-</td><td>3.45</td><td>2.46 × 10<sup>-3</sup></td></tr><tr><td>90</td><td>-</td><td>1.24 × 10<sup>1</sup></td><td>5.63 × 10<sup>-3</sup></td></tr><tr><td>128</td><td>-</td><td>2.52 × 10<sup>1</sup></td><td>1.66 × 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>Table 3
Average time in seconds for estimating $V^{h_2}$ with the three algorithms.</caption><thead><tr><th>N</th><th>Algorithm A</th><th>Algorithm B</th><th>Algorithm C</th></tr></thead><tbody><tr><td>8</td><td>3.40 × 10<sup>-2</sup></td><td>1.09 × 10<sup>-3</sup></td><td>4.63 × 10<sup>-4</sup></td></tr><tr><td>11</td><td>6.48 × 10<sup>-1</sup></td><td>4.26 × 10<sup>-3</sup></td><td>5.26 × 10<sup>-4</sup></td></tr><tr><td>16</td><td>8.71</td><td>1.15 × 10<sup>-2</sup></td><td>7.41 × 10<sup>-4</sup></td></tr><tr><td>22</td><td>1.25 × 10<sup>2</sup></td><td>4.57 × 10<sup>-2</sup></td><td>1.01 × 10<sup>-3</sup></td></tr><tr><td>32</td><td>-</td><td>1.94 × 10<sup>-1</sup></td><td>1.45 × 10<sup>-3</sup></td></tr><tr><td>45</td><td>-</td><td>8.99 × 10<sup>-1</sup></td><td>2.53 × 10<sup>-3</sup></td></tr><tr><td>64</td><td>-</td><td>3.47</td><td>4.38 × 10<sup>-3</sup></td></tr><tr><td>90</td><td>-</td><td>1.47 × 10<sup>1</sup></td><td>1.05 × 10<sup>-2</sup></td></tr><tr><td>128</td><td>-</td><td>2.81 × 10<sup>1</sup></td><td>3.28 × 10<sup>-2</sup></td></tr></tbody></table>
## Acknowledgments

The authors would like to thank the anonymous referees, the Associate Editor and the Editor for their constructive comments that improved the quality of this paper.

## Funding

This work was funded by the grant ANR-18-CE02-0010-01 of the French National Research Agency ANR (project EcoNet) and a grant from Région Île-de-France.

## References

[1] ALI, W., RITO, T., REINERT, G., SUN, F. and DEANE, C. M. (2014). Alignment-free protein interaction network comparison. *Bioinformatics* **30** i430–i437.

[2] ARVESEN, J. N. (1969). Jackknifing *U*-statistics. *The Annals of Mathematical Statistics* **40** 2076–2100. MR0264805

[3] ASTA, D. and SHALIZI, C. R. (2015). Geometric network comparison. *Uncertainty in Artificial Intelligence* 102–110.

[4] BAKER, N. J., KAARTINEN, R., ROSLIN, T. and STOUFFER, D. B. (2015). Species’ roles in food webs show fidelity across a highly variable oak forest. *Ecography* **38** 130–139.

[5] BARRAT, A., BARTHELEMY, M., PASTOR-SATORRAS, R. and VESPIGNANI, A. (2004). The architecture of complex weighted networks. *Proceedings of the National Academy of Sciences* **101** 3747–3752.

[6] BEARMAN, P. S., MOODY, J. and STOVEL, K. (2004). Chains of affection: The structure of adolescent romantic and sexual networks. *American Journal of Sociology* **110** 44–91.

[7] BHATTACHARYA, B. B., CHATTERJEE, A. and JANSON, S. (2023). Fluctuations of subgraph counts in graphon based random graphs. *Combinatorics, Probability and Computing* **32** 428–464. MR4574856

[8] BHATTACHARYA, B. B., DAS, S., MUKHERJEE, S. and MUKHERJEE, S. (2024). Fluctuations of quadratic chaos. *Communications in Mathematical Physics* **405** 237. MR4797748

[9] BHATTACHARYYA, S. and BICKEL, P. J. (2015). Subsampling bootstrap of count features of networks. *The Annals of Statistics* **43** 2384–2411. MR3405598

[10] BICKEL, P. J., CHEN, A. and LEVINA, E. (2011). The method of moments and degree distributions for network models. *The Annals of Statistics* **39** 2280–2301. MR2906868

[11] BICKEL, P. J. and FREEDMAN, D. A. (1981). Some asymptotic theory for the bootstrap. *The Annals of Statistics* **9** 1196–1217. MR0630103

[12] BILLINGSLEY, P. (1995). *Probability and Measure*, 3rd ed. John Wiley & Sons. MR1324786

[13] BLOM, G. (1976). Some properties of incomplete *U*-statistics. *Biometrika* **63** 573–580. MR0474582

[14] BROWN, B. M. and KILDEA, D. G. (1978). Reduced *U*-statistics and the Hodges-Lehmann estimator. *The Annals of Statistics* **6** 828–835. MR0491556
[15] CALLAERT, H. and VERAVERBEKE, N. (1981). The order of the normal approximation for a studentized U-statistic. *The Annals of Statistics* **9** 194–200. MR0600547

[16] CHUNG, F. and LU, L. (2002). The average distances in random graphs with given expected degrees. *Proceedings of the National Academy of Sciences* **99** 15879–15882. MR1944974

[17] COULSON, M., GAUNT, R. E. and REINERT, G. (2016). Poisson approximation of subgraph counts in stochastic block models and a graphon model. *ESAIM: Probability and Statistics* **20** 131–142. MR3528620

[18] DAVEZIES, L., D'HAULTFŒUILLE, X. and GUYONVARCH, Y. (2021). Empirical process results for exchangeable arrays. *The Annals of Statistics* **49** 845–862. MR4255110

[19] DIACONIS, P. and JANSON, S. (2008). Graph limits and exchangeable random graphs. *Rendiconti di Matematica e delle sue Applicazioni. Serie VII* **28** 33–61. MR2463439

[20] DOOB, J. L. (1953). *Stochastic Processes*. Wiley New York. MR0058896

[21] DORMANN, C. F., FRÜND, J., BLÜTHGEN, N. and GRUBER, B. (2009). Indices, graphs and null models: Analyzing bipartite ecological networks. *The Open Ecology Journal* **2**.

[22] EFRON, B. (1979). Bootstrap methods: Another look at the jackknife. *The Annals of Statistics* **7** 1–26. MR0515681

[23] EMMERT-STREIB, F., DEHMER, M. and SHI, Y. (2016). Fifty years of graph matching, network alignment and network comparison. *Information Sciences* **346** 180–197. MR3475436

[24] FRANK, O. and STRAUSS, D. (1986). Markov graphs. *Journal of the American Statistical Association* **81** 832–842. MR0860518

[25] GAO, C. and LAFFERTY, J. (2017). Testing network structure using relations between small subgraph probabilities. *arXiv preprint arXiv:1704.06742*.

[26] GOVAERT, G. and NADIF, M. (2003). Clustering with block mixture models. *Pattern Recognition* **36** 463–473.

[27] HOEFFDING, W. (1948). A class of statistics with asymptotically normal distribution. *The Annals of Mathematical Statistics* **19** 293–325. MR0026294

[28] HOEFFDING, W. (1961). The strong law of large numbers for U-statistics. Technical Report, North Carolina State University. Dept. of Statistics.

[29] HOFF, P. D., RAFTERY, A. E. and HANDCOCK, M. S. (2002). Latent space approaches to social network analysis. *Journal of the American Statistical Association* **97** 1090–1098. MR1951262

[30] JANSON, S. (1984). The asymptotic distributions of incomplete U-statistics. *Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete* **66** 495–505. MR0753810

[31] JANSON, S. and NOWICKI, K. (1991). The asymptotic distributions of generalized U-statistics with applications to random graphs. *Probability Theory and Related Fields* **90** 341–375. MR1133371

[32] KALLENBERG, O. (2005). *Probabilistic Symmetries and Invariance Principles* **9**. Springer. MR2161313

[33] KASHTAN, N., ITZKOVITZ, S., MILO, R. and ALON, U. (2004). Topological generalizations of network motifs. *Physical Review E* **70** 031909. MR2133833

[34] KAUR, G. and RÖLLIN, A. (2021). Higher-order fluctuations in dense random graph models. *Electronic Journal of Probability* **26** 1–36. MR4346670
[35] LANUZA, J. B., ALLEN-PERKINS, A. and BARTOMEUS, I. (2023). The non-random assembly of network motifs in plant-pollinator networks. *Journal of Animal Ecology* **92** 760–773.

[36] LAURITZEN, S., RINALDO, A. and SADEGHI, K. (2018). Random networks, graphical models and exchangeability. *Journal of the Royal Statistical Society Series B: Statistical Methodology* **80** 481–508. MR3798875

[37] LE MINH, T. (2023). *U-statistics on bipartite exchangeable networks*. *ESAIM: Probability and Statistics* **27** 576–620. MR4602388

[38] LE MINH, T. (2024). Characterization of the asymptotic behavior of *U*-statistics on row-column exchangeable matrices. *arXiv preprint arXiv:2401.07876*.

[39] MAUGIS, P. A. G., OLHEDE, S. C., PRIEBE, C. E. and WOLFE, P. J. (2020). Testing for equivalence of network distribution using subgraph counts. *Journal of Computational and Graphical Statistics* **29** 455–465. MR4153173

[40] MENZEL, K. (2021). Bootstrap with cluster-dependence in two or more dimensions. *Econometrica* **89** 2143–2188. MR4325248

[41] MICHALSKA-SMITH, M. J. and ALLESINA, S. (2019). Telling ecological networks apart by their structure: A computational challenge. *PLoS Computational Biology* **15** e1007076. Data available at github.com/mjsmith037/Classifying_Bipartite_Networks.

[42] NAULET, Z., ROY, D. M., SHARMA, E. and VEITCH, V. (2021). Bootstrap estimators for the tail-index and for the count statistics of graphex processes. *Electronic Journal of Statistics* **15** 282–325. MR4195774

[43] NEWMAN, M. E. J. (2001). The structure of scientific collaboration networks. *Proceedings of the National Academy of Sciences* **98** 404–409. MR1812610

[44] NOWICKI, K. and SNIJDERS, T. A. B. (2001). Estimation and prediction for stochastic blockstructures. *Journal of the American Statistical Association* **96** 1077–1087. MR1947255

[45] O'NEIL, K. A. and REDNER, R. A. (1993). Asymptotic distributions of weighted *U*-statistics of degree 2. *The Annals of Probability* **21** 1159–1169. MR1217584

[46] ORSINI, C., DANKULOV, M. M., COLOMER-DE SIMÓN, P., JAMAKOVIC, A., MAHADEVAN, P., VAHDAT, A., BASSLER, K. E., TOROCZKAI, Z., BOGUNÁ, M., CALDARELLI, G., FORTUNATO, S. and KRIOUKOV, D. (2015). Quantifying randomness in real networks. *Nature Communications* **6** 1–10.

[47] OUADAH, S., LATOUCHE, P. and ROBIN, S. (2022). Motif-based tests for bipartite networks. *Electronic Journal of Statistics* **16** 293–330. MR4360230

[48] PRŽULJ, N., CORNEIL, D. G. and JURISICA, I. (2004). Modeling interactome: Scale-free or geometric? *Bioinformatics* **20** 3508–3515.

[49] PRZYTYCKA, T. M. (2006). An important connection between network motifs and parsimony models. In *Annual International Conference on Research in Computational Molecular Biology* 321–335. Springer.

[50] REINERT, G. and RÖLLIN, A. (2010). Random subgraph counts and *U*-statistics: Multivariate normal approximation via exchangeable pairs and embedding. *Journal of Applied Probability* **47** 378–393. MR2668495

[51] SCHUCANY, W. R. and BANKSON, D. M. (1989). Small sample variance estimators for *U*-statistics. *Australian Journal of Statistics* **31** 417–426.

[52] SEN, P. K. (1960). On some convergence properties of *U*-statistics. *Calcutta Statistical Association Bulletin* **10** 1–18. MR0119286
[53] SEN, P. K. (1977). Some invariance principles relating to jackknifing and their role in sequential analysis. *The Annals of Statistics* **5** 316–329. MR0443175

[54] SHAO, M., XIA, D., ZHANG, Y., WU, Q. and CHEN, S. (2022). Higher-order accurate two-sample network inference and network hashing. *arXiv preprint arXiv:2208.07573*.

[55] SHAPIRO, C. P. and HUBERT, L. (1979). Asymptotic normality of permutation statistics derived from weighted sums of bivariate functions. *The Annals of Statistics* **7** 788–794. MR0532242

[56] SHEN-ORR, S. S., MILO, R., MANGAN, S. and ALON, U. (2002). Network motifs in the transcriptional regulation network of Escherichia coli. *Nature Genetics* **31** 64–68.

[57] SILVERMAN, B. W. (1976). Limit theorems for dissociated random variables. *Advances in Applied Probability* **8** 806–819. MR0431348

[58] SIMMONS, B. I., CIRTWILL, A. R., BAKER, N. J., WAUCHOPE, H. S., DICKS, L. V., STOUFFER, D. B. and SUTHERLAND, W. J. (2019). Motifs in bipartite ecological networks: Uncovering indirect interactions. *Oikos* **128** 154–170.

[59] SNIJDERS, T. A. B. (1981). The degree variance: An index of graph heterogeneity. *Social Networks* **3** 163–174. MR0639871

[60] SNIJDERS, T. A. B. and NOWICKI, K. (1997). Estimation and prediction for stochastic blockmodels for graphs with latent block structure. *Journal of Classification* **14** 75–100. MR1449742

[61] STOUFFER, D. B., CAMACHO, J., JIANG, W. and NUNES AMARAL, L. A. (2007). Evidence for the existence of a robust pattern of prey selection in food webs. *Proceedings of the Royal Society B: Biological Sciences* **274** 1931–1940. MR2306896

[62] TANTARDINI, M., IEVA, F., TAJOLI, L. and PICCARDI, C. (2019). Comparing methods for comparing networks. *Scientific Reports* **9** 1–19.

[63] VAN DER VAART, A. W. (2000). *Asymptotic Statistics*. Cambridge University Press. MR1652247

[64] WANG, Y., PELEKIS, C. and RAMON, J. (2015). U-statistics on network-structured data with kernels of degree larger than one. In *Workshop on Statistically Sound Data Mining* 37–48.

[65] YOUNG, S. J. and SCHEINERMAN, E. R. (2007). Random dot product graph models for social networks. In *Workshop on Algorithms and Models for the Web-Graph* 138–149. Springer. MR2504912

[66] ZHAO, L., BEVERLIN, B., NETOFF, T. and NYKAMP, D. Q. (2011). Synchronization from second order network connectivity statistics. *Frontiers in Computational Neuroscience* **5** 28.

[67] ZHOU, T., REN, J., MEDO, M. and ZHANG, Y. C. (2007). Bipartite network projection and personal recommendation. *Physical Review E* **76** 046115.