††thanks: Correspondence should be addressed to P.S.

]Department of Psychology & Tsinghua Laboratory of Brain and Intelligence,
Tsinghua University, Beijing, 100084, China. Also at ]Laboratory of Advanced
Computing and Storage, Central Research Institute, 2012 Laboratories, Huawei
Technologies Co. Ltd., Beijing, 100084, China.

]Institut de Mathématiques d’Orsay, Institut de Mathématiques d’Orsay, 14 Rue
Joliot-Curie, 91190 Gif-sur-Yvette, France.

]Department of Computer Science, University College London, London, WC1E 6AE,
UK.

]Laboratory of Advanced Computing and Storage, Central Research Institute,
2012 Laboratories, Huawei Technologies Co. Ltd., Beijing, 100084, China.

]Department of Psychology & Tsinghua Brain and Intelligence Lab, Tsinghua
University, Beijing, 100084, China.

# Encoding, decoding, and causality between complex networks

Yang Tian<EMAIL_ADDRESS>[  [  Hedong Hou
<EMAIL_ADDRESS>[  Guangzheng Xu<EMAIL_ADDRESS>[  Ziyang Zhang
<EMAIL_ADDRESS>[  Pei Sun<EMAIL_ADDRESS>[

###### Abstract

Complex networks are common in physics, biology, computer science, and social
science. Quantifying the relations (e.g., similarity or difference) between
complex networks paves the way for understanding the latent information shared
across networks. However, fundamental metrics of relations, such as
information divergence, mutual information, Fisher information, and causality,
are not well-defined between complex networks. As a compromise, commonly used
strategies (e.g., network embedding, matching, and kernel approaches) measure
network relations in data-driven ways. These approaches are computation-
oriented and inapplicable to analytic derivations in mathematics and physics.
To resolve these issues, we present a theory to derive an optimal
characterization of network topological properties. Our theory shows that a
network can be fully represented by a Gaussian variable defined by a function
of the discrete Schrödinger operator (i.e., the Laplacian), which
simultaneously satisfies network-topology-dependent smoothness and maximum
entropy properties. Based on this characterization, we can analytically
measure diverse relations between complex networks in terms of topology
properties. As illustrations, we primarily show how to define encoding (e.g.,
information divergence and mutual information), decoding (e.g., Fisher
information), and causality (e.g., transfer entropy and Granger causality)
between complex networks. We validate our framework on representative complex
networks (e.g., evolutionary random network models, protein-protein
interaction network, and chemical compound networks), and demonstrate that a
series of science and engineering challenges (e.g., network evolution,
clustering, and classification) can be tackled from a new perspective. A
computationally efficient implementation of our theory is released as an open-
source multi-platform toolbox.

††preprint: APS/123-QED

##  I Introduction

Complex networks exist universally across different disciplines [1]. Important
topics in physics (e.g., quantum system characterization [2, 3, 4] and non-
equilibrium dynamics analysis [5, 6, 7]), biology (e.g., brain [8, 9, 10, 11,
12], metabolic [13, 14, 15], and protein [16, 17, 18] networks analysis),
computer science (e.g., internet analysis [19, 20, 21] and information
tracking [22, 23]), and social science (e.g., scientific community analysis
[24, 25, 26] and opinion formation modelling [27, 28, 29]) would benefit from
studying complex networks [1] and, therefore, promote the rapid development of
network theories.

However, critical challenges to network theories persistently arise due to the
increasingly diverse application needs [1]. Among these challenges, a
fundamental yet intractable one concerns how to quantify the relations (e.g.,
similarity) between different complex networks [30]. To date, mainstream
metrics of network relations are developed in the contexts of network
embedding, matching, and kernel, three computation-oriented and data-driven
perspectives [30, 31]. Comprehensive reviews of these three perspectives can
be found in Refs. [30, 32], Refs. [33, 34, 35], and Refs. [36, 37],
respectively. In general, _embedding-based_ approaches embed networks into
low-dimensional metric spaces (e.g., future vectors) following preset rules to
calculate distances between networks or classify networks [30, 32, 31].
Therefore, these approaches critically depend on embedding rule designs and
may lack universal principles to guarantee their generalization capacities
[30]. _Matching-based_ approaches, such as exact [38] and inexact [36]
matching with and without edge-preserving properties, search for node mappings
between networks to realize optimal matching and measure similarity [34].
These approaches essentially deal with a kind of quadratic programming
problems [33, 39, 40] that are NP-hard [33] and require relaxations of problem
constraints to find approximate solutions (e.g., see Refs. [41, 42, 43] for
approximation). _Kernel-based_ methods evaluate the similarity between
networks by decomposing them into series of atomic substructures (e.g.,
graphlets [44], random walks [45], shortest paths [46], and cycles [47]) and
measuring kernel value among these substructures (i.e., counting the number of
shared substructures) [36, 37]. While these substructures generally reflect
network topology properties, they are essentially handcrafted [31], i.e.,
extracted by certain manually defined functions, and may imply extremely high-
dimensional, sparse, and non-smooth representations of large-scale networks
with poor generalization capacities [48]. In sum, while embedding-, matching-,
and kernel-based approaches have been extensively tested on empirical data
(e.g., on neural data [49, 50, 37, 51]), they are inevitably limited by
computational complexity (e.g., matching-based) or the dependence on empirical
choice of network features (e.g., embedding-based) and kernel functions (e.g.,
kernel-based) [49]. Even in cases where these methods are computationally
optimal, they may still be unsatisfactory because they do not derive intrinsic
relations between complex networks analytically and universally.

Analytic metrics of network relations are indispensable for developing
mathematics and statistic physics of complex networks [52] but remain largely
unknown yet. Certainly, one can simplify the distance between networks as the
Kolmogorov–Smirnov statistic between their degree distributions (e.g., see
discussions in [50]). However, two networks with the same degree distribution
can still have distinct topology properties because edge wring rules are not
completely determined by the degree distribution.

To suggest a possible way to define analytic and universal network relations,
here we develop a theory of an optimal characterization of complex networks
that simultaneously ensures smoothness (for optimal reflection of network
topology [53, 54, 55, 56]) and maximum entropy (for better support of
information-theoretical analysis [57]) properties in Secs. II-III. The derived
characterizations turn out to be specific Gaussian variables defined by the
functions of discrete Schrödinger operators of complex networks. Based on this
result, we can define analytic relation metrics (e.g., information divergence
[57], mutual information [57], Fisher information [57], and causality [58])
between complex networks in Sec. IV and explore their potential generalization
in Sec. V. In Secs. VI-VII, we demonstrate our approach on representative
networks (e.g., random network models, a protein-protein interaction network,
and chemical compound networks) to define encoding, decoding, and causal
analyses between them, which raise potential insights on related science and
engineering questions. In Sec. VIII, we discuss the advantages and
deficiencies of our theory and suggest future directions for improvements. An
open-source toolbox of our work is provided in Ref. [59].

##  II Question definition

An appropriate question definition is more important than the answer itself.
To suggest a potential direction, we consider following questions:

  * (I)

How to develop an analytic and universal characterization of network topology
that is free of subjective selection of topological properties and
computational optimization problems?

  * (II)

How to enable the characterization derived in question (I) to define analytic
metrics of network relation (e.g., information divergence [57], mutual
information [57], Fisher information [57], and causality [58]) without further
constraints?

Clearly, question (II) requires a probabilistic solution (e.g., define a
network as a random variable) of question (I). This idea inspires us to
consider a mapping $\phi:V\rightarrow\Omega$ between a network
$\mathcal{G}\left(V,E\right)$ without self-loops and a probability space
$\left(\Omega,\mathcal{F},\mathcal{P}\right)$ with $\Omega=\mathbb{R}$. Here
$V$ and $E$ denote the node and edge sets of network $\mathcal{G}$,
respectively. Function $\phi$ defines a random variable
$\mathcal{X}_{\phi}=\left(X_{\phi}\left(1\right),\ldots,X_{\phi}\left(n\right)\right)$
distributed on node set $V$, where
$X_{\phi}\left(i\right)=\phi\left(v_{i}\right)$ and $n=|V|$ (see Fig. 1 for
illustrations).

To properly reflect the network topology of $\mathcal{G}$ by
$\mathcal{X}_{\phi}$, we need to consider the _smoothness_ of mapping $\phi$
on $\mathcal{G}$ measured by the _discrete $2$-Dirichlet form_ of $\phi$ [54]

$\displaystyle\mathcal{S}\left(\phi\right)=\frac{1}{2}\sum_{v_{i}\in
V}\sum_{\left(v_{i},v_{j}\right)\in
E}\left(\frac{\partial\phi}{\partial\left(v_{i},v_{j}\right)}\bigg{|}_{v_{i}}\right)^{2},$
(1)

where $\left(v_{i},v_{j}\right)$ denotes the edge between nodes $v_{i}$ and
$v_{j}$. In Eq. (1), the edge derivative of $\phi$ with respect to edge
$\left(v_{i},v_{j}\right)$ at node $v_{i}$ is defined as [54]

$\displaystyle\frac{\partial\phi}{\partial\left(v_{i},v_{j}\right)}\bigg{|}_{v_{i}}=\sqrt{W_{i,j}}\left(X_{\phi}\left(j\right)-X_{\phi}\left(i\right)\right),$
(2)

where $W$ is the non-negative weighted adjacent matrix of $\mathcal{G}$. The
smaller $\mathcal{S}\left(\phi\right)$ in Eq. (1) is, the smoother $\phi$ is
on $\mathcal{G}$. To understand why the smoothness of $\phi$ matters in
defining $\mathcal{X}_{\phi}$ to reflect the network topology of
$\mathcal{G}$, we need to consider the combinatorial graph Laplacian $L$, a
kind of _discrete Schrödinger operator_ [60], of $\mathcal{G}$

$\displaystyle
L=\operatorname{diag}\left(\left[\operatorname{deg}\left(v_{1}\right),\ldots,\operatorname{deg}\left(v_{n}\right)\right]\right)-W,$
(3)

where $\operatorname{diag}\left(\cdot\right)$ generates a diagonal matrix and
operator $\operatorname{deg}\left(\cdot\right)$ measures node degree. Note
that $\operatorname{deg}\left(v_{i}\right)=\sum_{j=1}^{n}W_{ij}$ in a weighted
network. Laplacian $L$ captures key topology information of network
$\mathcal{G}$ (e.g., connected components, random walks on network, and latent
Laplace-Beltrami operator [61]), which has been extensively used in spectral
graph theory [61] and graph signal theory [54]. The first connection between
Laplacian $L$ and the smoothness of $\phi$ is well-known [54]

$\displaystyle\mathcal{S}\left(\phi\right)={\mathcal{X}_{\phi}}^{T}L\mathcal{X}_{\phi},$
(4)

which suggests that the smoothness of $\phi$ can be defined by Laplacian $L$
(see Fig. 1). The second connection is derived from the Courant-Fischer
theorem [62]. Letting $\\{\psi_{1},\ldots,\psi_{n}\\}$ and
$\\{\lambda_{1},\ldots,\lambda_{n}\\}$ be the eigenvectors and eigenvalues of
Laplacian $L$, one can use the Rayleigh quotient to derive that
$\lambda_{1}=\min\\{{\mathcal{X}_{\phi}}^{T}L\mathcal{X}_{\phi}\mid\mathcal{X}_{\phi}\in\mathbb{R}^{n},\;\|\mathcal{X}_{\phi}\|_{2}=1\\}$
and
$\lambda_{i}=\min\\{{\mathcal{X}_{\phi}}^{T}L\mathcal{X}_{\phi}\mid\mathcal{X}_{\phi}\in\mathbb{R}^{n},\;\|\mathcal{X}_{\phi}\|_{2}=1,\;\mathcal{X}_{\phi}\perp\operatorname{span}\left(\psi_{1},\ldots,\psi_{n}\right)\\}$
(here $2\leq i\leq n$) [54]. This property demonstrates that the smoothness of
$\phi$ is related to the eigenvectors and eigenvalues of Laplacian $L$.
Eigenvectors with smaller eigenvalues imply a smoother mapping $\phi$. Taken
together, the smoothness of $\phi$ matters in our analysis because it is
closely related to the topology information conveyed by Laplacian $L$. To
enable variable $\mathcal{X}_{\phi}$ to represent network $\mathcal{G}$, we
expect that the smoothness of $\phi$ is completely determined by the topology
properties of $\mathcal{G}$.

Figure 1:  General ideas of solving questions (I-II). We begin with defining a
random variable $\mathcal{X}_{\phi}$ on network $\mathcal{G}$ based on mapping
$\phi$. An ideal definition of $\mathcal{X}_{\phi}$ is expected to make the
smoothness of mapping $\phi$ be completely determined by network topology.
Meanwhile, it is expected to ensure the maximum entropy property of variable
$\mathcal{X}_{\phi}$. According to our derivations, we suggest
$\mathcal{X}_{\phi}\sim\mathcal{N}\left(\mathbf{0},L+\frac{1}{n}J\right)$ as a
possible solution.

To properly measure information quantities (e.g., mutual information) between
complex networks applying $\mathcal{X}_{\phi}$, the upper bounds of these
quantities in $\mathcal{X}_{\phi}$ should not be too small. Otherwise, these
quantities may be easily covered by noises in empirical data due to their
small orders of magnitude. To precisely understand why large upper bounds are
necessary, we consider the mutual information
$\mathcal{I}\left(\mathcal{X}_{\phi}^{i};\mathcal{X}_{\phi}^{j}\right)$
between two networks $\mathcal{G}_{i}$ and $\mathcal{G}_{j}$ as an instance.
It is trivial that
$\mathcal{I}\left(\mathcal{X}_{\phi}^{i};\mathcal{X}_{\phi}^{j}\right)\leq\min\\{\mathcal{H}\left(\mathcal{X}_{\phi}^{i}\right),\mathcal{H}\left(\mathcal{X}_{\phi}^{j}\right)\\}$,
where $\mathcal{H}\left(\cdot\right)$ denotes the Shannon entropy [57].
Therefore,
$\mathcal{I}\left(\mathcal{X}_{\phi}^{i};\mathcal{X}_{\phi}^{j}\right)$ may be
frequently clouded by noises or directly vanish given small quantities of
$\mathcal{H}\left(\mathcal{X}_{\phi}^{i}\right)$ and
$\mathcal{H}\left(\mathcal{X}_{\phi}^{j}\right)$. In the present study, we
primarily focus on the Shannon entropy because extensive information upper
bounds are related to it [57]. This idea inspires us to consider maximum
entropy distribution problem [57] while defining mapping $\phi$ (see Fig. 1).

In summary, an appropriate solution of questions (I-II) requires us to
consider the smoothness and maximum entropy properties of mapping $\phi$.
Below, we suggest a potential solution.

##  III Gaussian variable defined by the function of discrete Schrödinger
operator

To avoid that $\mathcal{S}\left(\phi\right)$ in Eq. (1) diverges, random
variable $\mathcal{X}_{\phi}$ is expected to have finite $1$-st and $2$-nd
moments on each dimension. In our theory, we suggest a possible scheme

$\displaystyle\mathbb{E}\left(\mathcal{X}_{\phi}\right)$
$\displaystyle=\mathbf{0},$ (5)
$\displaystyle\mathbb{D}\left(\mathcal{X}_{\phi}\right)$
$\displaystyle=\Sigma\in\mathbb{R}^{n\times n},$ (6)

where $\mathbf{0}=\left(0,\ldots,0\right)$ is a vector of zeros,
$\mathbb{E}\left(\cdot\right)$ denotes the $1$-st moment, and
$\mathbb{D}\left(\cdot\right)$ denotes the $2$-nd moment. Given Eqs. (5-6), we
can reformulate Eq. (1) as

$\displaystyle\mathbb{E}\left(\mathcal{S}\left(\phi\right)\right)$
$\displaystyle=\sum_{\left(v_{i},v_{j}\right)\in
E}\sqrt{W_{i,j}}\left(\Sigma_{ii}+\Sigma_{jj}-2\Sigma_{ij}\right),$ (7)

where $W$, the weighted adjacent matrix, is predetermined by a given network.
In Eq. (7), matrix $\Sigma$ is the only one adjustable term. To ensure the
effectiveness of $\mathcal{X}_{\phi}$ in reflecting the topology information
of $\mathcal{G}$, we suggest to choose matrix $\Sigma$ as

$\displaystyle\Sigma=L+\frac{1}{n}J,$ (8)

where $J$ is an all-one matrix. The motivation of the above definition lies in
four aspects. First, although $L$ is a singular matrix, previous studies have
proven that $L+\frac{1}{n}J$ is invertible [63, 64]

$\displaystyle\left(L+\frac{1}{n}J\right)^{-1}=L^{\dagger}+\frac{1}{n}J,$ (9)

where $\dagger$ denotes the Moore–Penrose pseudoinverse that satisfy
$LL^{\dagger}L=L$ [65]. This property ensures the possibility to calculate
numerous quantities defined with $\Sigma^{-1}$ in Sec. IV. Second, Eqs. (8-9)
relate $\Sigma$ with $L^{\dagger}$ directly. The pseudoinverse Laplacian
$L^{\dagger}$ is the _reproducing kernel_ of
$\mathbb{H}\left(\mathcal{G}\right)$, the Hilbert space of real-valued
functions over the node set $f_{\mathcal{G}}:V\rightarrow\mathbb{R}^{n}$ whose
inner product is $\langle
f_{\mathcal{G}},g_{\mathcal{G}}\rangle={f_{\mathcal{G}}}^{T}Lg_{\mathcal{G}}$
[66, 67]. Because $L^{\dagger}$ is unique for
$\mathbb{H}\left(\mathcal{G}\right)$, we can confirm a unique
$\mathbb{H}\left(\mathcal{G}\right)$ given $L^{\dagger}$. This property lays
foundations for kernel tricks [68, 69, 70] on network $\mathcal{G}$ when
future studies explore machine learning tasks on random variable
$\mathcal{X}_{\phi}$ (e.g., see kernel tricks in causality analysis [71, 72]).
Third, Laplacian $L$ and its pseudoinverse $L^{\dagger}$ directly determine
various topology properties of $\mathcal{G}$ (e.g., network coherence [73],
node importance [74], and the number of spanning trees [61]). Therefore, Eqs.
(8-9) ensure the expressive ability of $\mathcal{X}_{\phi}$ about network
topology. Fourth, we can apply Eq. (4) to derive the quadratic form

$\displaystyle\mathbb{E}\left(\mathcal{S}\left(\phi\right)\right)$
$\displaystyle=\mathbb{E}\left(\mathcal{X}_{\phi}\right)^{T}L\mathbb{E}\left(\mathcal{X}_{\phi}\right)+\operatorname{tr}\left(L\Sigma\right),$
(10)
$\displaystyle=\operatorname{tr}\left[L\left(L+\frac{1}{n}J\right)\right],$
(11)

if Eq. (8) holds. If network $\mathcal{G}$ is connected (i.e., there is only
one zero eigenvalue in $\\{\lambda_{1},\ldots,\lambda_{n}\\}$), we can further
apply $LL^{\dagger}=L^{\dagger}L=I-\frac{1}{n}J$ [75, 64, 74], where $I$ is
the unit matrix, to derive

$\displaystyle\mathbb{E}\left(\mathcal{S}\left(\phi\right)\right)$
$\displaystyle=\operatorname{tr}\left[LL+L\left(I-L^{\dagger}L\right)\right],$
(12) $\displaystyle=\operatorname{tr}\left(LL+L-LL^{\dagger}L\right),$ (13)
$\displaystyle=\operatorname{tr}\left(LL\right).$ (14)

Eq. (11) and Eq. (14) suggest a benefit of Eq. (8) that we can control the
expected smoothness of $\phi$ on a network completely by Laplacian $L$ (see
Fig. 1).

To ensure the maximum entropy property, we need to analyze the maximum entropy
distribution problem. Considering a random variable
$\mathcal{X}_{\phi}\in\mathbb{R}^{n}$ with finite $1$-st and $2$-nd moments
defined in Eqs. (5-6), we can know that

$\displaystyle\mathcal{H}\left(\mathcal{X}_{\phi}\right)\leq\mathcal{H}\left(\mathcal{N}\left(\mathbf{0},\Sigma\right)\right),$
(15)

where $\mathcal{N}\left(\mathbf{0},\Sigma\right)$ is the $n$-dimensional
Gaussian distribution. Eq. (15) is derived from the fact that the maximum
entropy distribution defined on $\mathbb{R}^{n}$ with given $1$-st and $2$-nd
moments in Eqs. (5-6) is the Gaussian distribution [57]. Therefore, we define
random variable $\mathcal{X}_{\phi}$ as

$\displaystyle\mathcal{X}_{\phi}\sim\mathcal{N}\left(\mathbf{0},L+\frac{1}{n}J\right)$
(16)

to reflect the topology information of network $\mathcal{G}$ (see Fig. 1). In
practice, we can readily derive the accurate entropy value

$\displaystyle\mathcal{H}\left(\mathcal{X}_{\phi}\right)=\frac{n}{2}+\frac{n}{2}\ln\left(2\pi\right)+\frac{1}{2}\ln\left[\operatorname{det}\left(L+\frac{1}{n}J\right)\right],$
(17)

where $\operatorname{det}\left(\cdot\right)$ denotes the determinant.

In sum, an appropriate solution of questions (I-II) requires us to represent
network $\mathcal{G}\left(V,E\right)$ by a Gaussian variable in Eq. (13),
which ensures the topology-dependent smoothness and maximum entropy properties
of mapping $\phi$ (see Fig. 1). Such a variable is completely characterized by
a function of the discrete Schrödinger operator of the network [60], whose
precision matrix (the inverse of covariance matrix) is
$L^{\dagger}+\frac{1}{n}J$.

Surprisingly, we notice that our result in Eq. (13) is similar with the graph
signal characterization [76, 77] derived by factor analysis [78] and low-rank
models [79, 80] in the field of engineering, which states that a Gaussian
Markov random field representation (a special type of Gaussian variable)
improves graph Laplacian learning in practice [76, 77]. The main difference
between Refs. [76, 77] and our work lies in that they pre-set the covariance
matrix as $\Sigma=L^{\dagger}$ while we define $\Sigma=L+\frac{1}{n}J$. This
similarity suggests the validity of our ideas from the perspective of
computation practice. In Sec. VIII.3, we present a more comprehensive
comparison between our results and previous studies [76, 77] from the
perspective of graph signal.

##  IV Analytic metrics of network relations

After representing two complex networks $\mathcal{G}_{a}$ and
$\mathcal{G}_{b}$ by Gaussian variables
$\mathcal{X}_{\phi}^{a}=\left(X_{\phi}^{a}\left(1\right),\ldots,X_{\phi}^{a}\left(n\right)\right)\sim\mathcal{N}\left(\mathbf{0},L\left(a\right)+\frac{1}{n}J\right)$
and
$\mathcal{X}_{\phi}^{b}=\left(X_{\phi}^{b}\left(1\right),\ldots,X_{\phi}^{b}\left(n\right)\right)\sim\mathcal{N}\left(\mathbf{0},L\left(b\right)+\frac{1}{n}J\right)$,
we can develop analytic metrics of network relations, such as information
divergence, mutual information, Fisher information, and casuality (e.g.,
transfer entropy and Granger causality). Because the all-one matrix $J$
frequently occurs in our derivations, we make a convention that $J$ has the
same size with matrix $A$ whenever we calculate $A+J$. The notions of matrix
sizes are omitted for convenience.

Figure 2:  Conceptual illustrations of information divergence and mutual
information between complex networks. The information divergence from network
$\mathcal{G}_{a}$ to network $\mathcal{G}_{b}$, defined by
$\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)$ in
Eq. (20), is not necessarily equal to the information divergence from network
$\mathcal{G}_{b}$ to network $\mathcal{G}_{a}$ (upper parallel). The mutual
information between networks $\mathcal{G}_{a}$ and $\mathcal{G}_{b}$, defined
by $\mathcal{I}\left(\mathcal{X}_{\phi}^{a};\mathcal{X}_{\phi}^{b}\right)$ in
Eq. (21), can be understood as the shared part of entropy quantities
$\mathcal{H}\left(\mathcal{X}_{\phi}^{a}\right)$ and
$\mathcal{H}\left(\mathcal{X}_{\phi}^{b}\right)$ (bottom parallel).

###  IV.1 Encoding: information divergence and mutual information

For information divergence (or referred to as the Kullback–Leibler divergence
[57]), we can formulate it in a conventional form

$\displaystyle\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)$
$\displaystyle=\mathbb{E}_{\rho_{a}}\left[\log\left(\rho_{a}\right)-\log\left(\rho_{b}\right)\right],$
(18)

where $\rho_{a}$ and $\rho_{b}$ are probability densities of
$\mathcal{X}_{\phi}^{a}$ and $\mathcal{X}_{\phi}^{b}$, respectively (see Fig.
2). Because $\mathcal{X}_{\phi}^{a}$ and $\mathcal{X}_{\phi}^{b}$ are Gaussian
variables, we can derive

$\displaystyle\log\left(\rho_{a}\right)=$
$\displaystyle-\frac{1}{2}\log\left[\left(2\pi\right)^{n}\operatorname{det}\left(L\left(a\right)+\frac{1}{n}J\right)\right]$
$\displaystyle-\frac{1}{2}\left[\mathcal{X}_{\phi}^{a}\right]^{T}\left(L^{\dagger}\left(a\right)+\frac{1}{n}J\right)\mathcal{X}_{\phi}^{a},$
(19)

which readily leads to

$\displaystyle\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)=$
$\displaystyle\frac{1}{2}\Bigg{[}\operatorname{tr}\left[\left(L^{\dagger}\left(b\right)+\frac{1}{n}J\right)\left(L\left(a\right)+\frac{1}{n}J\right)\right]$
$\displaystyle-n+\ln\frac{\operatorname{det}\left(L\left(b\right)+\frac{1}{n}J\right)}{\operatorname{det}\left(L\left(a\right)+\frac{1}{n}J\right)}\Bigg{]}.$
(20)

Eq. (20) measures the directional difference between the topology of
$\mathcal{G}_{a}$ and $\mathcal{G}_{b}$. The difference is directional since
$\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)\neq\mathcal{D}\left(\mathcal{X}_{\phi}^{b}\|\mathcal{X}_{\phi}^{a}\right)$
(see Fig. 2).

For mutual information
$\mathcal{I}\left(\mathcal{X}_{\phi}^{a};\mathcal{X}_{\phi}^{b}\right)$ that
quantifies the topology information of network $\mathcal{G}_{a}$ encoded by
network $\mathcal{G}_{b}$, we can calculate (see Fig. 2)

$\displaystyle\mathcal{I}\left(\mathcal{X}_{\phi}^{a};\mathcal{X}_{\phi}^{b}\right)$
$\displaystyle=\mathcal{H}\left(\mathcal{X}_{\phi}^{a}\right)+\mathcal{H}\left(\mathcal{X}_{\phi}^{b}\right)-\mathcal{H}\left(\mathcal{X}_{\phi}^{a},\mathcal{X}_{\phi}^{b}\right),$
(21)

where $\mathcal{H}\left(\mathcal{X}_{\phi}^{a}\right)$ and
$\mathcal{H}\left(\mathcal{X}_{\phi}^{b}\right)$ can be measured based on Eq.
(17). A challenge in Eq. (21) lies in that
$\mathcal{H}\left(\mathcal{X}_{\phi}^{a},\mathcal{X}_{\phi}^{b}\right)$ is
non-trivial for analytic derivations unless variables
$X_{\phi}^{a}\left(i\right)$ and $X_{\phi}^{b}\left(j\right)$ are jointly
Gaussian (this enables
$\mathcal{H}\left(\mathcal{X}_{\phi}^{a},\mathcal{X}_{\phi}^{b}\right)$ to be
defined by Eq. (17) as well). In practice, we can generate samples of
$\mathcal{X}_{\phi}^{a}\sim\mathcal{N}\left(\mathbf{0},L\left(a\right)+\frac{1}{n}J\right)$
and
$\mathcal{X}_{\phi}^{b}\sim\mathcal{N}\left(\mathbf{0},L\left(b\right)+\frac{1}{n}J\right)$
by inverse transform sampling [81] to estimate
$\mathcal{H}\left(\mathcal{X}_{\phi}^{a},\mathcal{X}_{\phi}^{b}\right)$ using
the Kozachenko-Leonenko estimator of Shannon entropy [82, 83].

Figure 3:  Conceptual illustrations of Fisher information between complex
networks. The $\left(i,j\right)$ element in the Fisher information matrix,
denoted by $\mathcal{F}_{ij}\left(\Theta\right)$ in Eq. (22), can be
understood as the information that network $\mathcal{G}_{a}$ carries about the
parameter vector $\Theta$ controlled by network $\mathcal{G}_{b}$. The
information is contained in $\rho_{a}\left(\cdot\mid\Theta\right)$, the
probability distribution of variable $\mathcal{X}_{\phi}^{a}$ that is affected
by $\Theta$. The information quantity can be understood as the mean
sensitivity of $\rho_{a}\left(\cdot\mid\Theta\right)$ towards the variation of
$\Theta$.

###  IV.2 Decoding: Fisher information

For Fisher information, we assume that a parameter vector
$\Theta=\left(\theta_{1},\ldots,\theta_{k}\right)\in\mathbb{R}^{k}$ controlled
by network $\mathcal{G}_{b}$ can affect the discrete Schrödinger operator
$L\left(a\right)$ of network $\mathcal{G}_{a}$. Fisher information measures
how precisely we can decode the topology information of $\mathcal{G}_{a}$ from
$\mathcal{G}_{b}$ according to parameter vector $\Theta$. We denote
$\mathcal{X}_{\phi}^{a}\sim\mathcal{N}\left(\mathbf{0},L^{\dagger}\left(a\mid\Theta\right)\right)$
as the Gaussian variable given parameter vector $\Theta$ (see Fig. 3). Then we
can have a special form of Fisher information matrix depending on the
covariance matrix [84, 85]

$\displaystyle\mathcal{F}_{ij}\left(\Theta\right)=$
$\displaystyle-\int\rho_{a}\left(\chi;\Theta\right)\left\\{\frac{\partial^{2}}{\partial\theta_{i}\theta_{j}}\log\left[\rho_{a}\left(\chi\mid\Theta^{\prime}\right)\right]\Bigg{|}_{\Theta^{\prime}=\Theta}\right\\}\mathsf{d}\chi,$
(22) $\displaystyle=$
$\displaystyle\frac{1}{2}\operatorname{tr}\Bigg{[}\left(L^{\dagger}\left(a\mid\Theta\right)+\frac{1}{n}J\right)\frac{\partial
L\left(a\mid\Theta\right)}{\partial\theta_{i}}$
$\displaystyle\times\left(L^{\dagger}\left(a\mid\Theta\right)+\frac{1}{n}J\right)\frac{\partial
L\left(a\mid\Theta\right)}{\partial\theta_{j}}\Bigg{]},$ (23)

where $\rho_{a}\left(\cdot\mid\Theta\right)$ is the probability density of
$\mathcal{X}_{\phi}^{a}$ given $\Theta$ (see Fig. 3). We define

$\displaystyle\dfrac{\partial L\left(a;\Theta\right)}{\partial\theta_{i}}$
$\displaystyle=\left[\begin{array}[]{ccc}\dfrac{\partial
L_{11}\left(a\mid\Theta\right)}{\partial\theta_{i}}&\ldots&\dfrac{\partial
L_{1n}\left(a\mid\Theta\right)}{\partial\theta_{i}}\\\ \vdots&\ddots&\vdots\\\
\dfrac{\partial
L_{n1}\left(a\mid\Theta\right)}{\partial\theta_{i}}&\ldots&\dfrac{\partial
L_{nn}\left(a\mid\Theta\right)}{\partial\theta_{i}}\end{array}\right].$ (27)

The expectation vector does not occur in Eq. (23) because
$\mathcal{X}_{\phi}^{a}$ always has zero expectation value on each dimension.
In application, one can further calculate Fisher information quantity,
$\operatorname{tr}\left[\mathcal{F}\left(\Theta\right)\right]$, as a metric of
decoding precision.

###  IV.3 Causality: Granger causality and transfer entropy

To this point, we have analytically derived information divergence, mutual
information, and Fisher information between complex networks. These metrics
lay the foundation of encoding and decoding analyses on network ensembles.
Compared with these analyses, causality is more difficult to study between
networks because it is initially analyzed on time series [58]. Although
dynamic networks feature time domain evolution [86, 87, 88], most networks
lack a well-defined concept of time (e.g., networks may be static). To develop
applicable causality metrics for arbitrary networks, we explore possible
generalization of the mainstream causality metrics, such as transfer entropy
[89, 90, 91] and Granger causality [92, 93, 94, 95, 95], from time domain to
graph domain.

Let us begin with a classic form of Granger causality analyzed by regression
models. Our basic idea is to consider a random partition on network
$\mathcal{G}_{a}$ to divide it into two sub-networks,
$\mathcal{G}_{a}^{\circ}$ and $\mathcal{G}_{a}^{\star}$. This is equivalent to
dividing random variable
$\mathcal{X}_{\phi}^{a}=\left(X_{\phi}^{a}\left(1\right),\ldots,X_{\phi}^{a}\left(n\right)\right)$
into two sub-vectors of multivariate Gaussian random variables
$\mathcal{X}_{\phi}^{a,\circ}$ and $\mathcal{X}_{\phi}^{a,\star}$. Without
loss of generality, we set
$\mathcal{X}_{\phi}^{a,\circ}=\left(X_{\phi}^{a}\left(1\right),\ldots,X_{\phi}^{a}\left(k\right)\right)\sim\mathcal{N}\left(\mathbf{0},L\left(a,\circ\right)+\frac{1}{n}J\right)$
and
$\mathcal{X}_{\phi}^{a,\star}=\left(X_{\phi}^{a}\left(k+1\right),\ldots,X_{\phi}^{a}\left(n\right)\right)\sim\mathcal{N}\left(\mathbf{0},L\left(a,\star\right)+\frac{1}{n}J\right)$,
where we know

$\displaystyle L\left(a,\circ\right)$
$\displaystyle=\left[\begin{array}[]{ccc}L_{11}\left(a\right)&\ldots&L_{1k}\left(a\right)\\\
\vdots&\ddots&\vdots\\\
L_{k1}\left(a\right)&\ldots&L_{kk}\left(a\right)\end{array}\right],$ (31)
$\displaystyle L\left(a,\star\right)$
$\displaystyle=\left[\begin{array}[]{ccc}L_{\left(k+1\right)\left(k+1\right)}\left(a\right)&\ldots&L_{\left(k+1\right)n}\left(a\right)\\\
\vdots&\ddots&\vdots\\\
L_{n\left(k+1\right)}\left(a\right)&\ldots&L_{nn}\left(a\right)\end{array}\right],$
(35)

and we can represent $L\left(a\right)$ in a block matrix form
$L\left(a\right)=\left[\begin{array}[]{cc}L\left(a,\circ\right)&L\left(a,\bigtriangleup\right)\\\
L\left(a,\bigtriangledown\right)&L\left(a,\star\right)\end{array}\right]$.
Then, we use $\mathcal{X}_{\phi}^{a,\circ}$ to predict
$\mathcal{X}_{\phi}^{a,\star}$ with a linear model

$\displaystyle\mathcal{X}_{\phi}^{a,\star}$
$\displaystyle=\beta+\mathcal{X}_{\phi}^{a,\circ}A+\varepsilon,$ (36)

where $A$ denotes the regression coefficient matrix, $\beta$ is a constant
vector, and $\varepsilon$ measures regression residuals. Meanwhile, we can
also use $\mathcal{X}_{\phi}^{a,\circ}$ and $\mathcal{X}_{\phi}^{b}$ to
predict $\mathcal{X}_{\phi}^{a,\star}$

$\displaystyle\mathcal{X}_{\phi}^{a,\star}$
$\displaystyle=\beta^{\prime}+\left(\mathcal{X}_{\phi}^{a,\circ}\oplus\mathcal{X}_{\phi}^{b}\right)A^{\prime}+\varepsilon^{\prime},$
(37)

where we have applied notion $\oplus$ to denote the concatenation of two
vectors, i.e.,
$\mathcal{X}_{\phi}^{a,\circ}\oplus\mathcal{X}_{\phi}^{b}=\left(X_{\phi}^{a}\left(1\right),\ldots,X_{\phi}^{a}\left(k\right),X_{\phi}^{b}\left(1\right),\ldots,X_{\phi}^{b}\left(n\right)\right)$
(see Fig. 4). According to Refs. [96, 97], the ordinary least squares
regression for Eqs. (36-37) is suggested to minimize the determinant of
covariance matrix of residuals (referred to as the generalized variance). The
covariance matrices of residuals for Eqs. (36-37) are

$\displaystyle\Sigma\left(\varepsilon\right)$
$\displaystyle=\Sigma\left(\mathcal{X}_{\phi}^{a,\star}\big{|}\mathcal{X}_{\phi}^{a,\circ}\right),$
(38) $\displaystyle\Sigma\left(\varepsilon^{\prime}\right)$
$\displaystyle=\Sigma\left(\mathcal{X}_{\phi}^{a,\star}\big{|}\mathcal{X}_{\phi}^{a,\circ}\oplus\mathcal{X}_{\phi}^{b}\right).$
(39)

The Granger causality can be defined as the natural logarithmic ratio between
the determinant values of Eqs. (38-39) [96]

$\displaystyle\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
$\displaystyle=\Bigg{\langle}\ln\left[\frac{\operatorname{det}\left(\Sigma\left(\varepsilon\right)\right)}{\operatorname{det}\left(\Sigma\left(\varepsilon^{\prime}\right)\right)}\right]\Bigg{\rangle},$
(40)

where the average $\langle\cdot\rangle$ is implemented across multiple
randomly generated configurations of $\mathcal{X}_{\phi}^{a,\star}$ and
$\mathcal{X}_{\phi}^{a,\circ}$ (i.e., we can randomly select $h$
configurations of $\mathcal{X}_{\phi}^{a,\star}$ and
$\mathcal{X}_{\phi}^{a,\circ}$ to calculate $h$ values of
$\ln\left[\frac{\Sigma\left(\varepsilon\right)}{\Sigma\left(\varepsilon^{\prime}\right)}\right]$
and average across them to derive Eq. (40)). Please see Fig. 4 for
illustrations. Because $\mathcal{X}_{\phi}^{a,\star}$ and
$\mathcal{X}_{\phi}^{a,\circ}$ are jointly Gaussian [97, 96], we can readily
derive

$\displaystyle\Sigma\left(\varepsilon\right)=$ $\displaystyle
L\left(a,\circ\right)+\frac{1}{n}J-\left(L\left(a,\bigtriangleup\right)+\frac{1}{n}J\right)$
$\displaystyle\times\left(L\left(a,\star\right)+\frac{1}{n}J\right)^{-1}\left(L\left(a,\bigtriangledown\right)+\frac{1}{n}J\right).$
(41)

However, $\Sigma\left(\varepsilon^{\prime}\right)$ can not be analytically
derived by the discrete Schrödinger operators of networks unless we relax
conditions (e.g., let $\mathcal{X}_{\phi}^{a,\star}$,
$\mathcal{X}_{\phi}^{a,\circ}$, and $\mathcal{X}_{\phi}^{b}$ be jointly
Gaussian as well). Similar to the situation in Eq. (21), we suggest that one
can generate samples of $\mathcal{X}_{\phi}^{a,\star}$,
$\mathcal{X}_{\phi}^{a,\circ}$, and $\mathcal{X}_{\phi}^{b}$ by inverse
transform sampling [81] to estimate $\Sigma\left(\varepsilon^{\prime}\right)$
in practice. For convenience, we usually write Eq. (40)) as
$\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)=\Big{\langle}\mathcal{T}_{g}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)\Big{\rangle}$,
where
$\mathcal{T}_{g}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)=\ln\left[\frac{\operatorname{det}\left(\Sigma\left(\varepsilon\right)\right)}{\operatorname{det}\left(\Sigma\left(\varepsilon^{\prime}\right)\right)}\right]$.

