# Heat diffusion distance processes : a statistically founded method to
analyze graph data sets

Etienne Lasalle,
Laboratoire de Mathématiques d’Orsay,
Université Paris-Saclay, Orsay, France
<EMAIL_ADDRESS>
###### Abstract

We propose two multiscale comparisons of graphs using heat diffusion, allowing
to compare graphs without node correspondence or even with different sizes.
These multiscale comparisons lead to the definition of Lipschitz-continuous
empirical processes indexed by a real parameter. The statistical properties of
empirical means of such processes are studied in the general case. Under mild
assumptions, we prove a functional Central Limit Theorem, as well as a
Gaussian approximation with a rate depending only on the sample size. Once
applied to our processes, these results allow to analyze data sets of pairs of
graphs. We design consistent confidence bands around empirical means and
consistent two-sample tests, using bootstrap methods. Their performances are
evaluated by simulations on synthetic data sets.

##  1 Introduction

Considering the current growth of available data and the modeling power of
networks, methods to analyze graph-structured data have gained interest over
the last few decades. Particular attention has been devoted to designing
notions of distance between graphs. The design of these notions is highly
constrained by the working framework. In particular, different types of
information can be used depending on whether the graphs are directed or
undirected, weighted or unweighted, have the same size or not. Another key
factor to define distances is whether a node correspondence (NC) is known or
not. In the case of a known NC, we can consider the graphs to be defined on
the same vertex set and comparisons can be made at the edge scale. In this
context, people have applied various metrics to compare adjacency matrices,
Laplacian matrices, heat kernels [15] and other matrices whose entries
represent quantities associated to pairs of nodes [19]. On the other hand,
when no NC is available, graphs are often compared at a mesoscopic or
macroscopic level using structural summaries. People have used global
statistics on graphs like degree distributions, network diameters, or
clustering coefficients [26]. Another well studied approach is to consider
graphlets [26], i.e. small given subgraphs that are counted in graphs. Then,
various methods have been developed to compare the graphlet counts [25, 33, 1,
12]. Some work has also been pursued to exploit the structural information
carried by spectra of operators [32, 14].

Another powerful way to encode structural information about graphs is to use
diffusion processes, like heat diffusion. When working with weighted graphs,
one can interpret weights as the thermal conductivity of edges, meaning that
heat diffuses faster along edges with higher weights. Note that unweighted
graphs can always be seen as weighted graphs with weights in $\\{0,1\\}$.
Given initial conditions, the way heat diffuses can be used to characterize
and compare graphs [8, 7, 15, 30]. This approach is appealing as it allows to
analyze graphs at different scales by looking at different diffusion times
$t$. For small values of $t$, the diffusion only concerns a small neighborhood
of the initially heated nodes, while for larger values it involves larger and
possibly more complex structures, taking into account topological properties
of the graph. Thus, the choice of relevant and informative diffusion times is
essential.

For more references on comparisons of graphs, we send the reader to [27, 10,
29] and references therein.

###  1.1 Our contributions

While a lot of the above notions of distances are often supported by
experimental results and applications to learning or data mining tasks, they
usually lack statistical foundations. In this context, we provide new tools to
analyze and compare graphs or even data sets of graphs, that benefit from
statistical guarantees. Our methods take advantage of the desirable multiscale
property of heat diffusion. Moreover, one of our methods can deal with graphs
without known NC or even graphs of different sizes, by using topological
descriptors from Topological Data Analysis (TDA). To circumvent the difficulty
of choosing a suitable diffusion time, we opt to take into account the whole
diffusion process. As a result, we define two real-valued processes, indexed
by all the diffusion times in $[0,T]$ for some $T>0$, representing comparisons
of heat distributions.

The first process, called Heat Kernel Distance (HKD) process is defined by
comparing heat kernels with the Frobenius norm. In this case, an NC between
graphs needs to be known for the entry-wise comparison of the heat kernels to
be meaningful.

The second process, called Heat Persistence Distance (HPD), is defined using
tools from TDA and can deal with graphs of different sizes. To do so, each
graph is equipped with a real-valued function, the Heat Kernel Signature (HKS)
[28, 17], defined on the vertex set. Then, graphs are converted into
topological descriptors called persistence diagrams. They are multisets of
points in $\mathbb{R}^{2}$, encoding how topological features, like connected
components and loops, evolve along with the families of sublevel and
superlevel subgraphs. The diagrams are then compared with the so-called
Bottleneck distance. Using persistence diagrams allows switching from node-
based representations of graphs to comparable topological summaries, hence
requiring no assumption on graph sizes and NC.

To statistically study the HKD and HPD processes, we prove general results on
Lipschitz-continuous real-valued empirical processes indexed by a real
parameter. Namely, we show that they verify a functional Central Limit Theorem
and admit Gaussian approximations with rates depending only on the sample
sizes. These results ensure the asymptotic validity of bootstrap methods to
design confidence bands around empirical mean processes, as well as consistent
two-sample tests. They are applied to the HKD and HPD processes and could be
applied to any other Lipschitz-continuous processes indexed by a real
parameter under mild assumptions. We illustrate these results on simulated
data sets of pairs of graphs, drawn from various models: Erdős-Rényi model,
stochastic block model, and random geometric graph models.

###  1.2 Organisation

The rest of the paper is organized as follows. Section 2 introduces the graph
framework and heat diffusion on graphs. We define the HKD and HPD processes,
and present their statistical properties. These results being actually not
restricted to the HKD and HPD processes, we generalize the study of such
processes. In Section 3, we introduce a framework for general continuous real-
valued empirical processes indexed by a real parameter. We prove their
statistical properties and present some consequences on bootstrap methods.
Finally, we illustrate the construction of confidence bands and two-sample
tests in Section 4 using several generative models of random graphs. All
Python codes are freely available at
https://github.com/elasalle/HeatDistanceProcess.

##  2 Study of graph data using heat distance processes

###  2.1 Background and definitions

####  2.1.1 Graphs

For $0\leq w_{\min}\leq w_{\max}$, we denote by
$\mathcal{G}_{n}(w_{\min},w_{\max})$ the set of undirected weighted graphs of
size $n$, without self-loop and whose weights are in
$\\{0\\}\cup[w_{\min},w_{\max}]$. The special case of unweighted graphs
correspond to $w_{\min}=w_{\max}=1$. For clarity in the notation, we remove
the $w_{\min}$ and $w_{\max}$, whenever there is no ambiguity. We also
consider $\mathcal{G}^{n}$ (with $n$ as an exponent) the set of graphs of size
at most $n$, i.e. $\mathcal{G}^{n}=\cup_{1\leq i\leq n}\mathcal{G}_{i}$. For a
graph $G$ in $\mathcal{G}_{n}$, we denote by $W(G)$ its weight matrix (or
adjacency matrix), i.e. the $n\times n$ symmetric matrix whose
$(i,j)$-coefficient is the weight $w_{i,j}$ of edge $\\{i,j\\}$. $D(G)$
denotes the diagonal matrix whose entry $D(G)_{i,i}$ is the degree of node $i$
defined by $\sum_{j}w_{i,j}(G)$. The combinatorial Laplacian $L(G)$ is defined
by $D(G)-W(G)$. Taking non-negative weights ensures that $L(G)$ is a real
symmetric positive-semidefinite matrix. From now on, we forget the dependence
in $G$ in the notation, whenever there is no ambiguity. Let
$\lambda_{1}\leq\dots\leq\lambda_{n}$ be the eigenvalues of $L$ and let
$(\phi_{1},\dots,\phi_{n})$ be a family of orthonormal eigenvectors. We denote
by $\Lambda$ the diagonal matrix containing the eigenvalues on the diagonal
and $\phi$ the matrix whose columns are the $\phi_{i}$’s so that $L$ admits
the following decomposition

$L=\phi\Lambda\phi^{T}=\sum\limits_{k=1}^{n}\lambda_{k}\phi_{k}\phi_{k}^{T}.$
(2.1)

Note that $\lambda_{1}=0$ and that $\phi_{1}$ can always be chosen to be the
vector whose entries are equal to $1/\sqrt{n}$. In the following, this choice
will always be made.

####  2.1.2 Persistence on graphs

We present here the basics of ordinary and extended persistence. We send the
reader to [6, 9, 23] for a complete description of these theories.

Persistence theory allows to study the topology of topological spaces in a
multiscale manner. Usually, given a topological space $X$ and a continuous
real-valued function $f:X\to\mathbb{R}$, one considers the family of sublevel
sets $X_{\alpha}:=\\{x\in X,\ f(x)\leq\alpha\\}$, for $\alpha$ varying from
$-\infty$ to $+\infty$. Ordinary persistence records the levels at which
topological features (connected components, loops, cavities, or higher
dimensional holes…) appear and disappear. For each feature, its birth and
death levels are stored as the coordinates $(b,d)$ of a point in
$\mathbb{R}^{2}$. The multiset of these points is called a persistence
diagram. This framework can be applied to graphs.

Let $G=(V,E)$ be a graph with vertex set $V$ and edge set $E$, and $f$ be a
real-valued function on $V$. Consider the family of sublevel subgraphs
$(G_{\alpha})_{\alpha\in\mathbb{R}}$, where
$G_{\alpha}=(V_{\alpha},E_{\alpha})$ with $V_{\alpha}=\\{v\in V,\
f(v)\leq\alpha\\}$ and $E_{\alpha}=\\{\\{v,v^{\prime}\\}\in E,\
v,v^{\prime}\in V_{\alpha}\\}$. Across the family of sublevel graphs, as
$\alpha$ increases, we can record birth and death levels of connected
components, and birth levels of loops. As a connected component dies when it
gets connected to an older connected component, remark that the connected
components of $G$ will never die. Similarly, as a loop dies when it gets
"filled-in" by a 2-dimensional object, loops appearing in $G_{\alpha}$ (an
object of maximal dimension 1) will never die. To prevent topological features
from having no death levels (or infinite death levels), the theory of extended
persistence suggests also considering the family of superlevel subgraphs.
Define $G^{\alpha}=(V^{\alpha},E^{\alpha})$ similarly to $G_{\alpha}$ with
$V^{\alpha}=\\{v\in V,\ f(v)\geq\alpha\\}$. A death level is now assigned to
connected components of $G$ and loops as the level at which they appear in the
family of superlevel subgraphs when $\alpha$ decreases. Additionally, we
record the birth and death of connected components in the family of superlevel
subgraphs. Hence, extended persistence is able to detect four types of
topological features and extract their birth and death levels $(b,d)$
corresponding to four types of points :

  * •

$Ord_{0}$ : birth and death of a connected component in $(G_{\alpha})$.

  * •

$Rel_{1}$ : birth and death of a connected component in $(G^{\alpha})$.

  * •

$Ext_{0}$ : birth and death of a connected component of $G$ when using both
$(G_{\alpha})$ and $(G^{\alpha})$.

  * •

$Ext_{1}$ : birth and death of a loop when using both $(G_{\alpha})$ and
$(G^{\alpha})$.

These four types of topological features can be seen as downward branches,
upward branches, connected components, and loops, respectively; with the
orientation being taken with respect to $f$. For a graph $G$ and a function
$f$ on its vertices, we will denote by $Ord_{0}(G,f)$, $Rel_{1}(G,f)$,
$Ext_{0}(G,f)$ and $Ext_{1}(G,f)$ the persistence diagrams containing the
corresponding points. In the following, $Dg(G,f)$ will generically denote any
of these four diagrams. We send the reader to [4, Section 2.1] for a more
precise and illustrative presentation of extended persistence diagrams on
graphs.

The space of diagrams can be equipped with the Bottleneck distance $d_{B}$. We
recall its definition. Let $\mu$ and $\nu$ be two diagrams, i.e. two multisets
of points in $\mathbb{R}^{2}$, and let $\Delta:=\\{(a,a),a\in\mathbb{R}\\}$ be
the diagonal. Denote by $\Pi(\mu,\nu)$ the set of bijections from
$\mu\cup\Delta$ to $\nu\cup\Delta$. $d_{B}$ is defined by

$d_{B}(\mu,\nu):=\underset{\pi\in\Pi(\mu,\nu)}{\inf}\
\underset{x\in\mu\cup\Delta}{\sup}\ \|x-\pi(x)\|_{\infty}.$ (2.2)

We state a stability result for extended persistence diagrams computed on
graphs. It is a consequence of a more general stability result for persistence
diagrams.

######  Theorem 2.1 ([5, 6]).

For all graph $G=(V,E)$, for all $f,f^{\prime}:V\to\mathbb{R}$, and for all
diagram construction $Dg$ among $Ord_{0}$, $Rel_{1}$, $Ext_{0}$ and $Ext_{1}$,

$d_{B}(Dg(G,f),Dg(G,f^{\prime}))\leq\|f-f^{\prime}\|_{\infty},$ (2.3)

with $\|f\|_{\infty}=\max\\{\left|f(v)\right|,\ v\in V\\}$.

###  2.2 Heat distance processes

Let $G$ be a graph in $\mathcal{G}_{n}$ and let $L$ be its Laplacian. For
$t\geq 0$, let $u_{t}\in\mathbb{R}^{n}$ be the vector whose $i$-th coefficient
represents the amount of heat of node $i$ at time $t$. Then, $u_{t}$ follows
the heat equation :

$\forall t\geq 0,\quad\frac{d}{dt}u_{t}=-Lu_{t}$ (2.4)

The solution is given by $u_{t}=e^{-tL}u_{0}$. The matrix $e^{-tL}$ is called
the heat kernel and describes how heat diffuses in the graph. The $i$-th
column of $e^{-tL}$ contains the amount of heat of each node at time $t$, when
a single unit of heat was placed at node $i$ at time $t=0$. From (2.1), the
heat kernel decomposes in

$e^{-tL}=\phi
e^{-t\Lambda}\phi^{T}=\sum\limits_{k=1}^{n}e^{-t\lambda_{k}}\phi_{k}\phi_{k}^{T}.$
(2.5)

####  2.2.1 Heat Kernel Distance

Assume that we know the NC for graphs in $\mathcal{G}_{n}$ and that we number
the nodes such that the identity mapping gives the correspondences. Hence,
comparing adjacency matrices, Laplacians, or heat kernels entry-wise becomes
meaningful. Here we compare graphs through their heat kernels. For two graphs
$G$ and $G^{\prime}$, define their Heat Kernel Distance (HKD) at time $t$ by

$D_{t}((G,G^{\prime}))=\|e^{-tL}-e^{-tL^{\prime}}\|_{F},$ (2.6)

where $L$ and $L^{\prime}$ are the laplacians of $G$ and $G^{\prime}$
respectively, and $\|\cdot\|_{F}$ denotes the Frobenius norm. This notion of
distance was introduced by [15]. To turn the HKD into a parameter-free notion
of distance, [15] define the Graph Diffusion Distance as
$\max_{t}D_{t}((G,G^{\prime}))$. This has the drawback of comparing different
pairs of graphs at different times. Instead, our approach consists in using
the whole function $t\to D_{t}((G,G^{\prime}))$. More precisely, considering a
probability distribution $P$ on $\mathcal{G}_{n}\times\mathcal{G}_{n}$ and a
random pair of graphs $(G,G^{\prime})\sim P$, we are interested in the
stochastic process $\\{D_{t}((G,G^{\prime})),\ t\in[0,T]\\}$ for some $T>0$.
That is, the process obtained by evaluating the functions of the family
$\mathcal{F}_{HKD}:=\\{D_{t},\ t\in[0,T]\\}$ on a random pair of graphs
$(G,G^{\prime})$. This framework corresponds to the general framework of
_empirical processes_. The properties of the process associated with
$\mathcal{F}_{HKD}$ are studied in Section 2.3, while a more general study of
such empirical processes is carried out in Section 3.

####  2.2.2 Heat Persistence Distance

In practice, the NC between graphs is not always known. Additionally one may
be interested in comparing graphs of different sizes. In these cases, HKD can
not be computed. To circumvent these issues and following ideas from [4], we
define the Heat Persistence Distance (HPD) by using extended persistence
diagrams computed with the Heat Kernel Signature (HKS). These persistence
diagrams can be compared with the Bottleneck distance $d_{B}$ without any
assumption on graph sizes and node identification.

The HKS was first introduced by [28] for the study of shapes. Here we restrict
ourselves to the definition of the HKS on graphs of [17]. For a graph $G$ of
size $n$ with vertex set $V=\\{1,\dots,n\\}$, the HKS at time $t$ is the
function $h_{t}(G):V\to\mathbb{R}$ such that

$h_{t}(G)(i)=\sum\limits_{k=1}^{n}e^{-t\lambda_{k}}\phi_{k}(i)^{2},\quad 1\leq
i\leq s.$ (2.7)

Intuitively, the image of $h_{t}(G)$ corresponds to the diagonal of the heat
kernel $e^{-tL}$. Hence, $h_{t}(G)(i)$ represents the remaining amount of heat
at node $i$ after a diffusion time $t$, when a single unit of heat was placed
at node $i$ at time $t=0$. For each value of $t$, the HKS provides a function
on the vertices of a graph, that we use to compute an extended persistence
diagram. The HPD at time $t$ between two graphs $G,G^{\prime}$ in
$\mathcal{G}^{n}$ is defined by

$H_{t}((G,G^{\prime}))=\underset{Dg}{\max}\
d_{B}(Dg(G,h_{t}(G)),Dg(G^{\prime},h_{t}(G^{\prime}))),$ (2.8)

where the maximum is taken over the four diagram constructions $Ord_{0}$,
$Rel_{1}$, $Ext_{0}$ and $Ext_{1}$. In our simulations, persistence diagrams
are computed by following the approach of [4] and using the Gudhi library
[21].

Similarly to $\mathcal{F}_{HKD}$, we define the family
$\mathcal{F}_{HPD}:=\\{H_{t},\ t\in[0,T]\\}$ in order to study the induced
stochastic process : $\\{H_{t}((G,G^{\prime})),\ t\in[0,T]\\}$, for some
random pair of graphs $(G,G^{\prime})\in\mathcal{G}^{n}\times\mathcal{G}^{n}$.

###  2.3 Results

The goal of this section is to show that the HKD and HPD processes admit a
functional version of the Central Limit Theorem, as well as a Gaussian
approximation. As we will show in Section 3, this is a special case of a more
general result on uniformly bounded Lipschitz-continuous processes. For the
sake of completeness of this section, we choose to state here the results on
the distance processes, before exposing and proving in Section 3 the general
results.

Let us start by proving a lemma on the extremal laplacian eigenvalues of
graphs in $\mathcal{G}_{n}$. We define $\Lambda_{\min}$ and $\Lambda_{\max}$
as, respectively, the minimal and maximal positive laplacian eigenvalues of
graphs in $\mathcal{G}_{n}$ :

$\displaystyle\Lambda_{\min}$ $\displaystyle:=\inf\\{\lambda>0,\
\textnormal{s.t. }\lambda\textnormal{ is an eigenvalue of
}L(G),G\in\mathcal{G}_{n}\\}$ (2.9) $\displaystyle\Lambda_{\max}$
$\displaystyle:=\sup\\{\lambda>0,\ \textnormal{s.t. }\lambda\textnormal{ is an
eigenvalue of }L(G),G\in\mathcal{G}_{n}\\}.$ (2.10)

######  Lemma 2.1.

$\Lambda_{\min}$ and $\Lambda_{\max}$ satisfy the following bounds :