Figure 4:  Conceptual illustrations of Granger causality and transfer entropy
between complex networks. To quantify the causal effects from network
$\mathcal{G}_{b}$ to network $\mathcal{G}_{a}$, we first do random partition
on network $\mathcal{G}_{a}$ to obtain two sub-networks. Then we predict the
topology properties of one sub-network based on another sub-network. The
prediction may involve with residuals (in terms of Granger causality) or
uncertainty (in terms of transfer entropy). The causal effects from network
$\mathcal{G}_{b}$ to network $\mathcal{G}_{a}$ are reflected by the reduction
of residuals and uncertainty after we include the information of
$\mathcal{G}_{b}$ into the prediction. By repeating random partition and
prediction, we can measure Granger causality and transfer entropy in terms of
the average reduction of residuals and uncertainty.

Then, we turn to defining transfer entropy. In general, we can formulate
transfer entropy based on conditional mutual information [57] (see Fig. 4)

$\displaystyle\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
$\displaystyle=$
$\displaystyle\Big{\langle}\mathcal{I}\left(\mathcal{X}_{\phi}^{a,\star};\mathcal{X}_{\phi}^{b}\big{|}\mathcal{X}_{\phi}^{a,\circ}\right)\Big{\rangle}$
(42) $\displaystyle=$
$\displaystyle\Big{\langle}\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\star}\big{|}\mathcal{X}_{\phi}^{a,\circ}\right)-\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\star}\big{|}\mathcal{X}_{\phi}^{a,\circ}\oplus\mathcal{X}_{\phi}^{b}\right)\Big{\rangle}.$
(43)

Similar to Eq. (41), we can derive

$\displaystyle\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\star}\big{|}\mathcal{X}_{\phi}^{a,\circ}\right)=$
$\displaystyle\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\star}\oplus\mathcal{X}_{\phi}^{a,\circ}\right)-\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\circ}\right),$
(44) $\displaystyle=$
$\displaystyle\frac{n-k}{2}+\frac{n-k}{2}\log\left(2\pi\right)$
$\displaystyle+\frac{1}{2}\log\left[\frac{\operatorname{det}\left(L\left(a\right)+\frac{1}{n}J\right)}{\operatorname{det}\left(L\left(a,\circ\right)+\frac{1}{n}J\right)}\right]$
(45)

because $\mathcal{X}_{\phi}^{a,\star}$ and $\mathcal{X}_{\phi}^{a,\circ}$ are
jointly Gaussian. Similar to Eq. (21) and
$\Sigma\left(\varepsilon^{\prime}\right)$, we can not derive a general
expression of
$\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\star}\big{|}\mathcal{X}_{\phi}^{a,\circ}\oplus\mathcal{X}_{\phi}^{b}\right)=\mathcal{H}\left(\mathcal{X}_{\phi}^{a},\mathcal{X}_{\phi}^{b}\right)-\mathcal{H}\left(\mathcal{X}_{\phi}^{a,\circ}\oplus\mathcal{X}_{\phi}^{b}\right)$.
In practice, we can resolve this challenge by inverse transform sampling [81]
the Kozachenko-Leonenko estimator of Shannon entropy [82, 83]. For
convenience, Eq. (43)) is usually written as
$\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)=\Big{\langle}\mathcal{T}_{t}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)\Big{\rangle}$,
where we define
$\mathcal{T}_{t}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)=\mathcal{I}\left(\mathcal{X}_{\phi}^{a,\star};\mathcal{X}_{\phi}^{b}\big{|}\mathcal{X}_{\phi}^{a,\circ}\right)$.

Please note that the causality metrics considered here, such as
$\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
and
$\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$,
should be referred to as apparent causality metrics according to Ref. [98]
(one can see similar apparent causality metrics for time series in Refs. [89,
90, 91, 99, 92, 93, 94, 95, 95]). To derive complete causality metrics that
reflect causal relations more precisely (e.g., enable $\mathcal{T}_{G}$ and
$\mathcal{T}_{T}$ approximate causal information flow [100]), one need to
consider
$\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\big{|}\mathcal{Y}\right)$
and
$\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\big{|}\mathcal{Y}\right)$
given a reference variable $\mathcal{Y}$. Because the details of introducing
$\mathcal{Y}$ have been comprehensively explored in Refs. [98, 101, 96, 102]
and these details do not imply critical challenges for mathematical
derivations, we no longer repeat their analyses here. One can combine Refs.
[98, 101, 96, 102] and our framework to readily derive
$\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\big{|}\mathcal{Y}\right)$
and
$\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\big{|}\mathcal{Y}\right)$
between networks.