$\displaystyle\Lambda_{\min}$ $\displaystyle\geq\frac{8w_{\min}}{n^{2}}$
(2.11) $\displaystyle\Lambda_{\max}$ $\displaystyle\leq nw_{\max}.$ (2.12)

Note that (2.11) will not be used in the rest of the paper. Still, we choose
to present it as we believe it could be of future use.

###### Proof of Lemma 2.1.

Take $G\in\mathcal{G}_{n}$, we want to prove that $\lambda_{min}$, the
smallest positive eigenvalue of $L(G)$, verifies $\lambda_{min}\geq
8w_{\min}n^{-2}$. To do so, we apply a Cheeger-type inequality. First assume
that $G$ is connected, hence $\lambda_{min}=\lambda_{2}(G)$. Let $V$ be the
set of vertices of $G$. Following [13, Section 3] we define the average
minimal cut of $G$ by

$\gamma(G)=\min_{\emptyset\neq U\subsetneq V}\sum\limits_{i\in U,\ j\in
V\backslash U}\frac{w_{i,j}}{\left|U\right|(n-\left|U\right|)}.$ (2.13)

As $G$ is connected, there exists at least one edge with a weight greater than
$w_{\min}$ joining $U$ and $V\backslash U$. Hence,

$\sum\limits_{i\in U,\ j\in V\backslash U}w_{i,j}\geq w_{\min}.$ (2.14)

Moreover, for all $U$, $\left|U\right|(n-\left|U\right|)\leq n^{2}/4$. This
yields to $\gamma(G)\geq 4w_{\min}/n^{2}$. From [13, Theorem 2], we have that
$\lambda_{2}(G)\geq 2\gamma(G)$, hence $\lambda_{2}(G)\geq 8w_{\min}n^{-2}$.
If now $G$ is not connected, one can check that there exists $G_{sub}$ a
connected subgraph of $G$ of size $n_{sub}$, such that
$\lambda_{min}=\lambda_{2}(G_{sub})$. This gives

$\lambda_{min}=\lambda_{2}(G_{sub})\geq\frac{8w_{\min}}{n_{sub}^{2}}\geq\frac{8w_{\min}}{n^{2}}.$
(2.15)

This finishes the proof of the first bound.

We now prove the second bound. The largest eigenvalue of L(G) is denoted by
$\lambda_{n}$. Letting $x$ be an eigenvector associated to $\lambda_{n}$
verifying $\|x\|_{2}=1$, we have

$\lambda_{n}=x^{T}L(G)x=\sum\limits_{i<j}w_{i,j}(x_{i}-x_{j})^{2}\leq
w_{\max}\sum\limits_{i<j}(x_{i}-x_{j})^{2}=w_{\max}\ x^{T}L(G_{comp})x$ (2.16)

where $G_{comp}$ is the complete undirected graph, hence

$L(G_{comp})=\begin{pmatrix}n-1&-1&\cdots&-1\\\ -1&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&-1\\\ -1&\cdots&-1&n-1\end{pmatrix}.$ (2.17)

One can check that $x^{T}L(G_{comp})x\leq\lambda_{n}(L(G_{comp}))=n$, so
$\lambda_{n}\leq nw_{\max}$. ∎

####  2.3.1 Heat Kernel Distance process

Let $P$ be a probability distribution on
$\mathcal{G}_{n}\times\mathcal{G}_{n}$. Let $T$ be a positive real number and
recall that $\mathcal{F}_{HKD}:=\\{D_{t},\ t\in[0,T]\\}$, with $D_{t}$ defined
in (2.6). We will study the centered and rescaled empirical process
$\\{G_{N}D_{t},\ t\in I\\}=\\{\sqrt{N}(P_{N}-P)D_{t},\ t\in I\\}$, where
$P_{N}$ is the empirical measure
$N^{-1}\sum_{i}\delta_{(G_{i},{G^{\prime}}_{i})}$ associated to a $N$-sample
$((G_{1},{G^{\prime}}_{1}),\dots,(G_{N},{G^{\prime}}_{N}))$ drawn under $P$.
We denote by $\langle\cdot,\cdot\rangle$ the standard scalar product in
$\mathbb{R}^{n}$.

Let us prove that the HKD process is uniformly bounded and Lipschitz-
continuous. Note that these properties will exactly be the conditions required
in Section 3 to prove the statistical results on general empirical processes.

######  Proposition 2.1.

For all $G,G^{\prime}$ in $\mathcal{G}_{n}$, denote by $(\lambda_{k})_{1\leq
k\leq n}$ and $(\phi_{k})_{1\leq k\leq n}$ (resp.
$(\lambda_{l}^{\prime})_{1\leq l\leq n}$ and $(\phi_{l}^{\prime})_{1\leq l\leq
n}$) the eigenvalues and orthonormal eigenvectors of $L(G)$ (resp.
$L(G^{\prime})$). Remember that we always choose $\phi_{1}$ and
$\phi_{1}^{\prime}$ equal to the vector with all entries equal to
$1/\sqrt{n}$. Then, the application ${t\to D_{t}((G,G^{\prime}))}$ verifies :

  1. 1.

For all $t\in[0,T]$, $D_{t}((G,G^{\prime}))$ can be written in terms of the
eigen-elements of $L(G)$ and $L(G^{\prime})$ :

$D_{t}((G,G^{\prime}))=\left(\sum\limits_{k,l=2}^{n}(e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}})^{2}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}\right)^{1/2}.$
(2.18)

  2. 2.

For all $t\in[0,T]$, $D_{t}((G,G^{\prime}))\leq\sqrt{n}$.

  3. 3.

$t\mapsto D_{t}((G,G^{\prime}))$ is $(n^{3/2}w_{\max})$-Lipschitz continuous
on $[0,T]$.

###### Proof of Proposition 2.1.

Let $G,G^{\prime}$ be in $\mathcal{G}_{n}$. We start by proving (2.18). This
is done through the following computation :

$\displaystyle\|e^{-tL}-e^{-tL^{\prime}}\|_{F}^{2}$ $\displaystyle=$
$\displaystyle\|e^{-tL}\|_{F}^{2}+\|e^{-tL^{\prime}}\|_{F}^{2}-2Tr\left(e^{-tL}e^{-tL^{\prime}}\right)$
$\displaystyle=$
$\displaystyle\sum\limits_{k=1}^{n}e^{-2t\lambda_{k}}\|\phi_{k}\phi_{k}^{T}\|_{F}^{2}+\sum\limits_{l=1}^{n}e^{-2t\lambda^{\prime}_{l}}\|\phi^{\prime}_{l}{\phi^{\prime}_{l}}^{T}\|_{F}^{2}-2\sum\limits_{k,l=1}^{n}e^{-t\lambda_{k}}e^{-t\lambda^{\prime}_{l}}Tr\left(\phi_{k}\phi_{k}^{T}\phi^{\prime}_{l}{\phi^{\prime}_{l}}^{T}\right)$

One can prove that
$Tr\left(\phi_{k}\phi_{k}^{T}\phi^{\prime}_{l}{\phi^{\prime}_{l}}^{T}\right)=\langle\phi_{k},\phi^{\prime}_{l}\rangle^{2}$.
Similarly
$\|\phi_{k}\phi_{k}^{T}\|_{F}^{2}=\|\phi_{k}\|_{2}^{2}=\sum\limits_{l=1}^{n}\langle\phi_{k},\phi^{\prime}_{l}\rangle^{2}$
and
$\|\phi^{\prime}_{l}{\phi^{\prime}_{l}}^{T}\|_{F}^{2}=\|\phi^{\prime}_{l}\|_{2}^{2}=\sum\limits_{k=1}^{n}\langle\phi_{k},\phi^{\prime}_{l}\rangle^{2}$.
So

$\displaystyle\|e^{-tL}-e^{-tL^{\prime}}\|_{F}^{2}=$
$\displaystyle\sum\limits_{k,l=1}^{n}(e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}})^{2}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}.$

The sums can start at $k=2$ and $l=2$ thanks to the facts that
$\lambda_{1}=\lambda^{\prime}_{1}=0$ and $\phi_{1}=\phi^{\prime}_{1}$ combined
with the orthogonality of the eigenvectors families. This finishes the proof
of (2.18).

Bounding all terms $(e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}})^{2}$ by 1
in (2.18), and using the orthonormality of the eigenvectors families yields to
$D_{t}((G,G^{\prime}))\leq\sqrt{n}$.

We now prove the Lipschitz result. One can check that $t\to
D_{t}((G,G^{\prime}))$ is $\mathcal{C}^{1}$ on $[0,T]$. The case
$G=G^{\prime}$ is easily dealt with. Assume now that $G\neq G^{\prime}$. For
all $t\in(0,T]$,

$\displaystyle\left|\frac{d}{dt}D_{t}((G,G^{\prime}))\right|=$
$\displaystyle\left|\frac{-\sum\limits_{k,l=2}^{n}(e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}})(\lambda_{k}e^{-t\lambda_{k}}-\lambda_{l}^{\prime}e^{-t\lambda_{l}^{\prime}})\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}}{\left(\sum\limits_{k,l=2}^{n}(e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}})^{2}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}\right)^{1/2}}\right|$
$\displaystyle\leq$
$\displaystyle\sum\limits_{k,l=2}^{n}\frac{\left|e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}}\right|\left|\lambda_{k}e^{-t\lambda_{k}}-\lambda_{l}^{\prime}e^{-t\lambda_{l}^{\prime}}\right|\left|\langle\phi_{k},\phi_{l}^{\prime}\rangle\right|^{2}}{\left(\sum\limits_{i,j=2}^{n}(e^{-t\lambda_{i}}-e^{-t\lambda_{j}^{\prime}})^{2}\langle\phi_{i},\phi_{j}^{\prime}\rangle^{2}\right)^{1/2}}$
$\displaystyle\leq$
$\displaystyle\sqrt{\sum\limits_{k,l=2}^{n}\left|\lambda_{k}e^{-t\lambda_{k}}-\lambda_{l}^{\prime}e^{-t\lambda_{l}^{\prime}}\right|^{2}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}}\sqrt{\frac{\sum\limits_{k,l=2}^{n}(e^{-t\lambda_{k}}-e^{-t\lambda_{l}^{\prime}})^{2}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}}{\sum\limits_{i,j=2}^{n}(e^{-t\lambda_{i}}-e^{-t\lambda_{j}^{\prime}})^{2}\langle\phi_{i},\phi_{j}^{\prime}\rangle^{2}}}$
$\displaystyle=$
$\displaystyle\sqrt{\sum\limits_{k,l=2}^{n}\left|\lambda_{k}e^{-t\lambda_{k}}-\lambda_{l}^{\prime}e^{-t\lambda_{l}^{\prime}}\right|^{2}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}},$