##  V Generalization of analytic metrics

One may notice that our derivations of all analytic metrics in Sec. IV are
shown in a special case where $\mathcal{X}_{\phi}^{a}$ and
$\mathcal{X}_{\phi}^{b}$ are both $n$-dimensional, meaning that networks
$\mathcal{G}_{a}$ and $\mathcal{G}_{b}$ both contain $n$ nodes. This
limitation arises from the fact that we need to calculate
$\left(L^{\dagger}\left(b\right)+\frac{1}{n}J\right)\left(L\left(a\right)+\frac{1}{n}J\right)$
in information divergence. The definitions of mutual information, Fisher
information, Granger causality, and transfer entropy have no such a limitation
and are generally applicable to arbitrary cases.

Figure 5:  Conceptual illustrations of the generalization of information
divergence. When networks $\mathcal{G}_{a}$ and $\mathcal{G}_{b}$ have
different sizes (e.g., network $\mathcal{G}_{a}$ contains more nodes), the
information divergence defined in Eq. (20) can not be calculated directly.
Consequently, we need to measure the importance of nodes in $\mathcal{G}_{a}$
from the perspective of network topology. In Eq. (47), we introduce
$\mathcal{A}\left(\cdot,\mathcal{G}_{a}\right)$, the importance measurement
based on Laplacian energy, as a possible approach. After assigning the
importance of each node in $\mathcal{G}_{a}$, we can exclude the nodes with
relatively lower importance to make the filtered $\mathcal{G}_{a}$ have the
same size as $\mathcal{G}_{b}$. Information divergence from the filtered
$\mathcal{G}_{a}$ to $\mathcal{G}_{b}$ can be measured and treated as an
approximation of the information divergence from the original
$\mathcal{G}_{a}$ to $\mathcal{G}_{b}$. The rationality $\gamma$ of this
approximation can be measured as the fraction of the lost Laplacian energy in
the original Laplacian energy.

In practice, we frequently need to analyze relations between complex networks
with distinct sizes (number of nodes). To make our information divergence
applicable to these networks, we suggest a practical solution based on
Laplacian energy. Let us consider a case where $\mathcal{X}_{\phi}^{a}$ and
$\mathcal{X}_{\phi}^{b}$ are $m$-dimensional and $n$-dimensional,
respectively. Without loss of generality, we primarily discuss the case where
$m>n$. The Laplacian energy of network
$\mathcal{G}_{a}\left(V_{a},E_{a}\right)$ is defined as [103, 104, 105, 106]

$\displaystyle\mathcal{L}\left(\mathcal{G}_{a}\right)$
$\displaystyle=\sum_{i=1}^{m}\lambda_{i}^{2}=\sum_{i=1}^{m}\operatorname{deg}\left(v_{i}\right)^{2}+\sum_{j\neq
i}W_{ij}^{2}\left(a\right),$ (46)

where $\left(\lambda_{1},\ldots,\lambda_{m}\right)$ are eigenvalues of
$L\left(a\right)$ and $W\left(a\right)$ denotes the weighted adjacent matrix
of $\mathcal{G}_{a}$. Note that
$\operatorname{deg}\left(v_{i}\right)=\sum_{j=1}^{n}W_{ij}$. Based on Eq.
(46), we can measure the importance of each node $v_{i}$ in maintaining
topology properties of $\mathcal{G}_{a}$ by Laplacian centrality [104, 106]

$\displaystyle\mathcal{A}\left(v_{i},\mathcal{G}_{a}\right)$
$\displaystyle=\frac{\mathcal{L}\left(\mathcal{G}_{a}\right)-\mathcal{L}\left(\mathcal{G}_{a}/\\{v_{i}\\}\right)}{\mathcal{L}\left(\mathcal{G}_{a}\right)},$
(47)

where $\mathcal{G}_{a}/\\{v_{i}\\}$ means deleting node $v_{i}$ from network
$\mathcal{G}_{a}$. Note that we have
$\mathcal{L}\left(\mathcal{G}_{a}\right)\geq\mathcal{L}\left(\mathcal{G}_{a}/\\{v_{i}\\}\right)$,
where the equality holds if and only if $v_{i}$ is an isolate node (has no
influence on main topology properties of $\mathcal{G}_{a}$) [106, 104]. In
general, the Laplacian centrality of node $v_{i}$ is determined by the number
of walks it participates in $\mathcal{G}_{a}$, i.e., the number of closed
walks $\left(v_{i},\ldots,v_{i}\right)$, the number of non-closed walks
$\left(v_{i},\ldots\right)$ and $\left(\ldots,v_{i}\right)$ where $v_{i}$ is
one of the end nodes, the number of non-closed walks
$\left(\ldots,v_{i},\ldots\right)$ containing $v_{i}$ as a middle node. In
Ref. [106], it is suggested that Eq. (47) can be reformulated by analyzing
walks of length $2$ (contain $3$ nodes). Specifically, one can derive Eq. (48)
following Ref. [106]

$\displaystyle\mathcal{L}\left(\mathcal{G}_{a}\right)-\mathcal{L}\left(\mathcal{G}_{a}/\\{v_{i}\\}\right)$
$\displaystyle=4\eta_{c}\left(v_{i}\right)+2\eta_{e}\left(v_{i}\right)+2\eta_{m}\left(v_{i}\right),$
(48)

where $\eta_{c}\left(v_{i}\right)$ measures the number of closed walks
$\left(v_{i},v_{j},v_{i}\right)$, $\eta_{e}\left(v_{i}\right)$ measures the
number of non-closed walks $\left(v_{i},v_{j},v_{k}\right)$ and
$\left(v_{j},v_{k},v_{i}\right)$, $\eta_{m}\left(v_{i}\right)$ measures the
number of non-closed walks $\left(v_{j},v_{i},v_{k}\right)$. Applying the
weight adjacent matrices $W\left(a\right)$ and $M\left(a\right)$ (here
$M\left(a\right)$ is the weight adjacent matrix of network
$\mathcal{G}_{a}/\\{v_{i}\\}$), we suggest that Eq. (48) can be readily
calculated by Eq. (49)

$\displaystyle\mathcal{L}\left(\mathcal{G}_{a}\right)-\mathcal{L}\left(\mathcal{G}_{a}/\\{v_{i}\\}\right)$
$\displaystyle=$ $\displaystyle
4\left[W\left(a\right)^{2}\right]_{ii}+2\left[\sum_{j\neq
k}\left[W\left(a\right)^{2}\right]_{jk}-\sum_{j\neq
k}\left[M\left(a\right)^{2}\right]_{jk}\right].$ (49)

Note that
$\left[W\left(a\right)^{2}\right]_{jk}=\left[W\left(a\right)W\left(a\right)\right]_{jk}$,
the $\left(i,j\right)$ element in matrix $W\left(a\right)^{2}$, is not
equivalent to $W_{ij}^{2}\left(a\right)$, the second power of
$\left(i,j\right)$ element in matrix $W\left(a\right)$. Combining Eqs. (46-47)
and Eq. (49), we can measure the Laplacian centrality of each node in network
$\mathcal{G}_{a}$ and only keep $n$ nodes with relatively large Laplacian
centrality values. In other words, $m-n$ nodes are filtered because their
effects on topology properties of $\mathcal{G}_{a}$ are less significant (see
Fig. 5). We refer to the network after filtering as
$\widehat{\mathcal{G}}_{a}$ and denote its Laplacian and Gaussian variable as
$\widehat{L}\left(a\right)$ and $\widehat{\mathcal{X}}_{\phi}^{a}$,
respectively. Then we can approximatively calculate information divergence in
Eq. (20) by replacing $L\left(a\right)$ and $\mathcal{X}_{\phi}^{a}$ as
$\widehat{L}\left(a\right)$ and $\widehat{\mathcal{X}}_{\phi}^{a}$ (see Fig.
5). In the case where $m<n$, we can similarly deal with network
$\mathcal{G}_{b}$ following the above approach.

The rationality $\gamma$ of the above approximation can be measured based on
the loss of Laplacian energy. Taking the case where $m>n$ as an instance, we
define the rationality of approximating $\mathcal{G}_{a}$ by
$\widehat{\mathcal{G}}_{a}$ as (see Fig. 5)

$\displaystyle\gamma\left(\mathcal{L}\left(\mathcal{G}_{a}\right),\mathcal{L}\left(\widehat{\mathcal{G}}_{a}\right)\right)$
$\displaystyle=\frac{\mathcal{L}\left(\widehat{\mathcal{G}}_{a}\right)}{\mathcal{L}\left(\mathcal{G}_{a}\right)}.$
(50)

To this point, we have presented analytic metrics of network relations from
the perspectives of encoding, decoding, and causal analyses in Sec. IV. We
have also explored their possible generalization in Sec. V. Below, we validate
our approach on representative complex networks to define encoding, decoding,
and causal analysis tasks.

##  VI Encoding, decoding, and causal analyses on random network models

We first consider encoding, decoding, and causal analyses on random network
models, such as Watts-Strogatz model (with small-world properties) [107],
Erdos-Renyi model (each pair of nodes are connected according to a probability
quantity) [108], and Barabási–Albert model (with scale-free properties) [109].
These random network models are important in statistical physics and
mathematics (e.g., for analyzing percolation on small-world networks [110],
Erdos-Renyi networks [111], and scale-free networks [112]). Meanwhile, they
are prototypes in analyzing social [113, 114], biological [115, 116, 117, 118,
119, 120], and chemical [121, 122] networks. Therefore, the encoding,
decoding, and causal analyses implemented on these models can be further
generalized to diverse real networks with corresponding topology properties.
The main motivation of our analyses on random network models is to validate
the proposed analytic metrics of network relations and suggest practical
solutions of potential limitations.

Figure 6:  Encoding, decoding, and causal analyses on random network models
during node deletion processes. (a) Watts-Strogatz (WS), Erdos-Renyi (ER), and
Barabási–Albert (BA) random networks and their adjacent matrices. (b)
Variables used for network approximation and encoding analysis are shown as
the functions of iterations. (c) Examples of Fisher information matrices in
the first iteration are illustrated (upper parallel). Fisher information
quantities, the traces of Fisher information matrices, are calculated in each
iteration (deep blue, purple, and light blue lines in bottom parallel).
Meanwhile, the mean Fisher information quantities, which are averaged across
different $\Theta$ conditions, are also shown (green lines in bottom
parallel). (d-e) Two causality metrics, Granger causality and transfer
entropy, are positively correlated with each other for
$\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ or
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$. (f) Causality
metrics are asymmetric, i.e., Granger causality and transfer entropy for
$\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ are not correlated
with those for $\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$.

###  VI.1 The node deletion process

In Fig. 6, we implement encoding, decoding, and causal analyses during the
node deletion process. In our experiment, a Watts-Strogatz model (each node
originally connects with $\alpha=15$ nodes, and edges are randomly re-wired
according a probability of $\beta=0.7$), an Erdos-Renyi model (each pair of
nodes are connected according to a probability of $\rho=0.15$), and a
Barabási–Albert model (a connected network of $\kappa=50$ nodes is set as the
initialization seed) are generated and initially contain $300$ nodes (see Fig.
6a for illustrations). Note that all the network parameters used in
initialization are set for convenience, and our analyses do not critically
relay on these parameters. The node deletion process consists of $200$
iterations. In each iteration, we randomly delete one node and all related
edges from these three networks. Encoding, decoding, and causal analyses are
implemented between $\mathcal{X}_{\phi}^{b}$, the networks in the $n$-th
iteration ($n\in\mathbb{Z}\cap\left[1,200\right]$) and
$\mathcal{X}_{\phi}^{a}$, the networks in their initialized forms. Certainly,
one may notice that the decoding analysis has not been explicitly defined by
the above settings. In real cases, the decoding analysis should be defined
according to research demands. In our research, we present an instance of the
decoding analysis based on a parameter vector
$\Theta=\left(\theta_{1},\ldots,\theta_{10}\right)$ controlled by
$\mathcal{X}_{\phi}^{b}$. Parameter vector $\Theta$ is designed to affect
$\mathcal{X}_{\phi}^{a}$ and make it become
$\mathcal{X}_{\phi}^{a}+\varepsilon$, where
$\varepsilon\sim\mathcal{N}\left(\mathbf{0},\Theta\right)$ denotes the effects
of $\Theta$ on $\mathcal{X}_{\phi}^{a}$. For convenience, we consider a case
where $\Theta$ is defined as the degree vector
$\left(\theta_{1}=\operatorname{deg}\left(v_{1}\right),\ldots,\theta_{10}=\operatorname{deg}\left(v_{10}\right)\right)$
of a set of nodes $\\{v_{1},\ldots,v_{10}\\}$ randomly selected from network
$\mathcal{G}_{b}$. By repeating random sampling, we can obtain a set of
observations
$\\{\Theta^{i}=\left(\theta_{1}^{i},\ldots,\theta_{10}^{i}\right)\\}$ of the
parameter vector, each of which corresponds to an effect on
$\mathcal{X}_{\phi}^{a}$ to create an observation
$\mathcal{X}_{\phi}^{a}+\varepsilon_{i}$. Based on these settings, the
decoding analysis can be implemented to measure the information of $\Theta$
contained in the probability distribution of
$\mathcal{X}_{\phi}^{a}+\varepsilon$.

| Direction | Variables | WS network | ER network | BA network | Statistical significance
---|---|---|---|---|---|---
Node deletion | $\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ |  $\mathcal{T}_{g}$ and $\mathcal{T}_{t}$ | $R=0.9608$ | $0.9614$ | $0.9731$ | $p<10^{-10}$
$\mathcal{T}_{G}$ and $\mathcal{T}_{T}$ | $R=0.9617$ | $0.9476$ | $0.96$ | $p<10^{-10}$
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$ |  $\mathcal{T}_{g}$ and $\mathcal{T}_{t}$ | $R=0.7960$ | $0.7756$ | $0.9269$ | $p<10^{-10}$
$\mathcal{T}_{G}$ and $\mathcal{T}_{T}$ | $R=0.4821$ | $0.5113$ | $0.7942$ | $p<10^{-10}$
Edge rewiring | $\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ |  $\mathcal{T}_{g}$ and $\mathcal{T}_{t}$ | $R=0.9672$ | $0.9355$ | $0.9581$ | $p<10^{-10}$
$\mathcal{T}_{G}$ and $\mathcal{T}_{T}$ | $R=0.9497$ | $0.7118$ | $0.8568$ | $p<10^{-10}$
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$ |  $\mathcal{T}_{g}$ and $\mathcal{T}_{t}$ | $R=0.9618$ | $0.6834$ | $0.9641$ | $p<10^{-10}$
$\mathcal{T}_{G}$ and $\mathcal{T}_{T}$ | $R=0.9145$ | $0.3809$ | $0.9618$ | $p<10^{-10}$
Node adding | $\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ |  $\mathcal{T}_{g}$ and $\mathcal{T}_{t}$ | $R=0.6996$ | $0.7362$ | $0.7116$ | $p<10^{-10}$
$\mathcal{T}_{G}$ and $\mathcal{T}_{T}$ | $R=0.3229$ | $0.3599$ | $0.3409$ | $p<10^{-10}$
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$ |  $\mathcal{T}_{g}$ and $\mathcal{T}_{t}$ | $R=0.8816$ | $0.8911$ | $0.8161$ | $p<10^{-10}$
$\mathcal{T}_{G}$ and $\mathcal{T}_{T}$ | $R=0.7893$ | $0.8105$ | $0.7447$ | $p<10^{-10}$
Table 1: Correlation $R$ between Granger causality and transfer entropy on
Watts-Strogatz (WS), Erdos-Renyi (ER), and Barabási–Albert (BA) random
networks.

As $n$ increases, more topology properties are changed due to node deletion.
Therefore, $\mathcal{X}_{\phi}^{a}$ and $\mathcal{X}_{\phi}^{b}$ are expected
to become increasingly different during node deletion. In Fig. 6b, the
enlarging difference between $\mathcal{X}_{\phi}^{a}$ and
$\mathcal{X}_{\phi}^{b}$ can be reflected by the decreasing mutual information
and the slightly increasing information divergence from
$\mathcal{X}_{\phi}^{a}$ to $\mathcal{X}_{\phi}^{b}$. However, we also observe
reductions in the information divergence from $\mathcal{X}_{\phi}^{b}$ to
$\mathcal{X}_{\phi}^{a}$. We hypothesize that the inconsistent variation
trends of
$\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)$ and
$\mathcal{D}\left(\mathcal{X}_{\phi}^{b}\|\mathcal{X}_{\phi}^{a}\right)$ may
arise from the asymmetry properties of information divergence (i.e.,
information divergence is a kind of pseudo-distance because
$\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)\neq\mathcal{D}\left(\mathcal{X}_{\phi}^{b}\|\mathcal{X}_{\phi}^{a}\right)$).
Meanwhile, the information divergence generalized by the Laplacian-energy-
based approach in Sec. V may fail to reflect the actual divergence between
networks with different sizes. In real situations, we suggest one to tackle
this problem by analyzing the mean value of
$\mathcal{D}\left(\mathcal{X}_{\phi}^{a}\|\mathcal{X}_{\phi}^{b}\right)$ and
$\mathcal{D}\left(\mathcal{X}_{\phi}^{b}\|\mathcal{X}_{\phi}^{a}\right)$. In
Fig. 6c, we illustrate several instances of Fisher information matrices in our
decoding analysis. Using the Fisher information matrices generated in each
iteration, we can further calculate Fisher information quantities, the traces
of Fisher information matrices, to reflect the decoding precision of
$\mathcal{X}_{\phi}^{b}$ from $\mathcal{X}_{\phi}^{a}$. In our experiments,
the observations of decoding precision and their mean values persistently
oscillate across different iterations. This phenomenon may arise from the
randomness during the sampling of parameter vector $\Theta$ on network
$\mathcal{G}_{b}$. In real situations, a more monotonous trend of Fisher
information quantities may be observed if $\Theta$ is controlled by specific
global properties of network $\mathcal{G}_{b}$. In Fig. 6d-e, we demonstrate
that our proposed causality metrics, Granger causality and transfer entropy,
are positively correlated with each other. In the direction of
$\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$, we calculate $10$
pairs of
$\mathcal{T}_{g}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
and
$\mathcal{T}_{t}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
in each iteration, whose mean values are the corresponding Granger causality
$\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
and transfer entropy
$\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$.
The positive correlations between
$\mathcal{T}_{g}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
and
$\mathcal{T}_{t}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
(they are analyzed in pairs) on three random network models and the positive
correlations
$\mathcal{T}_{G}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
and
$\mathcal{T}_{T}\left(\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}\right)$
on three models are shown in Table 1 and Fig. 6d, which turn out to be
significantly strong. Similarly, we implement the same correlation analysis in
the direction of $\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$.
The positive correlations between corresponding variables on three models are
relatively strong as well (see Table 1 and Fig. 6e). Meanwhile, we show that
these exists no significant correlation between the causality measured in the
direction of $\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$ and
that measured in the direction of
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$ because causality is
generally asymmetric (see Fig. 6f). These results demonstrate the validity of
our proposed analytic network relation metrics.

Figure 7:  Encoding, decoding, and causal analyses on random network models
during rewiring processes. (a) Variables used in encoding analysis are shown
as the functions of iterations. Note that $\gamma$, the rationality of network
approximation, is not shown since it is constant during rewiring. (b) Examples
of Fisher information matrices in the first iteration are illustrated in the
upper parallel. Fisher information quantities (deep blue, purple, and light
blue lines) and their mean values averaged across different $\Theta$
conditions (black lines) are shown in the bottom parallel. (c-d) The positive
correlations between Granger causality and transfer entropy are shown.

###  VI.2 The rewiring process

One may reasonably suspect that the reduction of mutual information in Fig. 6
is caused by the decreasing dimension number of $\mathcal{X}_{\phi}^{b}$
rather than the enlarging difference between $\mathcal{X}_{\phi}^{a}$ and
$\mathcal{X}_{\phi}^{b}$. To disprove this possibility, we replace the node
adding process by a rewiring process. All the settings used for random network
initialization keep the same as those in Fig. 6. After the initialization, we
randomly select one node and rewire its edges in each iteration. The edge
rewiring rules can be set in diverse forms but should be different from the
edge wiring rules in original random networks (e.g., the rewiring rules in the
Barabási–Albert random network should not be preferential attachment).
Otherwise, the generated network after rewiring may not become increasingly
different from the original one as the iteration number $n$ increases. To
ensure the enlarging difference, we define the rewiring processes of
initialized Watts-Strogatz, Erdos-Renyi, and Barabási–Albert random networks
following the wring rules of Erdos-Renyi ($\rho=0.5$), Watts-Strogatz
($\alpha=40$ and $\beta=0.1$), and Erdos-Renyi ($\rho=0.5$) models,
respectively. Please see Sec. VI.1 for the precise meaning of parameters
$\rho$, $\alpha$, and $\beta$. Note that our rewiring rules are set only for
convenience. One can consider other settings in practice. During the rewiring
process, we follow the same idea described in Sec. VI.1 to define encoding,
decoding, and causal analyses between the original networks,
$\mathcal{X}_{\phi}^{a}$, and the generated networks after each time of
rewiring, $\mathcal{X}_{\phi}^{b}$.

In Fig. 7, the results of encoding, decoding, and causal analyses are shown.
During rewiring processes, we can see the decreasing mutual information and
the increasing information divergence (in both directions), which keeps
consistency with the enlarging difference between $\mathcal{X}_{\phi}^{a}$ and
$\mathcal{X}_{\phi}^{b}$ (see Fig. 7a). In Fig. 7b, we present the instances
of Fisher information matrices, based on which, we can further calculate
Fisher information quantities to measure the decoding precision of
$\mathcal{X}_{\phi}^{b}$ from $\mathcal{X}_{\phi}^{a}$. Similar to the
situation in Fig. 6c, the observations of decoding precision and their mean
quantities oscillate due to the random sampling effects on parameter vector
$\Theta$. In Fig. 7c-d and Table 1, we can see the strong positive
correlations between Granger causality and transfer entropy in the direction
of $\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ and in the
direction of $\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$. All
these results suggest that our observations in Fig. 6 are not coincidentally
caused by the decreasing dimension number of $\mathcal{X}_{\phi}^{b}$ since
the rewiring process does not change the network size.

###  VI.3 The node adding process

Certainly, one may still doubt whether the decreasing mutual information in
Figs. 6-7 is coincidentally caused by the decreasing Shannon entropy of
$\mathcal{X}_{\phi}^{b}$. To disprove this possibility, we replace the
rewiring process in Fig. 7 by the node adding process to re-conduct our
encoding, decoding, and causal analyses. Similar to the rewiring process, the
node adding process is designed to ensure the enlarging difference between the
original networks, $\mathcal{X}_{\phi}^{a}$, and the generated networks after
each time of node adding, $\mathcal{X}_{\phi}^{b}$ as well. Specifically, we
add a node in each iteration and connect it with existing nodes according to
certain wiring rules. The wiring rules are set to be distinct from those in
the original networks. For convenience, we design the wiring of added nodes in
initialized Watts-Strogatz, Erdos-Renyi, and Barabási–Albert networks
following the rules in Erdos-Renyi ($\rho=0.5$), Watts-Strogatz ($\alpha=40$
and $\beta=0.1$), and Erdos-Renyi ($\rho=0.5$) models, respectively.

In Fig. 8, we present the results of encoding, decoding, and causal analyses
during node adding processes. One can see that the decreasing mutual
information and the increasing information divergence (in both directions)
coexist with the increasing Shannon entropy of $\mathcal{X}_{\phi}^{b}$ (Fig.
8a). These results disprove the possible validity issues concerned above. The
decoding analysis results are shown in Fig. 8b while the causal analysis
results are illustrated in Fig. 8c-d and Table 1. Because these results are
presented in the same form as Figs. 6-7, here we no longer discuss them in
details.

Figure 8:  Encoding, decoding, and causal analyses on random network models
during node adding processes. (a) Variables used for network approximation and
encoding analysis are presented as the functions of iterations. (b) Examples
of Fisher information matrices in the first iteration (upper parallel), Fisher
information quantities (deep blue, purple, and light blue lines in the bottom
parallel), and mean Fisher information quantities averaged across different
$\Theta$ conditions (black lines in the bottom parallel) are illustrated.
(c-d) The positive correlations between Granger causality and transfer entropy
are suggested.

To this point, we have validated our encoding, decoding, and causal analyses
on representative random network models. Our proposed network relation metrics
are principally suggest as able to reflect the difference between complex
networks. Although some limitations of these metrics are observed in our
experiments (e.g., information divergence quantities measured in different
directions may have inconsistent variation trends), there exist potential
tricks to deal with them (e.g., average all the asymmetric network relation
metrics across the direction of
$\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ and the direction of
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$). Below, we apply
our theory on real complex networks with these tricks. We show how our theory
may promote the analyses of protein network (a node similarity analysis task)
and chemical compound networks (a network similarity analysis task).

##  VII Encoding, decoding, and causal analyses on real networks

Among various fundamental tasks in network analysis, assigning node similarity
in a given network and measuring network similarity in network ensembles are
two important ones, which are closely related to node classification and
network classification tasks in machine learning. Here we implemented these
two kinds of similarity measurement tasks under our theoretical framework,
respectively.

###  VII.1 The analysis of protein network

The data set used for node similarity measurement is a protein-protein
interaction network [123]. In the network, nodes (there are 2361 nodes of 13
types) correspond to proteins and edges (7182 edges in total) represent
protein-protein interactions in the yeast [123]. To apply our theory in node
similarity analysis, we need to first represent each node by a sub-network
consisting of this node, its neighbors, and all relevant edges. To ensure the
quality of such a representation, we primarily focus on the nodes with more
than 5 neighbors. Then, we measure similarity between these nodes by
implementing encoding, decoding, and causal analyses between their
corresponding sub-networks. The decoding analysis is implemented under the
same framework used in Sec. VI, where we set the parameter vector as
$\Theta=\left(\theta_{1},\ldots,\theta_{3}\right)$ (the vector length of
$\Theta$ is set as smaller than the minimum size of node representation to
support random sampling). To avoid the influence of the potential inconsistent
variation trends of asymmetric network relation metrics in different
directions (similarity should be symmetric), we average information
divergence, Fisher information, Granger causality, and transfer entropy across
the direction of $\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ and
the direction of $\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$ to
make them non-directional.

Figure 9:  Encoding, decoding, and causal analyses on a protein network.
Encoding, decoding, and causal analyses are only implemented among the nodes
with no less than five neighbors. Except for mutual information, the other
variables are averaged between two analysis directions (i.e., averaged between
$\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ and
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$) to make them non-
directional. (a) The network of 13 types of proteins is presented (the first
network in the upper parallel). Based on these non-directional encoding,
decoding, and causal analyses results, the constrained node centrality (e.g.,
the closeness) is measured on the nodes with more than 5 neighbors.
Specifically, the costs of edge traversal are assigned according to the
information divergence (KL) as well as the inverses of the mutual information
(MI), Fisher information quantity (FI), Granger causality (GC), and transfer
entropy (TE), to measure the constrained closeness on these nodes. As for the
nodes with no more than 5 neighbors, their centrality values are assigned as
the minimum value of constrained closeness. The derived constrained node
centrality distributions are shown, in which the colors of nodes are defined
by their closeness quantities. (b) The mean graph distance (GD), KL, MI, FI,
GC, and TE averaged within a protein class or across two different classes.
(c) The mean GD, KL, MI, FI, GC, and TE averaged within a protein class (W) or
across two protein classes (A) are compared in the violin plots, in which blue
lines indicate the variation of expectation values. (d) The variables used in
encoding, decoding, and causal analyses are shown as the functions of graph
distance. Figure 10:  Encoding, decoding, and causal analyses on chemical
compound networks. (a) Two instances of the chemical compounds selected from
different classes are presented, where the weighted adjacent matrices (upper
parallel) and the corresponding topology structures (bottom parallel) are
illustrated. Note that edges are colored according to chemical bond labels.
(b) The matrices of non-directional encoding, decoding, and causal analyses
results among all chemical compounds. Specifically, the encoding, decoding,
and causal relations among the first class of chemical compounds are
surrounded by white dotted boxes (marked as “Wi”). Those relations among the
second class of chemical compounds as well as those relations between the
first and the second classes of chemical compounds are outside of white dotted
boxes (marked as “Wo”). (c) The mean KL, MI, FI, GC, and TE averaged within
“Wi” or within “Wo” are compared in the violin plots, in which blue lines
indicate the variation of expectation values. Note that the mean KL is shown
in a logarithmic plot. (d) The $t$-SNE analyses are implemented based on the
non-directional encoding, decoding, causal, and Graph2Vec analyses results in
(b). (e-f) The assessment results of classification performance of a naive
Bayes classifier and a logistic classifier are shown as the functions of
training proportion, i.e., the proportion of training data among all the data,
to show the robustness of performance.