where the last inequality comes from the Cauchy-Schwarz inequality.

According to the mean value theorem, for all $k$ and $l$

$\left|\lambda_{k}e^{-t\lambda_{k}}-\lambda_{l}^{\prime}e^{-t\lambda_{l}^{\prime}}\right|\leq\left|\lambda_{k}-\lambda_{l}^{\prime}\right|\leq\Lambda_{\max}.$

Hence, for all $t\in(0,T]$

$\left|\frac{d}{dt}D_{t}((G,G^{\prime}))\right|\leq\Lambda_{\max}\sqrt{\sum\limits_{k,l=2}^{n}\langle\phi_{k},\phi_{l}^{\prime}\rangle^{2}}\leq
nw_{\max}\sqrt{n}=n^{3/2}w_{\max}.$

This finishes the proof. ∎

We now state the statistical results concerning the HKD process.

######  Theorem 2.2.

For all probability distribution $P$ on
$\mathcal{G}_{n}\times\mathcal{G}_{n}$, the family $\mathcal{F}_{HKD}$ is
$P$-Donsker. That is, the process $\\{\sqrt{N}(P_{N}-P)D_{t},\ t\in[0,T]\\}$
converges weakly to $\mathbb{G}$ in $\mathcal{C}([0,T])$, where
$\mathbb{G}=\\{\mathbb{G}_{t},\ t\in[0,T]\\}$ is a zero mean Gaussian process
with covariance function $\kappa(t,s)=P(D_{t}D_{s})-PD_{t}.PD_{s}$.

This theorem can be seen as a functional Central Limit Theorem for the HKD
process. As Section 3.3 will show, it validates the construction of consistent
confidence bands and consistent two-sample tests. To strengthen this result,
we provide information about the speed of convergence happening in Theorem
2.2.

######  Theorem 2.3.

For all probability distribution $P$ on
$\mathcal{G}_{n}\times\mathcal{G}_{n}$, the process
$\\{G_{N}D_{t},\ t\in I\\}$ admits a Gaussian approximation with rate
$r_{N}=N^{-1/7}\log N^{9/14}$.

Note that this rate is independent of the graph size $n$.

These two theorems are direct applications of Theorem 3.1 and Theorem 3.2 from
Section 3, combined with Proposition 2.1.

####  2.3.2 Heat Persistence Distance process

In this section we work in $\mathcal{G}^{n}$. We denote by $n(G)$ the size of
a graph $G$. Let $P$ be a probability distribution on
$\mathcal{G}^{n}\times\mathcal{G}^{n}$. Let $T$ be a positive real number and
recall that $\mathcal{F}_{HPD}:=\\{H_{t},\ t\in[0,T]\\}$, with $H_{t}$ defined
in (2.8). Again, we study the empirical process $\\{G_{N}H_{t},\ t\in
I\\}=\\{\sqrt{N}(P_{N}-P)H_{t},\ t\in I\\}$, where $P_{N}$ is the empirical
measure $N^{-1}\sum_{i}\delta_{(G_{i},{G^{\prime}}_{i})}$ associated to a
$N$-sample $((G_{1},{G^{\prime}}_{1}),\dots,(G_{N},{G^{\prime}}_{N}))$ drawn
under $P$.

Before stating the statistical results on HPD processes, we prove that they
are uniformly bounded and Lipschitz-continuous.

######  Proposition 2.2.

For all $G,G^{\prime}$ in $\mathcal{G}^{n}$, the application $t\to
H_{t}((G,G^{\prime}))$ verifies :

  1. 1.

For all $t\in[0,T]$ , $0\leq H_{t}((G,G^{\prime}))\leq 1$.

  2. 2.

$t\mapsto H_{t}((G,G^{\prime}))$ is $(2nw_{\max})$-Lipschitz-continuous on
$[0,T]$.

###### Proof of Proposition 2.2.

Recall that for all $G\in\mathcal{G}^{n}$, and for a vertex $i$,

$h_{t}(G)(i)=\sum\limits_{k=1}^{n(G)}e^{-t\lambda_{k}}\phi_{k}(i)^{2},$

where $(\lambda_{k})_{1\leq k\leq n(G)}$ and $(\phi_{k})_{1\leq k\leq n(G)}$
are the eigenvalues and orthonormal eigenvectors of $L(G)$. Hence $0\leq
h_{t}(G)(i)\leq 1$ for all $i$, meaning that all points in the diagram
$Dg(G,h_{t}(G))$ are contained in $[0,1]^{2}$. So from the definition of the
Bottleneck distance, $0\leq H_{t}((G,G^{\prime}))\leq 1$, for all
$G,G^{\prime}$ in $\mathcal{G}^{n}$ and for all $t\in[0,T]$.

Let us now compute the first derivative of $h_{t}(G)(i)$ :

$\frac{d}{dt}h_{t}(G)(i)=-\sum\limits_{k=1}^{n(G)}\lambda_{k}e^{-t\lambda_{k}}\phi_{k}(i)^{2}.$

Its absolute value is upper-bounded by $\lambda_{n(G)}$, the largest
eigenvalue of $L(G)$. From Lemma 2.1, we have $\lambda_{n(G)}\leq
n(G)w_{\max}$. Hence, $t\to h_{t}(G)$ is $(nw_{\max})$-Lipschitz continuous on
$[0,T]$. To conclude, we come back to the definition of the HPD in terms of
distance between persistence diagrams. Applying the triangular inequality to
the Bottleneck distance gives for all $G,G^{\prime}\in\mathcal{G}^{n}$, for
all $t,t^{\prime}\in[0,T]$, and for all diagram construction $Dg$,

$\displaystyle\left|d_{B}(Dg(G,h_{t}(G)),Dg(G^{\prime},h_{t}(G^{\prime})))-d_{B}(Dg(G,h_{t^{\prime}}(G)),Dg(G^{\prime},h_{t^{\prime}}(G^{\prime})))\right|$
$\displaystyle\leq$ $\displaystyle
d_{B}(Dg(G,h_{t}(G)),Dg(G,h_{t^{\prime}}(G)))+d_{B}(Dg(G^{\prime},h_{t}(G^{\prime})),Dg(G^{\prime},h_{t^{\prime}}(G^{\prime}))).$

Similarly, the same inequality but with maxima over the diagram constructions
holds :

$\displaystyle\left|H_{t}((G,G^{\prime}))-H_{t^{\prime}}((G,G^{\prime}))\right|$
$\displaystyle\leq$ $\displaystyle\max
d_{B}(Dg(G,h_{t}(G)),Dg(G,h_{t^{\prime}}(G)))+\max
d_{B}(Dg(G^{\prime},h_{t}(G^{\prime})),Dg(G^{\prime},h_{t^{\prime}}(G^{\prime}))).$

Applying Theorem 2.1 and using the Lispchitz continuity of the HKS yields

$\displaystyle\left|H_{t}((G,G^{\prime}))-H_{t^{\prime}}((G,G^{\prime}))\right|$
$\displaystyle\leq$
$\displaystyle\|h_{t}(G)-h_{t^{\prime}}(G)\|_{\infty}+\|h_{t}(G^{\prime})-h_{t^{\prime}}(G^{\prime})\|_{\infty}$
$\displaystyle\leq$ $\displaystyle 2nw_{\max}\left|t-t^{\prime}\right|.$

∎

We now state the statistical results concerning the HPD process.

######  Theorem 2.4.

For all probability distribution $P$ on
$\mathcal{G}^{n}\times\mathcal{G}^{n}$, the family $\mathcal{F}_{HPD}$ is
$P$-Donsker. Thus, the process $\\{\sqrt{N}(P_{N}-P)H_{t},\ t\in[0,T]\\}$
converges weakly to $\mathbb{G}$ in $\mathcal{C}([0,T])$, where
$\mathbb{G}=\\{\mathbb{G}_{t},\ t\in[0,T]\\}$ is a zero mean Gaussian process
with covariance function $\kappa(t,s)=P(H_{t}H_{s})-PH_{t}.PH_{s}$.

######  Theorem 2.5.

For all probability distribution $P$ on
$\mathcal{G}^{n}\times\mathcal{G}^{n}$, the process
$\\{G_{N}H_{t},\ t\in I\\}$ admits a Gaussian approximation with rate
$r_{N}=N^{-1/7}\log N^{9/14}$.

As previously, these results are a direct consequence of Theorem 3.1 and
Theorem 3.2, as well as Proposition 2.2.

##  3 General continuous empirical processes

In this section, we generalize the previous results on empirical processes. We
properly introduce the general framework for continuous empirical processes,
then show that uniform boundedness and Lipschitz-continuity implies a
functional Central Limit Theorem, as well as a Gaussian approximation.
Finally, we derive consequences on the construction of confidence bands and
two-sample tests.

###  3.1 Background and definitions