In Fig. 9a, we first illustrate the original protein-protein interaction
network. As one of the explored applications, we measure the node importance
in term of the closeness [124, 125] constrained by non-directional encoding,
decoding, and causality metrics (the constraint is realized by assigning the
costs of edge traversal according to these metrics). The spatial distributions
of constrained closeness are shown in Fig. 9a, which may be considered as the
reflection of node importance from the perspectives of encoding, decoding, and
causality in future explorations. Given our proposed network relation metrics,
our primary motivation is to verify if they can reflect the category relations
between different nodes, i.e., whether two nodes correspond to the same type
of protein. Because the neighbouring relations represent protein-protein
interactions in the yeast, the graph distance (the length of the shortest path
between two nodes) can naturally reflect the difference between proteins.
Therefore, we also include the graph distance in our analysis. In Fig. 9b, we
calculate the mean graph distance, encoding, decoding, and causality metrics
averaged within a protein class or across two different classes. Based on
these results, we can further compare the mean values of the graph distance
and our proposed network relation metrics within the same class or across
different classes. As Fig. 9c suggests, most of the analyzed variables can
correctly reflect the fact that the proteins from distinct classes have larger
difference (e.g., the average mutual information quantities between the
proteins from the same class are larger than those between the proteins from
different classes). The only one exception is the Fisher information quantity
calculated under the preset decoding framework because the defined parameter
vector $\Theta$ may not establish an informative decoding relation between
$\mathcal{X}_{\phi}^{a}$ and $\mathcal{X}_{\phi}^{b}$. In Fig. 9d, we
represent the averaged quantities of encoding, decoding, and causality
variables as the functions of graph distance. All variables are observed to
decrease as the graph distance increases. Except for information divergence,
the variation trends of other variables suggest their capacities to reflect
the decreasing similarity between the proteins with a larger graph distance.
As for information divergence, its averaged quantities decrease with the graph
distance because extremely large information divergence quantities may coexist
with small graph distance values within a class (e.g., class $5$ in Fig. 9b).
Although information divergence principally reflects the difference between
proteins (e.g., see Fig. 9c), it fails on concrete protein classes (e.g.,
class $5$).

In sum, we have suggested that our encoding, decoding, and causal analyses may
principally support the measurement of node similarity on a protein-protein
interaction network. Nevertheless, some proposed network relation metrics
(e.g., information divergence) may fail to correctly reflect the within-class
and cross-class differences of specific proteins, which remain as limitations
for future explorations.

###  VII.2 The analysis of chemical compound networks

The data used for network similarity measurement is a collection of mutagenic
aromatic and heteroaromatic nitro compounds (188 compounds in total) [126].
Each chemical compound is represented by a network of atoms and is classified
into two classes according to their mutagenic effects on specific gram
negative bacteria [126]. The first class includes all the chemical compounds
with mutagenic effects while the second class consists of the chemical
compounds without mutagenic effects. To reflect the chemical properties of
each compound $\mathcal{X}_{\phi}^{a}$, we define its weighted adjacent matrix
as

$\displaystyle W_{i,j}\left(a\right)$
$\displaystyle=\beta_{i}\beta_{j}+\xi_{i,j},$ (51)

where $\beta_{i}$ denotes the label of atom $i$ and $\xi_{i,j}$ denotes the
label of chemical bond between atoms $i$ and $j$. Based on this definition, we
conduct encoding, decoding, and causal analyses between chemical compounds to
measure their similarity. Specifically, the decoding analysis is implemented
following the same way described in Sec. VI, where the parameter vector is set
as $\Theta=\left(\theta_{1},\ldots,\theta_{15}\right)$.

In Fig. 10a, we show several instances from two chemical compound classes.
According to the definition of these two classes, it is expected that the
first class of chemical compounds have higher similarity with each other
because they all have similar mutagenic effects on gram negative bacteria. As
for the second class of chemical compounds, we naturally speculate that they
have lower within-class and cross-class similarities because there may exist
various factors to make them irrelevant to gram negative bacteria (the
diversity of chemical compounds in the second class is higher). Fig. 10b-c
present the results of encoding, decoding, and causal analyses to reflect the
similarity between chemical compounds (all asymmetric network relation metrics
are averaged across the direction of
$\mathcal{X}_{\phi}^{b}\rightarrow\mathcal{X}_{\phi}^{a}$ and the direction of
$\mathcal{X}_{\phi}^{a}\rightarrow\mathcal{X}_{\phi}^{b}$). We discover that
information divergence quantities calculated within the first class are
smaller than those calculated in other situations. Other network relation
variables are consistently larger when they are calculated within the first
class. These results generally validate the above speculations and suggest the
possibility to distinguish between these two chemical compound classes based
on our proposed theory. To demonstrate this possibility, we first implement
the $t$-SNE analysis [127], a kind of unsupervised machine learning approach,
on the matrices shown in Fig. 10b. We treat each raw of these matrices as a
feature vector of the corresponding chemical compound, based on which we can
embed the compound into a space with high data separability. For comparison,
we also provide the $t$-SNE results of the feature vectors generated by the
Graph2Vec framework, a graph representation learning approach (note that we
set the feature vector length as $32$) [31]. Fig. 10d presents our $t$-SNE
analysis results, where two classes of chemical compounds generally form two
separable clusters after embedding. Some embedding results based on our theory
have better separability than those obtained from the Graph2Vec. In Fig.
10e-f, we further conduct supervised machine learning tasks on the feature
vectors defined above to classify chemical compounds. Based on the results of
encoding, decoding, and causal analyses, even simple learners such as a naive
Bayes classifier (Fig. 10e) and a logistic classifier (Fig. 10f) can achieve
high and robust learning performance (measured by the F1 score) on the task
originally prepared for advanced embedding frameworks (e.g., the Graph2Vec
[31]) or learners (e.g., graph neural networks [128]).

In sum, we have realized the similarity analysis of complex chemical compounds
by our encoding, decoding, and causal analyses, which further supports machine
learning tasks on chemical compounds with ideal performance.

##  VIII Discussion

###  VIII.1 Summary of our work

In this paper, we have introduced a possible way to measure analytic relations
between complex networks. Our theory explores a mapping $\phi$ to map an
arbitrary network $\mathcal{G}\left(V,E\right)$ to a random variable
$\mathcal{X}_{\phi}$ distributed on node set $V$. The random variable is
defined as
$\mathcal{X}_{\phi}\sim\mathcal{N}\left(\mathbf{0},L+\frac{1}{n}J\right)$, a
Gaussian variable characterized by a function of the discrete Schrödinger
operator $L$ (or referred to as the combinatorial graph Laplacian) of network
$\mathcal{G}$. On the one hand, such a definition ensures that the average
smoothness of mapping $\phi$ on network $\mathcal{G}$ is fully determined by
the information of network topology properties contained in discrete
Schrödinger operator $L$. On the other hand, this definition satisfies the
requirements of maximum entropy property of variable $\mathcal{X}_{\phi}$ to
promote the applicability in measuring information quantities. In Sec. VIII.3,
we shall suggest a dual definition of
$\mathcal{X}_{\phi}\sim\mathcal{N}\left(\mathbf{0},L+\frac{1}{n}J\right)$ and
discuss their differences from the perspective of network representation.

Based on the proposed $\mathcal{X}_{\phi}$, we further define encoding
(information divergence and mutual information), decoding (Fisher
information), and causal analyses (Granger causality and transfer entropy)
between any pair of complex networks. We have validated these analyses on
random network models (node deletion, edge rewiring, and node adding
processes), protein-protein interaction network, and chemical compound network
ensembles, suggesting their abilities to analytically measure network
similarity and robustly promote network classification tasks.

An open-source release of our algorithms can be seen in Ref. [59].

###  VIII.2 Progress compared with previous works

Although mainstream network relation analysis approaches developed in computer
science, such as embedding- [30, 32], matching- [33, 34, 35], and kernel-based
[36, 37] methods, have been extensively applied in various science and
engineering fields (e.g., in neuroscience [49, 50, 37, 51]), these approaches
either suffer from the high computational complexity (e.g., matching-based
[33]) or the dependence on empirical choice of network features (e.g.,
embedding-based [30]) and kernel functions (e.g., kernel-based [48]). Most
importantly, these approaches do not present explicit and analytic expressions
of universal network relations. Consequently, they are not applicable to
developing mathematics and statistical physics theories of complex networks,
as theoretical derivations usually require analytic and explainable
definitions rather than numerical solutions.

Compared with previous works, the main progress accomplished in our research
is to suggest a general way to represent the topology properties of an
arbitrary network by a Gaussian variable characterized by the function of the
discrete Schrödinger operator $L$. The effectiveness of such a representation
has been demonstrated by our experiments on random network models and real
networks. We notice that similar ideas can be seen in several studies of graph
signal, where researchers apply factor analysis [78] and low-rank models [79,
80] to suggest that graph Laplacian learning can be promoted by a Gaussian
Markov random field representation defined whose covariance matrix as
$\Sigma=L^{\dagger}$ [76, 77]. Recently, such an idea has been applied in
solving graph optimal transport problems to support graph comparison and
classification [53, 129]. The difference between these previous studies [76,
77, 53, 129] and our work lies in that they pre-set the covariance matrix of
Gaussian variable as $\Sigma=L^{\dagger}$ while we make progress in
analytically deriving $\Sigma=L+\frac{1}{n}J$ according to smoothness and
maximum entropy properties (see Sec. VIII.3 for a mathematical analysis of the
difference). As suggested by Sec. III, our derived results make the covariance
matrix invertible, ensuring the possibility to calculate numerous quantities
defined with $\Sigma^{-1}$ in Sec. IV.

Based on the proposed representation, we can directly generalize the metrics
originally developed for random variable relations to analyze complex
networks. In this work, we primarily illustrate representative metrics from
the perspectives of encoding, decoding, and causal analyses. One can further
consider Fisher-Rao metric in information geometry [130, 131] and
Wasserstein-2 metric in optimal transport [130, 131, 53, 129] in future, both
of which can be readily generalized to Gaussian variables.

###  VIII.3 Mathematical relations between our theory and related results

To understand the difference between our work and previous studies [76, 77,
53, 129], we begin with discussing the meaning of the following covariance
matrix

$\displaystyle\Sigma^{\heartsuit}=L^{\dagger}+\frac{1}{n}J,$ (52)

which is directly related to the covariance matrix in Refs. [76, 77, 53, 129].
As suggested by Eq. (9), the covariance matrix in Eq. (52) is exactly the
inverse of our result $\Sigma=L+\frac{1}{n}J$, i.e,
$\Sigma^{\heartsuit}=\Sigma^{-1}$. If one defines a Gaussian variable
$\mathcal{X}_{\phi}^{\heartsuit}\sim\mathcal{N}\left(\mathbf{0},\Sigma^{\heartsuit}\right)$,
then its precision matrix $Q^{\heartsuit}$ equals our proposed covariance
matrix

$\displaystyle Q^{\heartsuit}:={\Sigma^{\heartsuit}}^{-1}=\Sigma.$ (53)

Because the partial correlation between $X_{\phi}^{\heartsuit}\left(i\right)$
and $X_{\phi}^{\heartsuit}\left(j\right)$, the actual values of
$\mathcal{X}_{\phi}^{\heartsuit}$ on nodes $v_{i}$ and $v_{j}$, is fully
characterized by the precision matrix [132]

$\displaystyle\operatorname{corr}\left(X_{\phi}^{\heartsuit}\left(i\right),X_{\phi}^{\heartsuit}\left(j\right)\big{|}\mathcal{X}_{\phi}^{\heartsuit}\setminus\\{X_{\phi}^{\heartsuit}\left(i\right),X_{\phi}^{\heartsuit}\left(j\right)\\}\right)$
$\displaystyle:=$
$\displaystyle-\frac{Q^{\heartsuit}_{ij}}{\sqrt{Q^{\heartsuit}_{ii}Q^{\heartsuit}_{jj}}}$
(54) $\displaystyle=$
$\displaystyle-\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}},$ (55)

it is trivial to know that $X_{\phi}^{\heartsuit}\left(i\right)$ and
$X_{\phi}^{\heartsuit}\left(j\right)$ are expected to have a stronger partial
correlation if nodes $v_{i}$ and $v_{j}$ are connected by an edge with larger
weight (i.e., a larger value of $-\Sigma_{ij}$). Therefore, the Gaussian
variable defined by $\Sigma^{\heartsuit}$ is more applicable to the cases
where edge weights represent node similarity (e.g., see functional brain
networks where edges are functional correlations between different cortices
[49]). According to Eq. (10). the expected smoothness of mapping $\phi$ in
such a Gaussian variable is

$\displaystyle\mathbb{E}\left(\mathcal{S}\left(\phi\right)\right)$
$\displaystyle=\mathbb{E}\left(\mathcal{X}_{\phi}\right)^{T}L\mathbb{E}\left(\mathcal{X}_{\phi}\right)+\operatorname{tr}\left(L\Sigma^{\heartsuit}\right),$
(56)
$\displaystyle=\operatorname{tr}\left[L\left(L^{\dagger}+\frac{1}{n}J\right)\right],$
(57)
$\displaystyle=\operatorname{tr}\left(I\right)-\frac{1}{n}\operatorname{tr}\left(J\right)+\operatorname{tr}\left[L\left(I-L^{\dagger}L\right)\right],$
(58) $\displaystyle=n-1.$ (59)

Eqs. (58-59) are derived from the fact that $LL^{\dagger}L=L$ [65] and
$LL^{\dagger}=L^{\dagger}L=I-\frac{1}{n}J$ [75, 64, 74]. In general, Eqs.
(56-59) mean that the expected smoothness of mapping $\phi$ in a network
characterized by
$\mathcal{X}_{\phi}^{\heartsuit}\sim\mathcal{N}\left(\mathbf{0},\Sigma^{\heartsuit}\right)$
is independent of the network topology properties conveyed by $L$. On the
contrary, the expected smoothness is fully determined by the network size. We
speculate that this property may limit the capacity of
$\mathcal{X}_{\phi}^{\heartsuit}\sim\mathcal{N}\left(\mathbf{0},\Sigma^{\heartsuit}\right)$
to describe complex networks with high heterogeneity.

In our results, the Gaussian variable is
$\mathcal{X}_{\phi}\sim\mathcal{N}\left(\mathbf{0},L+\frac{1}{n}J\right)$,
where the expected smoothness of mapping $\phi$ is fully characterized by $L$
via
$\mathbb{E}\left(\mathcal{S}\left(\phi\right)\right)=\operatorname{tr}\left(LL\right)$
(see Eq. (14)). The covariance matrix $\Sigma$ of such a random variable
implies that $X_{\phi}\left(i\right)$ and $X_{\phi}\left(j\right)$ are
expected to evolve inversely (i.e., stronger negative covariance) if nodes
$v_{i}$ and $v_{j}$ are connected by an edge with larger weight. In other
words, the Gaussian variable defined by $\Sigma$ is more applicable to the
cases where edge weights measure node distance or difference (e.g., see
traffic networks where edge weights are distances between locations).
Meanwhile, it may have higher potential to characterize heterogeneous networks
where node difference matters.

In sum, the proposed covariance matrix $\Sigma=L+\frac{1}{n}J$ and its inverse
$\Sigma^{\heartsuit}=L^{\dagger}+\frac{1}{n}J$ are applicable to opposite
conditions, respectively. In application, one should choose between them
according to research demands. Although we primarily use $\Sigma$ to define
encoding, decoding, and causal analyses in our paper, all derived results can
be readily reformulated using $\Sigma^{\heartsuit}$ (one only needs to replace
$L$ by $L^{\dagger}$ and replace $L^{\dagger}$ by $L$ in our formulas).
Meanwhile, numerous computational experiments of $\Sigma^{\heartsuit}$ have
been done by existing studies (e.g., see Ref. [53]). Here we no longer repeat
these derivations and experiments. Our released toolbox [59] allows users to
choose between $\Sigma$ and $\Sigma^{\heartsuit}$ for a better network
characterization. Because real networks inevitably become heterogeneous in
certain degrees, we suggest that $\Sigma$ may also support a promising network
representation when edge weights measure node similarity (e.g., see our
experiments on real networks in Sec. VII). In the future, one can further
conduct more comparisons between these two covariance matrices on empirical
data or explore an ideal combination of them.

###  VIII.4 Limitations

As an initial attempt, there certainly remain limitations in our work for
further exploration. Here we suggest two limitations whose solutions may
significantly advance related fields.

The first limitation arises from the requirements of non-negative edge weights
in defining the discrete Schrödinger operator $L$. There exist numerous real
networks whose edge weights can be negative (e.g., in neural populations,
excitatory and inhibitory synapses have positive and negative weights,
respectively). Although the effects of negative weights on the eigenvalues of
$L$ have drawn increasing attention (e.g., see Refs. [133, 134, 135]), an
optimal definition of $L$ on networks with negative weights remains rather
exclusive. Similarly, the second limitation occurs when one attempts to
consider networks with directed edges. While notable progress has been
accomplished in defining $L$ on directed networks [136, 137], these
achievements can not completely address our problems because an asymmetric
version of $L$ do not support the definition of Gaussian variable developed in
Sec. III (the covariance matrix $\Sigma$ must be symmetric). Overall, these
two critical limitations should be primarily considered in the future.

## Acknowledgments

Author Y.T. conceptualizes the idea, develops theoretical frameworks, designs
computational tools, and writes the manuscript. Authors H.D.H. and G.Z.X
contribute equally to literature collection, mathematics proofreading, and
manuscript revision. Author Z.Y.Z. contributes to manuscript proofreading and
technical support. Author P.S. contributes to idea conceptualization,
manuscript writing, and project supervision. This project is supported by the
Artificial and General Intelligence Research Program of Guo Qiang Research
Institute at Tsinghua University (2020GQG1017) as well as the Tsinghua
University Initiative Scientific Research Program.