Let $I$ be a compact interval of $\mathbb{R}$ and $\mathcal{C}(I)$ the space
of continuous real-valued functions on $I$ endowed with the metric induced by
the uniform norm : $\left\|h\right\|_{\infty}=\sup_{t\in I}|h(t)|$. Consider a
measurable space $(\mathbb{X},\mathcal{X})$. For all measure $Q$ on
$(\mathbb{X},\mathcal{X})$ and all measurable function
$g:\mathbb{X}\to\mathbb{R}$, we denote the integral of $g$ with respect to $Q$
by $Qg:=\int_{\mathbb{X}}g(x)dQ(x)$. Consider a probability measure $P$ on
$(\mathbb{X},\mathcal{X})$ and $\mathcal{F}:=\\{f_{t},\ t\in I\\}$, a family
of measurable real-valued functions on $\mathbb{X}$ indexed by $I$. For all
$x\in\mathbb{X}$, define $f(x)$ as the function $t\to f_{t}(x)$, and assume
that $f(x)\in\mathcal{C}(I)$. Therefore, given a random variable $X$ with
distribution $P$, one can equivalently see $\\{f_{t}(X),\ t\in I\\}$ either as
a random process or as $f(X)$ a random variable in $\mathcal{C}(I)$.

Given an i.i.d sample $X_{1},\dots,X_{N}$ drawn under $P$, we are interested
in the statistical properties of the mean function $N^{-1}\sum_{i}f(X_{i})$
and its centered and scaled version
$N^{-1/2}\left(\sum_{i}f(X_{i})-Pf\right)$. Equivalently, one can study the
empirical processes $\\{P_{N}f_{t},\ t\in I\\}$ and $\\{G_{N}f_{t},\ t\in
I\\}$, where $P_{N}=N^{-1}\sum_{i}\delta_{X_{i}}$, and
$G_{N}=\sqrt{N}(P_{N}-P)$. In the following, we see random processes and
random functions as the same objects.

When studying the statistical properties of $\mathcal{F}$, one might be
interested in a functional version of the Central Limit Theorem. This
corresponds to the concept of Donsker families.

######  Definition 3.1 (P-Donsker).

For $P$ a probability measure on $(\mathbb{X},\mathcal{X})$, the family
$\mathcal{F}$ is called $P$-Donsker if the process $\\{G_{N}f_{t},\ t\in I\\}$
converges in distribution to the centered Gaussian process
$\\{\mathbb{G}_{t},\ t\in I\\}$ with covariance function $\kappa$ defined as
$\kappa_{s,t}:=Pf_{t}f_{s}-Pf_{t}Pf_{s}$ for all $t,s\in I$..

Here convergence in distribution means weak convergence in the space
$\mathcal{C}(I)$. That is, for all continuous bounded function
$h:\mathcal{C}(I)\to\mathbb{R}$,
$\lim_{N\to\infty}\mathbb{E}\left[h\left(G_{N}f\right)\right]=\mathbb{E}\left[h\left(\mathbb{G}\right)\right]$,
where the expectation on the left-hand side is taken over the distribution of
the sample $X_{1},\dots,X_{N}$, and the one on the right-hand side is taken
over the distribution of the Gaussian process $\mathbb{G}$.

Going further into the statistical analysis of $\mathcal{F}$, one might want
to assess the speed at which $\\{G_{N}f_{t},\ t\in I\\}$ converges in
distribution to $\\{\mathbb{G}_{t},\ t\in I\\}$. This can be done by proving
Gaussian approximation results.

######  Definition 3.2 (Gaussian Approximation).

Let $(r_{N})_{N\geq 1}$ be a vanishing sequence of positive real numbers. We
say that the process $\\{G_{N}f_{t},\ t\in I\\}$ admits a Gaussian
approximation with rate $r_{N}$, if for all $\lambda>1$, there exists a
constant $C$ such that for all $N\geq 1$ one can construct on the same
probability space both the sample $X_{1},\dots,X_{N}$ and a version
$\mathbb{G}^{(N)}$ of the Gaussian process $\mathbb{G}$ verifying

$\mathbb{P}\left(\left\|G_{N}f-\mathbb{G}^{(N)}\right\|_{\infty}>C.r_{N}\right)\leq
N^{-\lambda}.$ (3.1)

Note that if $\\{G_{N}f_{t},\ t\in I\\}$ admits a Gaussian approximation with
rate $r_{N}$, applying the Borel-Cantelli Lemma would yield to
$\left\|G_{N}f-\mathbb{G}^{(N)}\right\|_{\infty}=O(r_{N})$ almost surely.

#### Assumptions

A few assumptions on $\mathcal{F}$ will be needed in the following. We present
the main ones here.

  * (L) -

There exits $k>0$ such that for all $x\in\mathbb{X}$ the function $t\to
f_{t}(x)$ is $k$-Lipschitz continuous on $I$, meaning that for all $t,s\in I$

$\left|f_{t}(x)-f_{s}(x)\right|\leq k\left|t-s\right|.$

  * (B) -

$\mathcal{F}$ is uniformly bounded. That is, there exists a constant $M>0$
such that for all $x\in\mathbb{X}$ and for all $t\in I$,

$\left|f_{t}(x)\right|\leq M.$

Note that the assumption (B) implies the existence of the covariance function
$\kappa$. It also implies that all moments of $f_{t}(X)$ are finite.

###  3.2 Donsker theorem and Gaussian approximation

In this section, we prove that assumptions (L) and (B) are sufficient to
obtain a Donsker theorem and a Gaussian approximation result. This will allow
us to derive, in the next section, more practical consequences of these
results. Namely, we will prove the validity of bootstrap methods to construct
consistent confidence bands and consistent two-sample tests. Note that these
results are very general and assume very little about the family
$\mathcal{F}$.

######  Theorem 3.1.

Assume that $\mathcal{F}$ verifies assumptions (L) and (B) .
Then $\mathcal{F}$ is $P$-Donsker.

###### Proof of Theorem 3.1.

We need to prove the weak convergence of $\\{G_{N}f_{t},\ t\in I\\}$ to the
centered gaussian process $\mathbb{G}$. To do so, remark that assumption (B)
ensures that second moments of $f(X)$ are finite. This allows applying the
multidimensional version of the Central Limit Theorem to all finite-
dimensional marginals of the random function $G_{N}f$. Hence, these finite-
dimensional marginals converge in distribution to those of $\mathbb{G}$. To
conclude to the weak convergence of the whole process, the sequence of random
functions $\\{G_{N}f,\ N\geq 1\\}$ needs to be tight. According to [3, Theorem
12.3], tightness of $\\{G_{N}f,\ N\geq 1\\}$ is implied by the following two
conditions :

  1. 1.

There exists $t_{0}\in I$, such that $\\{G_{N}f_{t_{0}},\ N\geq 1\\}$ is
tight.

  2. 2.

There exist $\gamma\geq 0$, $\alpha>1$ and a non-decreasing function
$\psi:I\to\mathbb{R}$ such that $\forall t,s\in I$, $\forall N\geq 1$,

$\mathbb{E}\left[\left|G_{N}f_{t}-G_{N}f_{s}\right|^{\gamma}\right]\leq\left|\psi(t)-\psi(s)\right|^{\alpha}.$
(3.2)

Let us start by proving point 1. Let $t_{0}$ be any point in $I$. Recall the
definition of tightness : $\\{G_{N}f_{t_{0}},\ N\geq 1\\}$ is tight if for all
$\eta>0$, there exists $\alpha>0$ such that for all $N\geq 1$

$P\left(\left|G_{N}f_{t_{0}}\right|\leq\alpha\right)>1-\eta.$ (3.3)

If $Var(f_{t_{0}}(X))=0$, $\\{G_{N}f_{t_{0}},\ N\geq 1\\}$ is tight, since for
all $N$, $G_{N}f_{t_{0}}=0$ $P$-a.s.. Otherwise, fix $\eta>0$. As the left
hand side in (3.3) is non-decreasing with respect to $\alpha$, we may as well
show that there exists $\alpha$ such that for all $N$,
$P\left(-\alpha<G_{N}f_{t_{0}}\leq\alpha\right)>1-\eta$. Following from (B) ,
$f_{t_{0}}(X)$ admits a third moment. Combined with the positive variance, we
can apply the Berry-Essen Theorem : if $F_{N}$ and $\phi$ denote the
cumulative distribution functions of, respectively, $G_{N}f_{t_{0}}$ and a
centered Gaussian variable of variance $Var(f_{t_{0}}(X))$, then there exists
a constant $C$ such that for all $\alpha\in\mathbb{R}$ and all $N\geq 1$

$\left|F_{N}(\alpha)-\phi(\alpha)\right|\leq\frac{C}{\sqrt{N}}.$ (3.4)

Now we take $N_{\eta}$ such that for all $N>N_{\eta}$, $C/\sqrt{N}<\eta/4$,
and we choose $\alpha_{\eta}$ such that $\phi(\alpha_{\eta})>1-\eta/4$. Then,
for all $N>N_{\eta}$

$\displaystyle P\left(-\alpha_{\eta}<G_{N}f_{t_{0}}\leq\alpha_{\eta}\right)$
$\displaystyle=$ $\displaystyle F_{N}(\alpha_{\eta})-F_{N}(-\alpha_{\eta})$
$\displaystyle=$
$\displaystyle\phi(\alpha_{\eta})-\phi(-\alpha_{\eta})+\left(F_{N}(\alpha_{\eta})-\phi(\alpha_{\eta})\right)-\left(F_{N}(-\alpha_{\eta})-\phi(-\alpha_{\eta})\right)$
$\displaystyle\geq$
$\displaystyle\phi(\alpha_{\eta})-\phi(-\alpha_{\eta})-2\frac{C}{\sqrt{N}}$
$\displaystyle=$ $\displaystyle 2\phi(\alpha_{\eta})-1-2\frac{C}{\sqrt{N}}$
$\displaystyle>$ $\displaystyle
2\left(1-\frac{\eta}{4}\right)-1-2\frac{\eta}{4}=1-\eta$

One can easily choose $\alpha>\alpha_{\eta}$ to extend the inequality
$P\left(-\alpha_{\eta}<G_{N}f_{t_{0}}\leq\alpha_{\eta}\right)>1-\eta$ to all
$N\geq 1$. This finishes the proof of point 1.

The proof of point 2, is a consequence of assumption (L) . For all $t,s\in I$,
one has the following inequality

$\displaystyle\mathbb{E}\left[\left|G_{N}f_{t}-G_{N}f_{s}\right|^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left|(f_{t}(X)-Pf_{t})-(f_{s}(X)-Pf_{s})\right|^{2}\right]$
$\displaystyle\leq$ $\displaystyle(2k\left|t-s\right|)^{2}$ $\displaystyle=$
$\displaystyle(2kt-2ks)^{2}.$

This proves point 2. with $\gamma=2$, $\alpha=2$ and $\psi(u)=2ku$, and
finishes the proof of Theorem 3.1. ∎

######  Theorem 3.2.

Assume that $\mathcal{F}$ verifies assumptions (L) and (B) .
Then $\\{G_{N}f_{t},\ t\in I\\}$ admits a Gaussian approximation with rate
$r_{N}=N^{-1/7}\log N^{9/14}$.

The proof of Theorem 3.2 is based on a result from [2]. They derive rates for
the Gaussian approximation of more general processes. They do so by
approaching the process by finite-dimensional marginals and applying a
multidimensional Gaussian approximation result. Controlling the covering
number of the family (or equivalently its metric entropy) allows them to
derive a good trade-off between a good approximation by the marginals and
keeping the dimension small enough to obtain good rates in the
multidimensional Gaussian approximation. For the sake of completeness, let us
recall their result.

Let $\mathcal{M}$ be the set of all measurable real-valued functions on
$(\mathbb{X},\mathcal{X})$. The authors work with a centered and scaled
empirical process $\\{G_{N}\tilde{f},\ \tilde{f}\in\tilde{\mathcal{F}}\\}$
indexed by a general family $\tilde{\mathcal{F}}\subset\mathcal{M}$, not
necessarily indexed by $I$. Their result shows that, under mild assumptions,
this process can be approached by a centered Gaussian process
$\tilde{\mathbb{G}}$ indexed by $\tilde{\mathcal{F}}$ with covariance
$\mathbb{E}\left[\tilde{\mathbb{G}}(\tilde{f})\tilde{\mathbb{G}}(\tilde{g})\right]=P\tilde{f}\tilde{g}-P\tilde{f}P\tilde{g}$,
for $\tilde{f},\tilde{g}\in\tilde{\mathcal{F}}$. Let us define the covering
number of $\tilde{\mathcal{F}}$. First, assume that there exists
$\tilde{F}\in\mathcal{M}$ such that for all $\tilde{f}\in\tilde{\mathcal{F}}$
and all $x\in\mathbb{X}$, $|\tilde{f}(x)|\leq\tilde{F}(x)$. We say that
$\tilde{F}$ is an envelope of $\tilde{\mathcal{F}}$. For all probability
measure $Q$ on $(\mathbb{X},\mathcal{X})$, we consider the semi-metric
$d_{Q}(g,h)^{2}=\int(g-h)^{2}dQ$, for $g,h\in\mathcal{M}$. Under $d_{Q}$, we
define the ball of radius $\delta>0$ centered in $h\in\mathcal{M}$ by
$B_{Q}(h,\delta):=\\{g\in\mathcal{M},\ d_{Q}(g,h)<\delta\\}$. We also define
$\tilde{F}_{Q}^{2}:=\int\tilde{F}^{2}dQ$. For $\delta>0$, let
$N(\tilde{\mathcal{F}},d_{Q},\delta)$ be the size of the smallest finite
subset $K\subset\mathcal{M}$ verifying that the union of balls
$B_{Q}(h,\delta)$ for $h$ in $K$ covers $\tilde{\mathcal{F}}$. Finally, we set
the covering number of $\tilde{\mathcal{F}}$ to be
$N(\tilde{\mathcal{F}},\delta):=\sup_{Q}N(\tilde{\mathcal{F}},d_{Q},\delta.\tilde{F}_{Q})$,
where the supremum is taken over all $Q$ such that $0<\tilde{F}_{Q}<\infty$.

Before stating the result of [2], consider these two basic assumptions.

  * (F.i)

For some $M>0$ and for all $\tilde{f}\in\tilde{\mathcal{F}}$,
$\sup_{x\in\mathbb{X}}|\tilde{f}(x)|\leq M/2$.

  * (F.ii)

The class $\tilde{\mathcal{F}}$ is point-wise measurable, i.e. there exists a
countable subclass $\tilde{\mathcal{F}}_{\infty}$ of $\tilde{\mathcal{F}}$
such that we can find for any function $\tilde{f}\in\tilde{\mathcal{F}}$ a
sequence of functions $\\{\tilde{f}_{m}\\}$ in $\tilde{\mathcal{F}}_{\infty}$
for which $\lim_{m\rightarrow\infty}\tilde{f}_{m}(x)=\tilde{f}(x)$ for all
$x\in\mathbb{X}$.

######  Proposition 3.1.

[2, Proposition 1.] Assume that $\tilde{\mathcal{F}}$ verifies (F.i) and
(F.ii). Take $\tilde{F}:=M/2$ as the envelope of $\tilde{\mathcal{F}}$.
Moreover, assume that there exist positive constants $c_{0}$ and $v_{0}$ such
that $N(\tilde{\mathcal{F}},\delta)\leq c_{0}\delta^{-v_{0}}$. Then, for each
$\lambda>1$ there is a constant $\rho(\lambda)$ such that for each $N$, one
can construct $X_{1},\dots,X_{N}$ and $\tilde{\mathbb{G}}^{(N)}$ such that

$\mathbb{P}\left(\sup\limits_{\tilde{f}\in\tilde{\mathcal{F}}}\left|G_{N}\tilde{f}-\tilde{\mathbb{G}}^{(N)}(\tilde{f})\right|>\rho(\lambda)N^{-\frac{1}{2+5v_{0}}}\log
N^{\frac{4+5v_{0}}{4+10v_{0}}}\right)\leq N^{-\lambda}.$ (3.5)

###### Proof of Theorem 3.2.

The proof consists in applying Proposition 3.5 to $\mathcal{F}$ with
$v_{0}=1$. Clearly, assumption (B) implies (F.i). Since the paths $t\to
f_{t}(x)$ are continuous for all $x\in\mathbb{X}$, taking
$\mathcal{F}_{\infty}=\\{f_{t},\ t\in I\cap\mathbb{Q}\\}$ gives (F.ii), where
$\mathbb{Q}$ is the set of rational numbers.

Let us prove the upper-bound on the covering number. Let $L$ be the length of
$I$. From (L) , we know that there exists a constant $k$, such that $t\to
f_{t}(x)$ is $k$-lipschitz, for all $x$. Considering a regular grid on $I$ ,
$\min I=t_{0}<t_{1}<\dots<t_{q}=\max I$. For all $t\in I$ there exists a
integer $j$ such that for all $x\in\mathbb{X}$,
$\left|f_{t}(x)-f_{t_{j}}(x)\right|\leq k\left|t-t_{j}\right|\leq kL/q$.
Hence, for all probability measure $Q$, $d_{Q}(f_{t},f_{t_{j}})\leq kL/q$,
meaning that $N(\mathcal{F},d_{Q},kT/q)\leq q+1$. As this last inequality
stands for all $Q$, we can find a constant $c_{0}>0$ such that
$N(\mathcal{F},\delta)\leq c_{0}\delta^{-1}$, for all $\delta>0$. This
finishes the proof. ∎

###  3.3 Statistical consequences

#### Confidence Band

Let $c_{\alpha}:=\inf\\{u,\
\mathbb{P}\left(\left\|\mathbb{G}\right\|_{\infty}>u\right)\leq\alpha\\}$ be
the upper $\alpha$-quantile of the maximum of the Gaussian process. As a
consequence of Theorem 3.1, we have

$\lim_{N\to\infty}\mathbb{P}\left(\forall
t,Pf_{t}\in\left[P_{N}f_{t}-\frac{c_{\alpha}}{\sqrt{N}},P_{N}f_{t}+\frac{c_{\alpha}}{\sqrt{N}}\right]\right)\geq
1-\alpha.$ (3.6)

Unfortunately, as the distribution of $\mathbb{G}$ is unknown, $c_{\alpha}$
can not be directly computed. Instead, consider a bootstrap sample
$\hat{X}_{1},\dots,\hat{X}_{N}$ drawn under $P_{N}$ and let $\hat{P}_{N}$ be
its empirical probability measure. Consider the process
$\\{\hat{G}_{N}f_{t},t\in I\\}$ where $\hat{G}_{N}$ is the measure
$\sqrt{N}(\hat{P}_{N}-P_{N})$. Theorem 2.6 in [18] ensures that in the case
where $\mathcal{F}$ is $P$-Donsker, $\\{\hat{G}_{N}f_{t},t\in I\\}$ converges
weakly to $\mathbb{G}$, given the data. Let $\hat{c}_{\alpha}$ be the upper
$\alpha$-quantile of $\|\hat{G}_{N}f\|_{\infty}$ given the data. We can use
$\hat{c}_{\alpha}$ to design a consistent confidence band around $P_{N}f$ :