## References

  * Boccaletti  _et al._ [2006] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Complex networks: Structure and dynamics, Physics reports 424, 175 (2006).
  * Biamonte  _et al._ [2019] J. Biamonte, M. Faccin, and M. De Domenico, Complex networks from classical to quantum, Communications Physics 2, 1 (2019).
  * Mülken and Blumen [2011] O. Mülken and A. Blumen, Continuous-time quantum walks: Models for coherent transport on complex networks, Physics Reports 502, 37 (2011).
  * Bianconi [2002] G. Bianconi, Quantum statistics in complex networks, Physical Review E 66, 056123 (2002).
  * Roudi and Hertz [2011] Y. Roudi and J. Hertz, Mean field theory for nonequilibrium network reconstruction, Physical review letters 106, 048702 (2011).
  * Dorogovtsev  _et al._ [2008] S. N. Dorogovtsev, A. V. Goltsev, and J. F. Mendes, Critical phenomena in complex networks, Reviews of Modern Physics 80, 1275 (2008).
  * Sánchez  _et al._ [2002] A. D. Sánchez, J. M. López, and M. A. Rodriguez, Nonequilibrium phase transitions in directed small-world networks, Physical review letters 88, 048701 (2002).
  * Zhou  _et al._ [2006] C. Zhou, L. Zemanová, G. Zamora, C. C. Hilgetag, and J. Kurths, Hierarchical organization unveiled by functional connectivity in complex brain networks, Physical review letters 97, 238103 (2006).
  * Bullmore and Sporns [2009] E. Bullmore and O. Sporns, Complex brain networks: graph theoretical analysis of structural and functional systems, Nature reviews neuroscience 10, 186 (2009).
  * Li  _et al._ [2011] C. Li, H. Wang, W. De Haan, C. Stam, and P. Van Mieghem, The correlation of metrics in complex networks with applications in functional brain networks, Journal of Statistical Mechanics: Theory and Experiment 2011, P11018 (2011).
  * Bassett and Sporns [2017] D. S. Bassett and O. Sporns, Network neuroscience, Nature neuroscience 20, 353 (2017).
  * Sporns  _et al._ [2000] O. Sporns, G. Tononi, and G. M. Edelman, Theoretical neuroanatomy: relating anatomical and functional connectivity in graphs and cortical connection matrices, Cerebral cortex 10, 127 (2000).
  * Jeong  _et al._ [2000] H. Jeong, B. Tombor, R. Albert, Z. N. Oltvai, and A.-L. Barabási, The large-scale organization of metabolic networks, Nature 407, 651 (2000).
  * Wagner and Fell [2001] A. Wagner and D. A. Fell, The small world inside large metabolic networks, Proceedings of the Royal Society of London. Series B: Biological Sciences 268, 1803 (2001).
  * Tanaka [2005] R. Tanaka, Scale-rich metabolic networks, Physical review letters 94, 168101 (2005).
  * Jeong  _et al._ [2001] H. Jeong, S. P. Mason, A.-L. Barabási, and Z. N. Oltvai, Lethality and centrality in protein networks, Nature 411, 41 (2001).
  * Yook  _et al._ [2004] S.-H. Yook, Z. N. Oltvai, and A.-L. Barabási, Functional and topological characterization of protein interaction networks, Proteomics 4, 928 (2004).
  * Wagner [2001] A. Wagner, The yeast protein interaction network evolves rapidly and contains few redundant duplicate genes, Molecular biology and evolution 18, 1283 (2001).
  * Vázquez  _et al._ [2002] A. Vázquez, R. Pastor-Satorras, and A. Vespignani, Large-scale topological and dynamical properties of the internet, Physical Review E 65, 066130 (2002).
  * Pastor-Satorras and Vespignani [2004] R. Pastor-Satorras and A. Vespignani,  _Evolution and structure of the Internet: A statistical physics approach_ (Cambridge University Press, 2004).
  * Kahng  _et al._ [2002] B. Kahng, Y. Park, and H. Jeong, Robustness of the in-degree exponent for the world-wide web, Physical Review E 66, 046107 (2002).
  * De Menezes and Barabási [2004] M. A. De Menezes and A.-L. Barabási, Fluctuations in network dynamics, Physical review letters 92, 028701 (2004).
  * de Menezes and Barabási [2004] M. A. de Menezes and A.-L. Barabási, Separating internal and external dynamics of complex systems, Physical review letters 93, 068701 (2004).
  * Newman [2001a] M. E. Newman, The structure of scientific collaboration networks, Proceedings of the national academy of sciences 98, 404 (2001a).
  * Tsallis and de Albuquerque [2000] C. Tsallis and M. P. de Albuquerque, Are citations of scientific papers a case of nonextensivity?, The European Physical Journal B-Condensed Matter and Complex Systems 13, 777 (2000).
  * Menczer [2004] F. Menczer, Correlated topologies in citation networks and the web, The European Physical Journal B 38, 211 (2004).
  * Sznajd-Weron and Sznajd [2000] K. Sznajd-Weron and J. Sznajd, Opinion evolution in closed community, International Journal of Modern Physics C 11, 1157 (2000).
  * Pluchino  _et al._ [2005] A. Pluchino, V. Latora, and A. Rapisarda, Changing opinions in a changing world: A new perspective in sociophysics, International Journal of Modern Physics C 16, 515 (2005).
  * Stauffer and Meyer-Ortmanns [2004] D. Stauffer and H. Meyer-Ortmanns, Simulation of consensus model of deffuant et al. on a barabasi–albert network, International Journal of Modern Physics C 15, 241 (2004).
  * Soundarajan  _et al._ [2014] S. Soundarajan, T. Eliassi-Rad, and B. Gallagher, A guide to selecting a network similarity method, in  _Proceedings of the 2014 Siam international conference on data mining_ (SIAM, 2014) pp. 1037–1045.
  * Narayanan  _et al._ [2017] A. Narayanan, M. Chandramohan, R. Venkatesan, L. Chen, Y. Liu, and S. Jaiswal, graph2vec: Learning distributed representations of graphs, arXiv preprint arXiv:1707.05005 (2017).
  * Attar and Aliakbary [2017] N. Attar and S. Aliakbary, Classification of complex networks based on similarity of topological network features, Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 091102 (2017).
  * Conte  _et al._ [2004] D. Conte, P. Foggia, C. Sansone, and M. Vento, Thirty years of graph matching in pattern recognition, International journal of pattern recognition and artificial intelligence 18, 265 (2004).
  * Foggia  _et al._ [2014] P. Foggia, G. Percannella, and M. Vento, Graph matching and learning in pattern recognition in the last 10 years, International Journal of Pattern Recognition and Artificial Intelligence 28, 1450001 (2014).
  * Yan  _et al._ [2016] J. Yan, X.-C. Yin, W. Lin, C. Deng, H. Zha, and X. Yang, A short survey of recent advances in graph matching, in  _Proceedings of the 2016 ACM on International Conference on Multimedia Retrieval_ (2016) pp. 167–174.
  * Gao  _et al._ [2010] X. Gao, B. Xiao, D. Tao, and X. Li, A survey of graph edit distance, Pattern Analysis and applications 13, 113 (2010).
  * Borgwardt  _et al._ [2020] K. Borgwardt, E. Ghisu, F. Llinares-López, L. O’Bray, B. Rieck, _et al._ , Graph kernels: State-of-the-art and future challenges, Foundations and Trends® in Machine Learning 13, 531 (2020).
  * De Santo  _et al._ [2003] M. De Santo, P. Foggia, C. Sansone, and M. Vento, A large database of graphs and its use for benchmarking graph isomorphism algorithms, Pattern Recognition Letters 24, 1067 (2003).
  * Cour  _et al._ [2006] T. Cour, P. Srinivasan, and J. Shi, Balanced graph matching, Advances in neural information processing systems 19 (2006).
  * Jiang  _et al._ [2017] B. Jiang, J. Tang, C. Ding, Y. Gong, and B. Luo, Graph matching via multiplicative update algorithm, Advances in neural information processing systems 30 (2017).
  * Enqvist  _et al._ [2009] O. Enqvist, K. Josephson, and F. Kahl, Optimal correspondences from pairwise constraints, in  _2009 IEEE 12th international conference on computer vision_ (IEEE, 2009) pp. 1295–1302.
  * Leordeanu  _et al._ [2009] M. Leordeanu, M. Hebert, and R. Sukthankar, An integer projected fixed point method for graph matching and map inference (2009).
  * Zaslavskiy  _et al._ [2008] M. Zaslavskiy, F. Bach, and J.-P. Vert, A path following algorithm for the graph matching problem, IEEE Transactions on Pattern Analysis and Machine Intelligence 31, 2227 (2008).
  * Shervashidze  _et al._ [2009] N. Shervashidze, S. Vishwanathan, T. Petri, K. Mehlhorn, and K. Borgwardt, Efficient graphlet kernels for large graph comparison, in  _Artificial intelligence and statistics_ (PMLR, 2009) pp. 488–495.
  * Kashima  _et al._ [2004] H. Kashima, K. Tsuda, and A. Inokuchi, Kernels for graphs, in  _Kernel methods in computational biology_ (MIT Press, 2004) pp. 155–170.
  * Borgwardt and Kriegel [2005] K. M. Borgwardt and H.-P. Kriegel, Shortest-path kernels on graphs, in  _Fifth IEEE international conference on data mining (ICDM’05)_ (IEEE, 2005) pp. 8–pp.
  * Horváth  _et al._ [2004] T. Horváth, T. Gärtner, and S. Wrobel, Cyclic pattern kernels for predictive graph mining, in  _Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining_ (2004) pp. 158–167.
  * Yanardag and Vishwanathan [2015] P. Yanardag and S. Vishwanathan, Deep graph kernels, in  _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_ (2015) pp. 1365–1374.
  * Mheich  _et al._ [2020] A. Mheich, F. Wendling, and M. Hassan, Brain network similarity: methods and applications, Network Neuroscience 4, 507 (2020).
  * Tomlinson  _et al._ [2022] C. E. Tomlinson, P. J. Laurienti, R. G. Lyday, and S. L. Simpson, A regression framework for brain network distance metrics, Network Neuroscience 6, 49 (2022).
  * Abbas  _et al._ [2020] K. Abbas, E. Amico, D. O. Svaldi, U. Tipnis, D. A. Duong-Tran, M. Liu, M. Rajapandian, J. Harezlak, B. M. Ances, and J. Goñi, Geff: Graph embedding for functional fingerprinting, NeuroImage 221, 117181 (2020).
  * Cimini  _et al._ [2019] G. Cimini, T. Squartini, F. Saracco, D. Garlaschelli, A. Gabrielli, and G. Caldarelli, The statistical physics of real-world networks, Nature Reviews Physics 1, 58 (2019).
  * Petric Maretic  _et al._ [2019] H. Petric Maretic, M. El Gheche, G. Chierchia, and P. Frossard, Got: an optimal transport framework for graph comparison, Advances in Neural Information Processing Systems 32 (2019).
  * Shuman  _et al._ [2013] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains, IEEE signal processing magazine 30, 83 (2013).
  * Sandryhaila and Moura [2013] A. Sandryhaila and J. M. Moura, Discrete signal processing on graphs, IEEE transactions on signal processing 61, 1644 (2013).
  * Sandryhaila and Moura [2014] A. Sandryhaila and J. M. Moura, Discrete signal processing on graphs: Frequency analysis, IEEE Transactions on Signal Processing 62, 3042 (2014).
  * Cover [1999] T. M. Cover,  _Elements of information theory_ (John Wiley & Sons, 1999).
  * Hlaváčková-Schindler  _et al._ [2007] K. Hlaváčková-Schindler, M. Paluš, M. Vejmelka, and J. Bhattacharya, Causality detection based on information-theoretic approaches in time series analysis, Physics Reports 441, 1 (2007).
  * Tian  _et al._ [2022] Y. Tian, H. Hou, G. Xu, Y. Wang, Z. Zhang, and P. Sun, A toolbox for analytic relations between complex networks: encoding, decoding, and causality (2022), open source codes available at https://github.com/doloMing/Analytic-relations-between-complex-networks-encoding-decoding-and-causality.
  * Biyikoglu  _et al._ [2007] T. Biyikoglu, J. Leydold, and P. F. Stadler,  _Laplacian eigenvectors of graphs: Perron-Frobenius and Faber-Krahn type theorems_ (Springer, 2007).
  * Chung and Graham [1997] F. R. Chung and F. C. Graham,  _Spectral graph theory_ , 92 (American Mathematical Soc., 1997).
  * Horn and Johnson [2012] R. A. Horn and C. R. Johnson,  _Matrix analysis_ (Cambridge university press, 2012).
  * Xiao and Gutman [2003] W. Xiao and I. Gutman, Resistance distance and laplacian spectrum, Theoretical chemistry accounts 110, 284 (2003).
  * Chebotarev and Shamis [2006] P. Chebotarev and E. Shamis, On proximity measures for graph vertices, arXiv preprint math/0602073 (2006).
  * Barata and Hussein [2012] J. C. A. Barata and M. S. Hussein, The moore–penrose pseudoinverse: A tutorial review of the theory, Brazilian Journal of Physics 42, 146 (2012).
  * Herbster  _et al._ [2005] M. Herbster, M. Pontil, and L. Wainer, Online learning over graphs, in  _Proceedings of the 22nd international conference on Machine learning_ (2005) pp. 305–312.
  * Wahba [1990] G. Wahba,  _Spline models for observational data_ (SIAM, 1990).
  * Zhang  _et al._ [2009] H. Zhang, Y. Xu, and J. Zhang, Reproducing kernel banach spaces for machine learning., Journal of Machine Learning Research 10 (2009).
  * Fukumizu  _et al._ [2004] K. Fukumizu, F. R. Bach, and M. I. Jordan, Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces, Journal of Machine Learning Research 5, 73 (2004).
  * Zhou [2003] D.-X. Zhou, Capacity of reproducing kernel spaces in learning theory, IEEE Transactions on Information Theory 49, 1743 (2003).
  * Chen  _et al._ [2014] Z. Chen, K. Zhang, L. Chan, and B. Schölkopf, Causal discovery via reproducing kernel hilbert space embeddings, Neural computation 26, 1484 (2014).
  * Brodu and Crutchfield [2022] N. Brodu and J. P. Crutchfield, Discovering causal structure with reproducing-kernel hilbert space $\varepsilon$-machines, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 023103 (2022).
  * Summers  _et al._ [2015] T. Summers, I. Shames, J. Lygeros, and F. Dörfler, Topology design for optimal network coherence, in  _2015 European Control Conference (ECC)_ (IEEE, 2015) pp. 575–580.
  * Van Mieghem  _et al._ [2017] P. Van Mieghem, K. Devriendt, and H. Cetinay, Pseudoinverse of the laplacian and best spreader node in a network, Physical Review E 96, 032311 (2017).
  * Gutman and Xiao [2004] I. Gutman and W. Xiao, Generalized inverse of the laplacian matrix and some applications, Bulletin (Académie serbe des sciences et des arts. Classe des sciences mathématiques et naturelles. Sciences mathématiques) , 15 (2004).
  * Dong  _et al._ [2016] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, Learning laplacian matrix in smooth graph signal representations, IEEE Transactions on Signal Processing 64, 6160 (2016).
  * Kalofolias [2016] V. Kalofolias, How to learn a graph from smooth signals, in  _Artificial Intelligence and Statistics_ (PMLR, 2016) pp. 920–929.
  * Bartholomew  _et al._ [2011] D. J. Bartholomew, M. Knott, and I. Moustaki,  _Latent variable models and factor analysis: A unified approach_ , Vol. 904 (John Wiley & Sons, 2011).
  * Tipping and Bishop [1999] M. E. Tipping and C. M. Bishop, Probabilistic principal component analysis, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 61, 611 (1999).
  * Roweis [1997] S. Roweis, Em algorithms for pca and spca, Advances in neural information processing systems 10 (1997).
  * Olver and Townsend [2013] S. Olver and A. Townsend, Fast inverse transform sampling in one and two dimensions, arXiv preprint arXiv:1307.1223 (2013).
  * Kozachenko and Leonenko [1987] L. F. Kozachenko and N. N. Leonenko, Sample estimate of the entropy of a random vector, Problemy Peredachi Informatsii 23, 9 (1987).
  * Kraskov  _et al._ [2004] A. Kraskov, H. Stögbauer, and P. Grassberger, Estimating mutual information, Physical review E 69, 066138 (2004).
  * Mardia and Marshall [1984] K. V. Mardia and R. J. Marshall, Maximum likelihood estimation of models for residual covariance in spatial regression, Biometrika 71, 135 (1984).
  * Malagò and Pistone [2015] L. Malagò and G. Pistone, Information geometry of the gaussian distribution in view of stochastic optimization, in  _Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII_ (2015) pp. 150–162.
  * Casteigts  _et al._ [2012] A. Casteigts, P. Flocchini, W. Quattrociocchi, and N. Santoro, Time-varying graphs and dynamic networks, International Journal of Parallel, Emergent and Distributed Systems 27, 387 (2012).
  * Zimmermann  _et al._ [2004] M. G. Zimmermann, V. M. Eguíluz, and M. San Miguel, Coevolution of dynamical states and interactions in dynamic networks, Physical Review E 69, 065102 (2004).
  * Hill and Braha [2010] S. A. Hill and D. Braha, Dynamic model of time-dependent complex networks, Physical Review E 82, 046105 (2010).
  * Schreiber [2000] T. Schreiber, Measuring information transfer, Physical review letters 85, 461 (2000).
  * Staniek and Lehnertz [2008] M. Staniek and K. Lehnertz, Symbolic transfer entropy, Physical review letters 100, 158101 (2008).
  * Tian  _et al._ [2021] Y. Tian, Y. Wang, Z. Zhang, and P. Sun, Fourier-domain transfer entropy spectrum, Physical Review Research 3, L042040 (2021).
  * Shojaie and Fox [2022] A. Shojaie and E. B. Fox, Granger causality: A review and recent advances, Annual Review of Statistics and Its Application 9, 289 (2022).
  * Friston  _et al._ [2014] K. J. Friston, A. M. Bastos, A. Oswal, B. van Wijk, C. Richter, and V. Litvak, Granger causality revisited, Neuroimage 101, 796 (2014).
  * Bueso  _et al._ [2020] D. Bueso, M. Piles, and G. Camps-Valls, Explicit granger causality in kernel hilbert spaces, Physical Review E 102, 062201 (2020).
  * Marinazzo  _et al._ [2008] D. Marinazzo, M. Pellicoro, and S. Stramaglia, Kernel-granger causality and the analysis of dynamical networks, Physical review E 77, 056215 (2008).
  * Barnett  _et al._ [2009] L. Barnett, A. B. Barrett, and A. K. Seth, Granger causality and transfer entropy are equivalent for gaussian variables, Physical review letters 103, 238701 (2009).
  * Davidson  _et al._ [2004] R. Davidson, J. G. MacKinnon, _et al._ ,  _Econometric theory and methods_ , Vol. 5 (Oxford University Press New York, 2004).
  * Lizier and Prokopenko [2010] J. T. Lizier and M. Prokopenko, Differentiating information transfer and causal effect, The European Physical Journal B 73, 605 (2010).
  * Barnett and Bossomaier [2012] L. Barnett and T. Bossomaier, Transfer entropy as a log-likelihood ratio, Physical review letters 109, 138105 (2012).
  * Ay and Polani [2008] N. Ay and D. Polani, Information flows in causal networks, Advances in complex systems 11, 17 (2008).
  * Hlavácková-Schindler [2011] K. Hlavácková-Schindler, Equivalence of granger causality and transfer entropy: A generalization, Applied Mathematical Sciences 5, 3637 (2011).
  * Lizier  _et al._ [2008] J. T. Lizier, M. Prokopenko, and A. Y. Zomaya, Local information transfer as a spatiotemporal filter for complex systems, Physical Review E 77, 026110 (2008).
  * Gutman and Zhou [2006] I. Gutman and B. Zhou, Laplacian energy of a graph, Linear Algebra and its applications 414, 29 (2006).
  * Baruah and Bharali [2017] D. Baruah and A. Bharali, A comparative study of vertex deleted centrality measures, Ann Pure Appl Math 14, 199 (2017).
  * Yang  _et al._ [2017] X.-H. Yang, Q.-P. Zhu, Y.-J. Huang, J. Xiao, L. Wang, and F.-C. Tong, Parameter-free laplacian centrality peaks clustering, Pattern Recognition Letters 100, 167 (2017).
  * Qi  _et al._ [2012] X. Qi, E. Fuller, Q. Wu, Y. Wu, and C.-Q. Zhang, Laplacian centrality: A new centrality measure for weighted networks, Information Sciences 194, 240 (2012).
  * Watts and Strogatz [1998] D. J. Watts and S. H. Strogatz, Collective dynamics of ‘small-world’networks, nature 393, 440 (1998).
  * Erdős  _et al._ [1960] P. Erdős, A. Rényi, _et al._ , On the evolution of random graphs, Publ. Math. Inst. Hung. Acad. Sci 5, 17 (1960).
  * Barabási and Albert [1999] A.-L. Barabási and R. Albert, Emergence of scaling in random networks, science 286, 509 (1999).
  * Moore and Newman [2000] C. Moore and M. E. Newman, Epidemics and percolation in small-world networks, Physical Review E 61, 5678 (2000).
  * Almeira  _et al._ [2020] N. Almeira, O. V. Billoni, and J. I. Perotti, Scaling of percolation transitions on erdös-rényi networks under centrality-based attacks, Physical Review E 101, 012306 (2020).
  * Xulvi-Brunet  _et al._ [2003] R. Xulvi-Brunet, W. Pietsch, and I. Sokolov, Correlations in scale-free networks: Tomography and percolation, Physical review E 68, 036119 (2003).
  * Csányi and Szendrői [2004] G. Csányi and B. Szendrői, Structure of a large social network, Physical Review E 69, 036131 (2004).
  * Zekri and Clerc [2001] N. Zekri and J. P. Clerc, Statistical and dynamical study of disease propagation in a small world network, Physical Review E 64, 056115 (2001).
  * Hallquist and Hillary [2018] M. N. Hallquist and F. G. Hillary, Graph theory approaches to functional network organization in brain disorders: A critique for a brave new small-world, Network neuroscience 3, 1 (2018).
  * Zhao  _et al._ [2021] K. Zhao, Q. Zheng, T. Che, M. Dyrba, Q. Li, Y. Ding, Y. Zheng, Y. Liu, and S. Li, Regional radiomics similarity networks (r2sns) in the human brain: Reproducibility, small-world properties and a biological basis, Network Neuroscience 5, 783 (2021).
  * Benigni  _et al._ [2021] B. Benigni, A. Ghavasieh, A. Corso, V. d’Andrea, and M. De Domenico, Persistence of information flow: A multiscale characterization of human brain, Network Neuroscience 5, 831 (2021).
  * Sotero  _et al._ [2020] R. C. Sotero, L. M. Sanchez-Rodriguez, N. Moradi, and M. Dousty, Estimation of global and local complexities of brain networks: A random walks approach, Network Neuroscience 4, 575 (2020).
  * Shirai  _et al._ [2020] S. Shirai, S. K. Acharya, S. K. Bose, J. B. Mallinson, E. Galli, M. D. Pike, M. D. Arnold, and S. A. Brown, Long-range temporal correlations in scale-free neuromorphic networks, Network Neuroscience 4, 432 (2020).
  * Tian and Sun [2022] Y. Tian and P. Sun, Percolation may explain efficiency, robustness, and economy of the brain, Network Neuroscience , 1 (2022).
  * Martınez-Martınez  _et al._ [2021] C. Martınez-Martınez, J. Mendez-Bermudez, J. M. Rodrıguez, and J. M. Sigarreta, Computational and analytical studies of the harmonic index in erdös–rényi models, MATCH Commun. Math. Comput. Chem 85, 395 (2021).
  * Emmerich  _et al._ [2012] T. Emmerich, A. Bunde, and S. Havlin, Diffusion, annihilation, and chemical reactions in complex networks with spatial constraints, Physical Review E 86, 046103 (2012).
  * Batagelj and Mrvar [2006] V. Batagelj and A. Mrvar, Pajek datasets (2006).
  * Du  _et al._ [2015] Y. Du, C. Gao, X. Chen, Y. Hu, R. Sadiq, and Y. Deng, A new closeness centrality measure via effective distance in complex networks, Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 033112 (2015).
  * Newman [2001b] M. E. Newman, Scientific collaboration networks. ii. shortest paths, weighted networks, and centrality, Physical review E 64, 016132 (2001b).
  * Debnath  _et al._ [1991] A. K. Debnath, R. L. Lopez de Compadre, G. Debnath, A. J. Shusterman, and C. Hansch, Structure-activity relationship of mutagenic aromatic and heteroaromatic nitro compounds. correlation with molecular orbital energies and hydrophobicity, Journal of medicinal chemistry 34, 786 (1991).
  * Van der Maaten and Hinton [2008] L. Van der Maaten and G. Hinton, Visualizing data using t-sne., Journal of machine learning research 9 (2008).
  * Xu  _et al._ [2018] K. Xu, W. Hu, J. Leskovec, and S. Jegelka, How powerful are graph neural networks?, arXiv preprint arXiv:1810.00826 (2018).
  * Dong and Sawin [2020] Y. Dong and W. Sawin, Copt: Coordinated optimal transport on graphs, Advances in Neural Information Processing Systems 33, 19327 (2020).
  * Ning [2018] L. Ning, Smooth interpolation of covariance matrices and brain network estimation, IEEE transactions on automatic control 64, 3184 (2018).
  * Ning [2019] L. Ning, Smooth interpolation of covariance matrices and brain network estimation: Part ii, IEEE transactions on automatic control 65, 1901 (2019).
  * Rue and Held [2005] H. Rue and L. Held,  _Gaussian Markov random fields: theory and applications_ (Chapman and Hall/CRC, 2005).
  * Ahmadizadeh  _et al._ [2017] S. Ahmadizadeh, I. Shames, S. Martin, and D. Nešić, On eigenvalues of laplacian matrix for a class of directed signed graphs, Linear Algebra and its Applications 523, 281 (2017).
  * Bronski  _et al._ [2016] J. C. Bronski, L. DeVille, and T. Ferguson, Graph homology and stability of coupled oscillator networks, SIAM Journal on Applied Mathematics 76, 1126 (2016).
  * Bronski  _et al._ [2015] J. Bronski, L. DeVille, and K. P. Koutsaki, The spectral index of signed laplacians and their structural stability, arXiv preprint arXiv:1503.01069 (2015).
  * Bauer [2012] F. Bauer, Normalized graph laplacians for directed graphs, Linear Algebra and its Applications 436, 4193 (2012).
  * Agaev and Chebotarev [2005] R. Agaev and P. Chebotarev, On the spectra of nonsymmetric laplacian matrices, Linear Algebra and its Applications 399, 157 (2005).