$\lim_{N\to\infty}\mathbb{P}\left(\forall
t,Pf_{t}\in\left[P_{N}f_{t}-\frac{\hat{c}_{\alpha}}{\sqrt{N}},P_{N}f_{t}+\frac{\hat{c}_{\alpha}}{\sqrt{N}}\right]\right))\geq
1-\alpha.$ (3.7)

Note that $\hat{c}_{\alpha}$ can be estimated with Monte-Carlo simulations, by
drawing as many bootstrap samples as we want.

#### Two Sample Tests

Consider the following setup. Let $P$ and $Q$ be two probability distributions
on $\mathbb{X}$, and assume we are given two independent iid samples
$(X_{1},\dots,X_{M})$ and $(Y_{1},\dots Y_{N})$, drawn under $P$ and $Q$,
respectively. We denote by $P_{M}$ and $Q_{N}$ the empirical measures. We
would like to test the null hypothesis $H_{0}:P=Q$ against the alternatives
$H_{1}:P\neq Q$, by using the family $\mathcal{F}$, assuming that it is
Donsker with respect to both $P$ and $Q$. We follow the approach described in
Section 3.7 of [31], and consider the following test statistic :

$D_{M,N}:=\sqrt{\frac{MN}{M+N}}\left\|P_{M}f-Q_{N}f\right\|_{\infty}.$ (3.8)

The strategy is to define a data-dependent threshold $\hat{c}_{M,N}(\alpha)$
and reject the null hypothesis whenever $D_{M,N}>\hat{c}_{M,N}(\alpha)$.
Consider the pooled data $(Z_{1},\dots
Z_{M+N})=(X_{1},\dots,X_{M},Y_{1},\dots,Y_{N})$, and its empirical measure
$H_{N+M}$. Let $(\hat{Z}_{1},\dots\hat{Z}_{M+N})$ be a bootstrap sample drawn
from $H_{M+N}$ and consider the bootstrap empirical measures

$\hat{P}_{M}=\frac{1}{M}\sum\limits_{i=1}^{M}\delta_{\hat{Z}_{i}}\qquad\textnormal{and}\qquad\hat{Q}_{N}=\frac{1}{N}\sum\limits_{i=1}^{N}\delta_{\hat{Z}_{M+i}}.$
(3.9)

We can define

$\hat{D}_{M,N}:=\sqrt{\frac{MN}{M+N}}\left\|\hat{P}_{M}f-\hat{Q}_{N}f\right\|_{\infty},$
(3.10)

as well as

$\hat{c}_{M,N}(\alpha)=\inf\left\\{t,\
\mathbb{P}\left(\left\|\hat{D}_{M,N}f\right\|_{\infty}>t\ \Big{|}\
Z_{1},\dots,Z_{N+M}\right)\leq\alpha\right\\},$ (3.11)

for $\alpha\in(0,1)$. Note that $\hat{c}_{M,N}(\alpha)$ can be estimated with
Monte-Carlo simulations. Using $\hat{c}_{M,N}(\alpha)$ as the threshold to
accept or reject $H_{0}$ leads to a consistent test.

######  Theorem 3.3 (Section 3.7.2, [31]).

Assume that $\mathcal{F}$ is Donsker with respect to both $P$ and $Q$ and that
$\left\|Pf\right\|_{\infty}$ and $\left\|Qf\right\|_{\infty}$ are finite.
Furthermore assume that $M/(M+N)\to\lambda\in(0,1)$. Then the test that
rejects $H_{0}$ whenever $D_{M,N}>\hat{c}_{M,N}(\alpha)$ is consistent, in the
sense that the asymptotic level is $\alpha$ and under any alternative
verifying $\left\|Pf-Qf\right\|_{\infty}>0$,
$\mathbb{P}\left(D_{M,N}>\hat{c}_{M,N}(\alpha)\right)\to 1$.

##  4 Experiments

We illustrate the construction of confidence bands and two-sample tests on
synthetic data sets of pairs of graphs. For that, we consider different random
graph models and combine them to create independent pairs of graphs. Let us
present the models.

###  4.1 Random graph models

#### Erdős-Rényi Model (ER)

[11] This model generates random graphs where each edge appears with
probability $p$, independently from all the others. It requires two
parameters: $n$ the graph size and $p$ the edge probability. Because of the
independence and their homogeneity, ER graphs are considered to have no
structure.

In our simulations, we take $n=50$ and $p=0.5$. Weights may be added by
assigning a uniform weight between 0 and 2 to each existing edge,
independently from all the others.

#### Stochastic Block Model (SBM)

[16] This model is a generalization of the ER model that introduces a block
structure. The $n$ nodes are clustered in $K$ groups $C_{1},\dots,C_{K}$, of
respective sizes $n_{1},\dots,n_{K}$. Edges appear independently from the
others, but with a probability depending on the groups : edge $\\{i,j\\}$
appears with probability $p_{k,l}$ when $i\in C_{k}$ and $j\in C_{l}$.

In our simulations, we take $K=2$, $n_{1}=n_{2}=25$, and
$p_{1,1}=p_{2,2}=0.75$, $p_{1,2}=p_{2,1}=0.25$. So graphs are composed of two
dense clusters, with few edges between them. Similarly to the ER model, we may
add random weights following the uniform distribution between 0 and 2.

#### Geometric Model (GM)

[24] Given a compact domain $U$ of $\mathbb{R}^{d}$ for some $d$, a graph is
generated from the GM by drawing $n$ points uniformly on $U$ and creating an
edge between two points if their Euclidean distance is smaller than a given
threshold. Here we choose a slight variation of this model by considering a
number $p\in[0,1]$ and creating the edges corresponding to the $\lfloor
p\binom{n}{2}\rfloor$ pairs of points with the smallest euclidean distances.

In our simulations, the compact domain $U$ is either $A_{\varepsilon}$ the
annulus in $\mathbb{R}^{2}$ with outer radius 1 and inner radius
$\varepsilon>0$ or $A_{0}$ the unit disk. We either fix $n=50$ or for each
graph, $n$ is drawn from a Poisson distribution of parameter 50. We take
$p=0.5$. To obtain weighted graphs, we may assign the weight $e^{-2d}$ to an
edge, where $d$ is the distance between the two points forming the edge.

We combine these models to generate pairs of independent graphs on which we
can compute HKD and HPD processes. We consider the pairs of independent ER
graphs (ER-ER) and the pairs containing one ER graph and one SBM graph (ER-
SBM). For these distributions, the groups’ composition is known and nodes are
treated independently among groups. Thus, we can consider that we know an NC
between graphs. As a result, we can compute both HKD and HPD processes.
Similarly, we consider pairs of independent geometric random graphs: Disk-Disk
and Disk-Annulus. However, as nodes in these models correspond to random
points, there is no default NC. Hence, only HPD processes are computed.

###  4.2 Simulation results

####  4.2.1 Confidence bands

In this section, we compute confidence bands under the different models of
pairs of graphs defined above. For each model, we draw a sample
$(G_{1},G^{\prime}_{1}),\dots,(G_{N},G^{\prime}_{N})$ with $N=100$. We compute
the mean process, that is $t\to N^{-1}\sum_{i}D_{t}((G_{i},G^{\prime}_{i}))$
or $t\to N^{-1}\sum_{i}H_{t}((G_{i},G^{\prime}_{i}))$ and compute a confidence
band of level 99% around this empirical mean using the bootstrap method
presented in Section 3.3. Computations are done with 1000 bootstrap samples.
Results are shown in Figure 1, 2 and 3, where solid lines represent empirical
means and transparent areas represent confidence bands.

Remark that confidence bands around HKD processes (Figure 1) seem to be
narrower than those around HPD processes (Figure 2). Therefore, users should
rather use HKD processes whenever NC’s are available. Nonetheless, the
versatility of HPD processes does not totally reduce their efficiency. As
Figure 3 indicates, HPD empirical means seem to be able to discriminate
between the different distributions. These observations will be confirmed in
the next section, where the performances of the two-sample tests are
investigated.

(a) Unweighted.

(b) Weighted.

Figure 1: Confidence band around the mean HKD processes with ER-ER (red) and
ER-SBM (blue) distributions. Figure 2: Confidence band around the mean HPD
processes with unweighted ER-ER (red) and ER-SBM (blue) distributions.

(a) Unweighted,
fixed size.

(b) Weighted,
fixed size.

(c) Unweighted,
random size.

Figure 3: Confidence band around the mean HPD processes with Disk-Disk (red)
and Disk-Annulus (blue) distributions.

####  4.2.2 Two-sample tests

##### Levels and Powers

Simulations are run to evaluate the performances of the two-sample tests using
HKD and HPD processes, see Figure 4 and Figure 5. In all tests, the desired
level is set to 0.05, and computations are done with 1000 bootstrap samples.
Figures 4(a) and 5(a) illustrate that the asymptotic levels of the tests
correspond to the set level. On the other hand, Figures 4(b) and 5(b)
illustrate that the powers tend to 1 when sample sizes increase, indicating
that the tests manage to distinguish between the different distributions.

(a) Levels estimations.

(b) Powers estimations.

Figure 4: Performances of the two-sample test using HKD processes. The level
estimations are made with weighted ER-ER distributions, while the power
estimations are made with weighted ER-ER and ER-SBM distributions. The samples
size varies in [25,50,100,200]. Estimations are done by repeating 400
independent tests. Vertical lines represent $95\%$-confidence intervals of the
estimations.

(a) Levels estimations.

(b) Powers estimations.

Figure 5: Performances of the two-sample test using HPD processes. The level
estimations are made with Disk-Disk distributions, while the power estimations
are made with Disk-Disk and Disk-Annulus distributions. Graphs are unweighted
and with fixed size. The samples size varies in [20,40,80,120]. Estimations
are done by repeating 100 independent tests. Vertical lines represent
$95\%$-confidence intervals of the estimations.

##### Comparison with the Neyman-Pearson Test

The Neyman-Pearson test [22] is an optimal testing procedure, in the sense
that it is the test with the highest power for a given level. But being based
on likelihood ratios, it rarely is computable. Its performances are determined
by the total variation distance (TV) between the two distributions. If we
consider the case of distributions that depend on the sample size, the speed
at which their TV tends to 0 determines if the Neyman-Pearson test will
asymptotically be able to distinguish between the two distributions. For ER
models, independence of the edges allows to theoretically compute bounds of
the TV and determine the phase transition. Let us take a real number $p$ such
that $0<p<1$. And for each sample size $N\geq 1$ consider the parameters
$p_{0}(N)$ and $p_{1}(N)$, such that both converge to $p$. Following [20,
Section 13.1.], we can show that as long as $\left|p_{0}(N)-p_{1}(N)\right|\gg
N^{-1/2}$, the Neyman-Pearson test asymptotically distinguishes between the
distributions of independent pairs of $ER(n,p_{0}(N))$ and independent pairs
of $ER(n,p_{1}(N))$ when using $N$-samples. Figure 6 shows that for large
enough sample sizes, our two-sample test based on HKD processes distinguishes
between ER models up to $\left|p_{0}(N)-p_{1}(N)\right|=C\log N/\sqrt{N}$,
with $C=0.01$. The tests are computed with a set level of 0.05 and with 1000
bootstrap samples.

Figure 6: Evolution of the power of the HKD two-sample test when
$\left|p_{0}(N)-p_{1}(N)\right|=C\log N/\sqrt{N}$. The samples size $N$ varies
in [20,50,100,200,300]. Estimations are done by repeating 100 independent
tests. Vertical lines represent $95\%$-confidence power intervals of the
estimations.

##  5 Conclusion

We proposed two multiscale comparisons of graphs using heat diffusion
processes, namely the HKD and HPD. The first one requires the assumption of
equal graph sizes and a known NC, while the second one is free of these
assumptions. The multiscale approach solves the problem of choosing an
informative diffusion time. We proposed to use these processes to analyze data
sets of pairs of graphs and were able to design consistent confidence bands
and two-sample tests. The methods are supported by theoretical results: the
HKD and HPD families are Donsker, meaning that the processes verify a
functional Central Limit Theorem. Moreover, the processes admit Gaussian
approximations with rates that are independent of the graph sizes. These
results are very general and can be applied to other processes under mild
assumptions. Essentially, the processes are required to be uniformly bounded
and Lipschitz-continuous. Moreover, the performances of our methods were
evaluated by simulations on synthetic data sets. We showed that the two-sample
tests were able to distinguish between Erdős-Rényi and SBM graphs, as well as
between geometric graphs sampled on different domains. On Erdős-Rényi models
with parameters depending on the sample size, the tests were still
distinguishing between the different distributions, even when working close to
the phase transition of Neyman-Pearson tests.

As future work, we would like to apply these methods to real-world data sets
and extend them to be able to perform learning tasks, e.g. clustering,
classification, or change point detection. Extensions to be able to deal with
data sets of graphs by opposition to pairs of graphs should also be developed
to broaden the application spectrum. On the theoretical side, studying the
interplay between graph sizes and sample sizes could be a first step toward
non-asymptotic methods to analyze data sets of graphs. Nonetheless, we believe
that the introduction of the HKD and HPD processes has the potential to bring
innovative and statistically founded ways to analyze data sets of graphs.
Moreover, the theoretical results presented in this work being very general,
our methods could be extended to other fields.

## Acknowledgments

I am thankful to Frédéric Chazal111Inria Saclay and Pascal
Massart222Université Paris-Saclay for the valuable discussions and advice on
this work. I also want to thank the people from Datashapefootnote 1 for the
enriching conversations.

## References

  * [1] Ali, W., Rito, T., Reinert, G., Sun, F., and Deane, C. M. Alignment-free protein interaction network comparison.  Bioinformatics 30, 17 (2014), i430–i437.
  * [2] Berthet, P., and Mason, D. M. Revisiting two strong approximation results of dudley and philipp.  In High dimensional probability. Institute of Mathematical Statistics, 2006, pp. 155–172.
  * [3] Billingsley, P. Convergence of probability measures.
  * [4] Carrière, M., Chazal, F., Ike, Y., Lacombe, T., Royer, M., and Umeda, Y. Perslay: a neural network layer for persistence diagrams and new graph topological signatures.  In International Conference on Artificial Intelligence and Statistics (2020), PMLR, pp. 2786–2796.
  * [5] Chazal, F., De Silva, V., Glisse, M., and Oudot, S. The structure and stability of persistence modules.  Springer, 2016.
  * [6] Cohen-Steiner, D., Edelsbrunner, H., and Harer, J. Extending persistence using
poincaré and lefschetz duality.  Foundations of Computational Mathematics 9, 1
(2009), 79–103.

  * [7] Coifman, R. R., and Hirn, M. J. Diffusion maps for changing data.  Applied and computational harmonic analysis 36, 1 (2014), 79–107.
  * [8] Coifman, R. R., and Lafon, S. Diffusion maps.  Applied and computational harmonic analysis 21, 1 (2006), 5–30.
  * [9] Edelsbrunner, H., and Harer, J. Computational topology: an introduction.  American Mathematical Soc., 2010.
  * [10] Emmert-Streib, F., Dehmer, M., and Shi, Y. Fifty years of graph matching, network alignment and network comparison.  Information sciences 346 (2016), 180–197.
  * [11] Erdos, P., Rényi, A., et al. On the evolution of random graphs.  Publ. Math. Inst. Hung. Acad. Sci 5, 1 (1960), 17–60.
  * [12] Faisal, F. E., Newaz, K., Chaney, J. L., Li, J., Emrich, S. J., Clark, P. L., and Milenković, T. Grafene: Graphlet-based alignment-free network approach integrates 3d structural and sequence (residue order) data to improve protein structural comparison.  Scientific reports 7, 1 (2017), 1–15.
  * [13] Fiedler, M. An estimate for the nonstochastic eigenvalues of doubly stochastic matrices.  Linear algebra and its applications 214 (1995), 133–143.
  * [14] Gera, R., Alonso, L., Crawford, B., House, J., Mendez-Bermudez, J., Knuth, T., and Miller, R. Identifying network structure similarity using spectral graph theory.  Applied network science 3, 1 (2018), 1–15.
  * [15] Hammond, D. K., Gur, Y., and Johnson, C. R. Graph diffusion distance: A difference measure for weighted graphs based on the graph laplacian exponential kernel.  In 2013 IEEE Global Conference on Signal and Information Processing (2013), IEEE, pp. 419–422.
  * [16] Holland, P. W., Laskey, K. B., and Leinhardt, S. Stochastic blockmodels: First steps.  Social networks 5, 2 (1983), 109–137.
  * [17] Hu, N., Rustamov, R. M., and Guibas, L. Stable and informative spectral signatures for graph matching.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (2014), pp. 2305–2312.
  * [18] Kosorok, M. R. Introduction to empirical processes.  Introduction to Empirical Processes and Semiparametric Inference (2008).
  * [19] Koutra, D., Vogelstein, J. T., and Faloutsos, C. Deltacon: A principled massive-graph similarity function.  In Proceedings of the 2013 SIAM International Conference on Data Mining (2013), SIAM, pp. 162–170.
  * [20] Lehmann, E. L., and Romano, J. P. Testing statistical hypotheses.  Springer Science & Business Media, 2006.
  * [21] Maria, C., Boissonnat, J.-D., Glisse, M., and Yvinec, M. The gudhi library: Simplicial complexes and persistent homology.  In International congress on mathematical software (2014), Springer, pp. 167–174.
  * [22] Neyman, J., and Pearson, E. S. Ix. on the problem of the most efficient tests of statistical hypotheses.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 231, 694-706 (1933), 289–337.
  * [23] Oudot, S. Y. Persistence theory: from quiver representations to data analysis, vol. 209.  American Mathematical Society Providence, 2015.
  * [24] Penrose, M., et al. Random geometric graphs, vol. 5.  Oxford university press, 2003.
  * [25] Pržulj, N. Biological network comparison using graphlet degree distribution.  Bioinformatics 23, 2 (2007), e177–e183.
  * [26] Pržulj, N., Corneil, D. G., and Jurisica, I. Modeling interactome: scale-free or geometric?  Bioinformatics 20, 18 (2004), 3508–3515.
  * [27] Soundarajan, S., Eliassi-Rad, T., and Gallagher, B. A guide to selecting a network similarity method.  In Proceedings of the 2014 Siam international conference on data mining (2014), SIAM, pp. 1037–1045.
  * [28] Sun, J., Ovsjanikov, M., and Guibas, L. A concise and provably informative multi-scale signature based on heat diffusion.  In Computer graphics forum (2009), vol. 28, Wiley Online Library, pp. 1383–1392.
  * [29] Tantardini, M., Ieva, F., Tajoli, L., and Piccardi, C. Comparing methods for comparing networks.  Scientific reports 9, 1 (2019), 1–19.
  * [30] Tsitsulin, A., Mottin, D., Karras, P., Bronstein, A., and Müller, E. Netlsd: hearing the shape of a graph.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (2018), pp. 2347–2356.
  * [31] Van Der Vaart, A. W., and Wellner, J. A. Weak convergence.  In Weak convergence and empirical processes. Springer, 1996.
  * [32] Wilson, R. C., and Zhu, P. A study of graph spectra for comparing graphs and trees.  Pattern Recognition 41, 9 (2008), 2833–2841.
  * [33] Yaveroğlu, Ö. N., Malod-Dognin, N., Davis, D., Levnajic, Z., Janjic, V., Karapandza, R., Stojmirovic, A., and Pržulj, N. Revealing the hidden language of complex networks.  Scientific reports 4, 1 (2014), 1–9.