# Learning low-rank latent mesoscale structures in networks

Hanbaek Lyu† $\dagger$ Department of Mathematics, University of Wisconsin-
Madison, WI 53706, USA  ,  Yacoub H. Kureh∗ ,  Joshua Vendrow∗ and  Mason A.
Porter∗,∗∗ $*$ Department of Mathematics, University of California, Los
Angeles, CA 90095, USA  $**$ Santa Fe Institute, Sante FE, NM 87501, USA
<EMAIL_ADDRESS>{ykureh, jvendrow<EMAIL_ADDRESS>
###### Abstract.

It is common to use networks to encode the architecture of interactions
between entities in complex systems in applications in the physical,
biological, social, and information sciences. To study the large-scale
behavior of complex systems, it is useful to study mesoscale structures in
networks as building blocks that influence such behavior [41, 54]. We present
a new approach for describing low-rank mesoscale structure in networks, and we
illustrate our approach using several synthetic network models and empirical
friendship, collaboration, and protein–protein interaction (PPI) networks. We
find that these networks possess a relatively small number of ‘latent motifs’
that together can successfully approximate most subgraphs of a network at a
fixed mesoscale. We use an algorithm that we call ‘network dictionary
learning’ (NDL) [64], which combines a network-sampling method [59] and
nonnegative matrix factorization [7, 64], to learn the latent motifs of a
given network. The ability to encode a network using a set of latent motifs
has a wide variety of applications to network-analysis tasks, such as
comparison, denoising, and edge inference. Additionally, using our new network
denoising and reconstruction (NDR) algorithm, we demonstrate how to denoise a
corrupted network by using only the latent motifs that one learns directly
from the corrupted networks.

Our code for the main algorithms and simulations is available at
https://github.com/HanbaekLyu/NDL_paper. At
https://github.com/jvendrow/Network-Dictionary-Learning, we provide a user-
friendly version as a Python package ndlearn.

It is often insightful to examine structures in networks [55] at intermediate
scales (i.e., at ‘mesoscales’) that lie between the microscale of nodes and
edges but below macroscale distributions of local network properties.
Researchers have considered subgraph patterns (i.e., the connection patterns
of subsets of nodes) as building blocks of network structure at various
mesoscales [69]. In particular, it is common to identify $k$-node (where $k$
is typically between $3$ and $5$) subgraph patterns in a network that are
unexpectedly common in comparison to some random-graph null model as ‘motifs’
of that network [9]. In the past two decades, the study of motifs has been an
important approach for analyzing various networked systems in many areas,
including biology [10, 31, 14, 47, 18], sociology [45, 24], and economics [30,
56]. However, to the best of our knowledge, researchers have not examined how
to use such motifs (or related mesoscale structures), after their discovery,
as building blocks to reconstruct a network. In this paper, we provide this
missing computational framework to bridge inferred subgraph-based mesoscale
structures and the global structure of networks. To do this, we propose (1) a
‘network dictionary learning’ (NDL) algorithm that learns ‘latent motifs’ from
samples of certain random $k$-node subgraphs and (2) a complementary algorithm
for ‘network denoising and reconstruction’ (NDR) that constructs a best
‘mesoscale linear approximation’ of a given network using the learned latent
motifs. We also provide a rigorous theoretical analysis of the proposed
algorithms. This analysis includes a novel result in which we prove that one
can accurately reconstruct an entire network if one has a dictionary of latent
motifs that can accurately approximate mesoscale structures of the network. We
compare the present work to related prior work [64] in the Methods section and
in our Supplementary Information (SI).

Using our approach, we find that various real-world networks (such as Facebook
friendship networks, Coronavirus and Homo sapiens protein–protein interaction
(PPI) networks, and an arXiv collaboration network) have low-rank subgraph
patterns, in the sense that one can successfully approximate their $k$-node
subgraph patterns by a weighted sum of a small number of latent motifs. The
latent motifs of these networks thereby reveal their low-rank mesoscale
structures. Our claim of the low-rank nature of a network’s mesoscale
structure concerns the space of certain subgraph patterns, rather than the
embedding of an entire network into a low-dimensional Euclidean space (as
considered in spectral-embedding and graph-embedding methods [46, 50]). It is
impossible to obtain such a low-dimensional graph embedding for networks with
small mean degrees and large clustering coefficients [65]. Additionally, as we
demonstrate in this paper, the ability to encode a network using a set of
latent motifs has a wide variety of applications in network analysis. Such
applications include network comparison, denoising, and edge inference.

###  $k$-path motif sampling

Computing all $k$-node subgraphs of a network is computationally expensive and
is the main computational bottleneck of traditional motif analysis [9]. Our
approach, which bypasses this issue, is to learn latent motifs by drawing
random samples of a particular class of $k$-node connected subgraphs. We
consider random $k$-node subgraphs that we obtain by uniformly randomly
sampling a ‘$k$-path’ from a network and including all edges between the
sampled nodes of the network. A $k$-path is a subgraph that consists of $k$
distinct nodes, with the $i$th node adjacent to the $(i+1)$th node for all
$i\in\\{1,\ldots,k-1\\}$. Sampling a $k$-path first serves two purposes: (1)
It ensures that the sampled $k$-node induced subgraph is connected with the
minimum number of imposed edges; and (2) it induces a natural node ordering of
the $k$-node induced subgraph. (Such an ordering is important for computations
that involve subgraphs.) By using the ‘$k$-walk’ motif-sampling algorithm in
[59] in conjunction with rejection sampling, we can sample a large number of
$k$-paths and obtain their associated induced subgraphs.

Figure 1.  Six subgraphs that are induced by uniformly sampled $k$-paths with
$k=20$ (red edges) from Caltech and UCLA Facebook networks, an Erdős–Rényi
(ER) random graph (ER1), a Barabási–Albert (BA) random graph (BA2), a
Watts–Strogatz (WS) small-world network (WS2), and a stochastic-block-model
(SBM) network (SBM1). See the Methods section for more details about these
networks.

The $k$-node subgraph that is induced by such a uniformly random $k$-path is
guaranteed to be connected and can exhibit diverse connection patterns (see
Figure 1), depending on the architecture of the original network. These
$k$-node subgraphs are the mesoscale structures that we consider in the
present paper. We study the connection patterns of such a random $k$-node
subgraph by decomposing it as a weighted sum of more elementary subgraph
patterns (possibly with continuous-valued edge weights), which we call latent
motifs. To study mesoscale structures in networks, we investigate several
questions. How many distinct latent motifs does one need to successfully
approximate all such $k$-node subgraph patterns? What do they look like? How
do these latent motifs differ for different networks?

### Dictionary learning and latent motifs

Dictionary-learning algorithms are machine-learning techniques that learn
interpretable latent structures of complex data sets. They are employed
regularly in the data analysis of text and images [15, 21, 26]. Dictionary-
learning algorithms usually consist of two steps. First, one samples a large
number of structured subsets of a data set (e.g., square patches of an image
or collections of a few sentences of a text); we refer to such a subset as a
mesoscale patch of a data set. Second, one finds a set of basis elements such
that taking a nonnegative linear combination of them can successfully
approximate each of the sampled mesoscale patches. Such a set of basis
elements is called a dictionary, and one can interpret each basis element as a
latent structure of the data set.

Figure 2.  Illustration of mesoscale structures that we learn from (a) images
and (b,c) networks. In each experiment in this figure, we form a matrix $X$ of
size $d\times n$ by sampling $n$ mesoscale patches of size $d=21\times 21$
from the corresponding object. For the image in panel a, the columns of $X$
are square patches of $21\times 21$ pixels. In panels b and c, we show both
networks and adjacency matrices. We take the columns of $X$ to be the $k\times
k$ adjacency matrices of the connected subgraphs that are induced by a path of
$k=21$ nodes, where a path of $k$ nodes consists of $k$ distinct nodes
$x_{1},\ldots,x_{k}$ such that $x_{i}$ and $x_{i+1}$ are adjacent for all
$i\in\\{1,\ldots,k-1\\}$. Using nonnegative matrix factorization (NMF), we
compute an approximate factorization $X\approx WH$ into nonnegative matrices
$W$ and $H$, where $W$ is called a dictionary and has $r=25$ columns. Because
of this factorization, we can approximate any sampled mesoscale patches (i.e.,
the columns of $X$) of an object by a nonnegative linear combination of the
columns of $W$, which we interpret as latent shapes for the images and latent
motifs (i.e., subgraphs) for the networks, respectively. The columns of $H$
give the coefficients in these linear combinations. The network dictionaries
of latent motifs that we learn from the (b) UCLA and (c) Caltech Facebook
networks have distinctive social structures. In the heat map of the UCLA
network, we show only the first 3000 nodes (according to the node labeling in
the data set). The image in panel a is from the collection Die Graphik Ernst
Ludwig Kirchners bis 1924, von Gustav Schiefler Band I bis 1916 (Accession
Number 2007.141.9, Ernst Ludwig Kirchner, 1926). We use the image with
permission from the National Gallery of Art in Washington, DC, USA.]

As an example, consider the artwork image in Figure 2a. We first uniformly
randomly sample 10,000 square patches of $21\times 21$ pixels and vectorize
them to get a $21^{2}\times\text{10,000}$ matrix $X$. The choice of
vectorization $\mathbb{R}^{k\times k}\rightarrow\mathbb{R}^{k^{2}}$ is
arbitrary; we use the column-wise vectorization in Algorithm A4 in the
Supplementary Information (SI). We then use a nonnegative matrix factorization
(NMF) [7] algorithm to find an approximate factorization $X\approx WH$, where
$W$ and $H$ are nonnegative matrices of sizes $21^{2}\times 25$ and
$25\times\text{10,000}$, respectively. Reshaping the columns of $W$ into
$21\times 21$ square images yields an image dictionary that describes ‘latent
shapes’ in the image.

Our network dictionary learning (NDL) algorithm for computing a ‘network
dictionary’ that consists of latent motifs is based on a similar idea. As
mesoscale patches of a network, we use the $k\times k$ binary (i.e.,
unweighted) matrices that encode connectivity patterns between the nodes that
form a uniformly random $k$-path. After obtaining sufficiently many mesoscale
patches of a network (e.g., by using a motif-sampling algorithm [59] with
rejection sampling), we apply a dictionary-learning algorithm (e.g., NMF [7])
to obtain latent motifs of the network.  A latent motif is a $k$-node weighted
network with nodes $\\{1,\ldots,k\\}$ and edges that have weights between $0$
and $1$. We use the term ‘on-chain edges’ for the edges of a latent motif
between nodes $i$ and $i+1$ for $i\in\\{1,\ldots,k-1\\}$; we use the term
‘off-chain edges’ for all other edges. We give more background about our NDL
algorithm in the Methods section and provide a complete implementation of our
approach in Algorithm NDL of the SI. We give theoretical guarantees for
Algorithm NDL in Theorems G.4 and G.7 in the SI.

In Figure 2, we compare 25 latent motifs with $k=21$ nodes of Facebook
friendship networks (which were collected on one day in fall 2005) from UCLA
(‘UCLA’) and Caltech (‘Caltech’) [34, 42]. Each node in one of these networks
is a Facebook account of an individual, and each edge encodes a Facebook
friendship between two individuals. The latent motifs reveal striking
differences between these networks in the connection patterns in the subgraphs
that are induced by $k$-paths with $k=21$. For example, the latent motifs in
UCLA’s dictionary (see Figure 2b) have sparse off-chain connections with a few
clusters, whereas Caltech’s dictionary (see Figure 2c) has relatively dense
off-chain connections. Most of Caltech’s latent motifs have ‘hub’ nodes (which
are adjacent to many other nodes in the latent motif) or communities with six
or more nodes. (See Figure 3 for community-size statistics.) An important
property of $k$-node latent motifs is that any network structure (e.g., hub
nodes, communities, and so on) in the latent motifs must also exist in actual
$k$-node subgraphs.  Indeed, we observe hubs and communities in the subgraph
samples from Caltech in Figure 1. By contrast, most of UCLA’s latent motifs do
not have such structures, as is also the case for the subgraph samples from
UCLA in Figure 1.

### Example networks

We demonstrate our approach using 16 example networks; 8 of them are real-
world networks and 8 of them synthetic networks. The 8 real-world networks are
Coronavirus PPI (for which we use the shorthand Coronavirus) [60, 66, 61] and
Homo sapiens PPI (for which we use the shorthand H. sapiens) [60, 50];
Facebook networks from Caltech, UCLA, Harvard, and MIT [34, 42]; SNAP Facebook
(for which we use the shorthand SNAP FB) [37, 50]; and arXiv ASTRO-PH (for
which we use the shorthand arXiv) [63, 50]. The first network is a
protein–protein interaction (PPI) network of proteins that are related to the
coronaviruses that cause Coronavirus disease 2019 (COVID-19), Severe Acute
Respiratory Syndrome (SARS), and Middle Eastern Respiratory Syndrome (MERS)
[66]. The second network is a PPI network of proteins that are related to
_Homo sapiens_ [60]. The third network is a 2012 Facebook network that was
collected from participants in a survey [37]. The fourth network is a
collaboration network based on coauthorship of preprints that were posted in
the astrophysics category of the arXiv preprint server. The last four real-
world networks are 2005 Facebook networks from four universities from the
Facebook100 data set [42]. In each Facebook network, nodes represent accounts
and edges encode Facebook ‘friendships’ between these accounts.

For the eight synthetic networks, we generate two instantiations each of
Erdős–Rényi (ER) $G(N,p)$ networks [1], Watts–Strogatz (WS) networks [5],
Barabási–Albert (BA) networks [6], and stochastic-block-model (SBM) networks
[3].  These are four types of well-studied random-graph models [55]. Each of
the ER networks has 5,000 nodes, and we independently connect each pair of
nodes with probabilities of $p=0.01$ (in the network that we call ER1) and
$p=0.02$ (in ER2). For the WS networks, we use rewiring probabilities of
$p=0.05$ (in WS1) and $p=0.1$ (in WS2) starting from a 5,000-node ring network
in which each node is adjacent to its $50$ nearest neighbors. For the BA
networks, we use $m=25$ (in BA1) and $m=50$ (in BA2), where $m$ denotes the
number of edges of each new node when it connects (via linear preferential
attachment) to the existing network, which we grow from an initial network of
$m$ isolated nodes (i.e., none of them are adjacent to any other node) until
it has 5,000 nodes.  The SBM networks SBM1 and SBM2 have three planted
1,000-node communities and two nodes in the $i_{0}$th and the $j_{0}$th
communities are connected by an edge independently with probability $0.5$ if
$i_{0}=j_{0}$ (i.e., if they are in the same community) and $0.001$ for SBM1
and $0.1$ for SBM2 if $i_{0}\neq j_{0}$ (i.e., if they are in different
communities).  See the Methods section for more details.

### Community sizes in subgraph samples and latent motifs

By comparing the $21$-node latent motifs of UCLA and Caltech (we extract
$r=25$ of each) in Figure 2, we observe that most latent motifs of Caltech
have larger communities than those of UCLA. To what extent does the community
structure of the latent motifs carries over to the subgraph samples in these
networks? Latent motifs are $k$-node networks with nonnegative edge weights,
so we can examine this question quantitatively by performing community
detection using a standard approach on these $k$-node networks. We will use
the locally greedy Louvain method for modularity maximization [22]. For most
of our example networks, we observe that the community-size statistics of the
learned latent motifs are close approximations of the corresponding statistics
for the subgraph samples of the networks. In Figure 3, we compare box plots of
community sizes of 10,000-node subgraphs that are induced by uniformly
randomly sampled $k$-paths to the corresponding box plots from community
detection of $r=25$ latent motifs of various networks. In our calculations,
the subgraphs and latent motifs from the networks UCLA and Caltech have the
same median community sizes. Moreover, the box plots for subgraphs and latent
motifs in Figure 3 are in good agreement for most of our example networks.

Figure 3.  A comparison of box plots of community sizes for 10,000 randomly
sampled subgraphs that are induced by uniformly random paths with $k=21$ nodes
(in red) to corresponding box plots from $r=25$ latent motifs of $k=21$ nodes
for various real-world and synthetic networks. We obtain communities of the
subgraphs and latent motifs by using the Louvain modularity-maximization
algorithm [22].

### Latent motifs of networks at various mesoscales

Because $k$-node latent motifs encode basic connectivity patterns of $k$ nodes
that are at most $k-1$ edges apart, one can interpret $k$ as a scale
parameter. With each latent motif, we associate a scalar ‘dominance score’,
which measures its total contribution in our reconstruction of the sampled
$k$-node subgraphs. (See Section D.2 of our SI for details.)  In Figure 4, we
show the two most-dominant latent motifs (i.e., the two with the largest
dominance scores)  that we learn from each of 15 example networks (i.e., for
all of them except Coronavirus PPI) at various scales (specifically, for
$k=6$, $k=11$, $k=21$, and $k=51$) when we use a dictionary with $r=25$ latent
motifs.

Figure 4.  Latent motifs that we learn from 16 networks (eight real-world
networks and eight synthetic networks, which include two distinct
instantiations of each of four random-graph models) at four different scales
(specifically, for $k=6$, $k=11$, $k=21$, and $k=51$), have distinct mesoscale
structures in the networks. Using NDL, we learn network dictionaries of $r=25$
latent motifs with $k$ nodes for each of the 16 networks. For each network at
each scale, we show the (top) first and (bottom) second most-dominant latent
motif from each dictionary. See Section D.2 of our SI for details about how we
measure latent-motif dominance.

As we increase the scale parameter $k$ for each network, various mesoscale
structures emerge in the latent motifs in Figure 4.  For instance, SNAP FB,
arXiv, and WS1 all have fully connected top (i.e., most dominant) latent
motifs at scale $k=6$ but their second latent motifs are distinct. In SNAP FB,
Caltech, and MIT at scales $k\in\\{6,11,21\\}$ and for the BA networks at
scales $k\in\\{11,21,51\\}$, the two most dominant latent motifs in Figure 4
have nodes that are adjacent to many other nodes in the latent motif. Hubs
(i.e., nodes that are adjacent to many other nodes) are characteristic of both
BA networks (which have heavy-tailed degree distributions) [6] and most social
networks (which typically have heavy-tailed degree distributions) [55]. We
also observe hubs in the network dictionaries of the latent motifs of the
Facebook networks UCLA and Harvard. (See Figure 14 in our SI.)

In Figure 3, we saw that community sizes in latent motifs reflect community
sizes of actual subgraphs in a network. Importantly, the type of community
structure that we examine is very different from typical network community
structure. For example, consider the WS networks. The top latent motif of the
network WS1 at scale $k=6$ is fully connected, but the top latent motif of WS2
is not fully connected because of its larger rewiring probability. At larger
scales (i.e., for larger $k$), both WS networks have latent motifs with
multiple communities. The WS networks have locally densely connected nodes on
a ring of nodes and random ‘shortcut’ edges that can connect distant nodes of
the ring. Therefore, when one samples a $k$-path uniformly at random, it is
very likely to use at least one shortcut edge. When a $k$-path uses a shortcut
edge, we expect the resulting sampled induced subgraph to have two distinct
densely connected communities. This local ‘community structure’ in the WS
networks is rather different from standard types of community structure [27,
49]. While we do observe such community structures on subgraphs that are
induced by $k$-paths as in Figure 1, it does not imply that the entire node
set of the WS networks is partitioned into a few communities.We can also see
the difference between our mesoscale structures and community structure by
examining the latent motifs of the SBM networks at different scales in Figure
4. The SBM networks have three (equal-sized) communities by construction, but
its latent motifs at none of the scales have three communities; uniformly
sampled $k$-paths do not always intersect with all three communities. For
example, the six $20$-paths from SBM1 in Figure 1 intersect with only one or
two of the network’s planted communities.

### Network denoising and reconstruction using latent motifs

Suppose that we are given a network $G$ and a dictionary $W$ of latent motifs
at scale $k$, where we may or may not learn $W$ from $G$. Consider the
following two scenarios. In one scenario, we suppose that we know $G$ exactly,
and we ask how to measure the ‘effectiveness’ of $W$ at approximating
mesoscale patches of $G$ at scale $k$. In the other scenario, we suppose that
$G$ is a noisy version of some true network $G_{\textup{true}}$ and that $W$
is ‘faithful’ in the sense that it can well-approximate mesoscale patches of
$G_{\textup{true}}$ at scale $k$. (See (7) in the SI for a precise formulation
of this idea.) We then ask how to infer the true network $G_{\textup{true}}$
from $G$ and $W$.

To examine the above questions, we develop an algorithm that we call network
denoising and reconstruction (NDR) (see Algorithm NDR) that takes a network
$G$ and network dictionary $W$ as input and outputs a weighted network
$G_{\textup{recons}}$ that has the same node set as $G$. Until it converges,
the NDR algorithm repeatedly samples mesoscale patches of $G$ at scale $k$ and
obtains a nonnegative linear approximation of those patches using the latent
motifs in $W$. Because each edge $e$ of $G$ can appear in multiple mesoscale
patches of $G$, this procedure can yield multiple reconstructed weights for
$e$. To obtain the final weight of $e$ in $G_{\textup{recons}}$, we calculate
their mean (see Algorithm NDR). The reconstructed network
$G_{\textup{recons}}$ is a version of the original network $G$ from the ‘lens’
of the latent motifs in the dictionary $W$. The latent motifs thereby dictate
the types of subgraph structures of $G$ that are captured in
$G_{\textup{recons}}$. This network-reconstruction process is inspired by
patch-based image reconstruction [15, 23]. We illustrate both patch-based
image reconstruction and our NDR algorithm in the Methods section.

### Network reconstruction using latent motifs

To measure the effectiveness of $W$ for an unweighted network $G$ (i.e., edges
are either present or absent and there are no multi-edges), one can threshold
the weighted edges of $G_{\textup{recons}}$ at some fixed value
$\theta\in[0,1]$ to obtain an undirected reconstructed network
$G_{\textup{recons}}(\theta)$ with binary edge weights (of either $0$ or $1$),
which one can then compare directly with the original unweighted network $G$.
We regard $W$ as effective at describing $G$ at scale $k$ if
$G_{\text{recons}}(\theta)$ is close to $G$ for some $\theta$. (We will
quantify our notion of ‘closeness’ in the next paragraph.) We interpret the
edge weights in $G_{\textup{recons}}$ as measures of confidence of the
corresponding edges in $G$ with respect to $W$. For example, if an edge $e$
has the smallest weight in $G_{\textup{recons}}$, we regard it as the most
‘outlying’ with respect to the latent motifs in $W$. See Theorems G.10 and
G.14 in the SI for theoretical guarantees and error bounds for NDR.

Figure 5.  Self-reconstruction and cross-reconstruction accuracy of several
real-world and synthetic networks versus the edge threshold $\theta$ and the
number $r$ of latent motifs in a network dictionary. The label $X\leftarrow Y$
indicates that we reconstruct network $X$ using a network dictionary that we
learn from network $Y$. The reconstruction process produces a weighted network
that we turn into an unweighted network by thresholding the edge weights at a
threshold value $\theta$; we keep only edges whose weights are strictly larger
than $\theta$. We measure reconstruction accuracy by calculating the Jaccard
index of an original network’s edge set and an associated reconstructed
network’s edge set. In panel a, we plot accuracies versus $\theta$ (with the
number of latent motifs fixed at $r=25$), where $X$ is one of five real-world
networks (two PPI networks, two Facebook networks, and one collaboration
network). In panels b–e, we reconstruct each of the four Facebook networks
using network dictionaries with $r\in\\{9,16,25,36,64,81,100\\}$ latent motifs
that we learn from one of eight networks (with the threshold value fixed at
$\theta=0.4$).

In Figure 5, we show the results of several network-reconstruction experiments
using several real-world networks and synthetic networks. We perform these
experiments for various values of the edge threshold $\theta\in[0,1]$ and
$r\in\\{9,16,25,36,49,81,100\\}$ latent motifs in a single dictionary. Each
network dictionary in Figure 5 has $k=21$ nodes, for which the dimension of
the space of all possible mesoscale patches (i.e., the adjacency matrices of
the induced subgraphs) is $\binom{21}{2}-20=190$. An important observation is
that one can reconstruct a given network using an arbitrary network
dictionary, which one can even learn from a different network. Such a ‘cross-
reconstruction’ allows us to quantitatively compare the learned mesoscale
structures of different networks. We label each subplot of Figure 5 with
$Y\leftarrow X$ to indicate that we are reconstructing network $Y$ using a
network dictionary that we learn from network $X$. (The reconstruction also
uses the original network $Y$.) We turn the weighted reconstructed networks
into unweighted networks by thresholding their edges using a threshold
$\theta$. We measure the reconstruction accuracy by calculating the Jaccard
index between the original network’s edge set and the reconstructed network’s
edge set. That is, to measure the similarity of two edge sets, we calculate
the number of edges in the intersection of these sets divided by the number of
edges in the union of these sets. This gives a ‘global’ measure of
reconstruction accuracy in the sense that if this measure is 1, then the
reconstructed network is precisely the same as the original one. We obtain the
same qualitative results as in Figure 5 if we instead measure similarity using
the Rand index [2]).

In Figure 5a, we plot the accuracy for ‘self-reconstruction’ $X\leftarrow X$
versus the threshold $\theta$ (with $r=25$), where $X$ is one of the real-
world networks Coronavirus, H. sapiens, SNAP FB, Caltech, and arXiv.  The
accuracies for Coronavirus, H. sapiens, and Caltech peak above $95\%$ when
$\theta\approx{0.4}$; the accuracies for arXiv and SNAP FB peak above $88\%$
and $70\%$, respectively, for $\theta\approx 0.6$. We choose $\theta=0.4$ for
the cross-reconstruction experiments for the Facebook networks Caltech,
Harvard, MIT, and UCLA in Figures 5b,c. These four Facebook networks have
self-reconstruction accuracies above $80\%$ for $r=25$ motifs with a threshold
of $\theta=0.4$. The total number of dimensions when using mesoscale patches
at scale $k=21$ is $190$, so this result suggests that all eight of these
real-world networks have low-rank mesoscale structures at scale $k=21$.

Figure 6.  A comparison of the degree distributions (which we show as
histograms) and the mean local clustering coefficients (which we write in the
legends as ‘mlcc’) of the original network and the reconstructed networks
using $r$ latent motifs at scale $k=21$ for the five networks in Figure 5a. In
panel a, we use the four unweighted reconstructed networks for Caltech with
$r\in\\{9,16,25,64\\}$ latent motifs that we used to compute the self-
reconstruction accuracy in 5b. As we increase $r$, we observe that the mean
local clustering coefficient increases towards its value in the original
network and that the degree distributions of the reconstructed networks
converge to that of the original network. (By increasing $r$, we are able to
include nodes with progressively larger degrees in the latent motifs.) In
panels b–e, we show the mean local clustering coefficients and degree
distributions of the unweighted reconstructed networks for Coronavirus, H.
sapiens, SNAP FB, and arXiv with $r=25$ latent motifs that we used to compute
the self-reconstruction accuracy in Figures 5b–e.

We gain further insights into our self-reconstruction experiments by comparing
the degree distributions and the mean local clustering coefficients of the
original and the unweighted reconstructed networks with threshold $\theta=0.4$
(see Figure 6). The mean local clustering coefficients of all reconstructed
networks are close to those of the corresponding original networks. In Figure
6a, we show that the degree distributions of the reconstructed networks for
Caltech with $r\in\\{9,16,25,64\\}$ latent motifs converge toward that of the
original network as we increase $r$. Reconstructing Caltech with larger values
of $r$ appears to gain accuracy by including nodes with a larger degree than
is possible for smaller values of $r$. In other words, ‘low-rank’
reconstructions (i.e., those with small values of $r$) of Caltech seem to
recover only a small number of the edges of each node, even though it is able
to achieve a large reconstruction accuracy (e.g., over $81\%$ for $r=9$).

We consider cross-reconstruction accuracies $Y\leftarrow X$ in Figures 5b,c,
where $Y$ is one of the Facebook networks Caltech, Harvard, MIT, and UCLA and
$X$ (with $X\neq Y$) is one of the these four networks or one of the four
synthetic networks ER2, WS2, BA2, and SBM2. From the cross-reconstruction
accuracies and examining the network structures of the the latent motifs (see
Section B.4 of the SI) in Figures 2 and 4 (also see Figures 14, 16, and 17 in
the SI), we draw the following conclusions at scale $k=21$. First, the
mesoscale structure of Caltech is distinct from those of Harvard, UCLA, and
MIT. This is consistent with prior studies of these networks (see, e.g., [42,
48]). Second, Caltech’s mesoscale structure at scale $k=21$ has a higher
dimension than those of the other three universities’ Facebook networks.
Third, Caltech has a lot more communities with at least $10$ nodes than the
other three universities’ Facebook networks (also see Figure 3). Fourth, the
BA network BA2 captures the mesoscale structures of MIT, Harvard, and UCLA at
scale $k=21$ better than the synthetic networks that we generate from the ER,
WS, and SBM models;  However, SBM2 captures the mesoscale structure of Caltech
better than all other networks in Figure 5b except for Caltech itself for
$r\in\\{9,15,25,49\\}$. See Section F.5 of the SI for further discussion.

Figure 7.  We prove that a Jaccard reconstruction error of weighted
reconstructed networks (i.e., without thresholding edge weights as in Figure
5) is upper-bounded by the mean approximation error of $k$-node subgraphs by
$k$-path latent motifs divided by $2(k-1)$. We state this result concisely as
the inequality in panel a. Given a network $G$ and a network dictionary $W$ of
$k$-path latent motifs, let $G_{\textup{recons};W}$ denote the weighted
reconstructed network that we obtain using our NDR algorithm (see Algorithm
NDR in our SI). We define the Jaccard distance in (82) in the SI. On the
right-hand side of the inequality, $\mathbf{x}$ is a uniformly random $k$-path
in $G$, the matrix $A_{\mathbf{x}}$ is the $k\times k$ adjacency matrix of the
subgraph that is induced by the node set of $\mathbf{x}$, and
$\hat{A}_{\mathbf{x};W}$ is the best nonnegative linear approximation of
$\mathbf{A}_{\mathbf{x}}$ that we obtain using the latent motifs in $W$. (See
Section G.3 and Theorem G.4 in our SI for precise statements of relevant
definitions and the mathematical result.) Using a crude network dictionary
that consists of a single $k$-path, this upper bound is the mean number of
off-chain edges in mesoscale patches of a network divided by $k-1$. For an ER
network with expected edge density $p$, this quantity equals $kp/2$.
Subtracting both sides of the inequality in panel a from $1$ yields a lower
bound for the Jaccard reconstruction accuracy. We plot this lower bound in
panels b–e for the same parameters as in Figures 5b–e (i.e., $k=21$ and
$r\in\\{6,16,25,36,49,64,81,100\\}$).

We also comment briefly about the cross-reconstruction experiments in Figure 5
that use latent motifs that we learn from ER networks. For instance, when
reconstructing MIT, Harvard, and UCLA using latent motifs that we learn from
ER2, we obtain a reconstruction accuracy of at least 72%. This may seem
unreasonable at first glance because the latent motifs that we learn from ER2
should not have any information about the Facebook networks. However, all of
these networks are sparse (with edge densities of at most $0.02$) and we are
sampling subgraphs using $k$-paths. The $k$-node subgraphs that are induced by
uniformly random $k$-paths in these sparse networks have only a few off-chain
edges (see Figure 1). For example, the $k$-node subgraphs that one samples
from the sparse ER network ER2 tend to have $k$-node paths and a few extra
off-chain edges (see Figure 1). A similarly sparse or sparser network, such as
UCLA (whose edge density is about $0.0036$) has similar subgraph patterns
(even though, unlike the subgraphs of ER2, the off-chain edges are not
independent). This is the reason that we can reconstruct some networks with
high accuracy by using latent motifs that we learn from a completely unrelated
network.

Our discussion in the previous paragraph suggests an important question. When
reconstructing a sparse network using a network dictionary, how much
reconstruction accuracy should we expect as a baseline? More generally, if we
have a network dictionary that is effective at approximating the mesoscale
patches of a network, how much reconstruction accuracy do we expect? In the
present paper, we state and prove a novel theorem that answers this question.
Specifically, we prove mathematically that a global Jaccard reconstruction
error of a weighted reconstructed network (i.e., without thresholding edge
weights as in Figure 5) is upper-bounded by the mean approximation error of
the mesoscale patches (i.e., $k$-node subgraphs) of a network by the $k$-path
latent motifs divided by $2(k-1)$. See Section G.3 and Theorem G.4 in our SI
for precise statements of the relevant definitions and the mathematical
result.

Our NDL algorithm (see Algorithm NDL) finds a network dictionary that
approximately minimizes the upper bound in the inequality in Figure 7a. Ssee
Theorem G.4 in the SI.) Using an arbitrary network dictionary is likely to
yield larger values of the upper bound. Therefore, according to the theorem in
Figure 7a, it is likely to yield a less accurate weighted reconstructed
network. For instance, consider using a network dictionary that consists of a
single $k$-path, the aforementioned upper bound is the mean number of off-
chain edges in mesoscale patches of a network divided by $k-1$. (See the right
inequality in Figure 7a.) For any ER network with expected edge density $p$,
this quantity equals $kp/2$ in expectation. At scale $k=20$ for ER2, this
value is $0.5$. Consequently, we expect a reconstruction accuracy of at least
$50\%$ when reconstructing ER2 using latent motifs (such as the ones for UCLA
in Figure 2) that have large on-chain entries and small off-chain entries.
Substituting the edge density of UCLA for $p$, we expect a reconstruction
accuracy of at least 94%. However, according to the first inequality in Figure
7a, the lower bound for the reconstruction accuracy that we obtain using
latent motifs of UCLA is about 80%. Therefore, for UCLA, we expect to obtain
many more off-chain edges in mesoscale patches than what we expect from an ER
network with the same edge density. We plot the lower bound of the
reconstruction accuracy in Figures 7b–e for the parameters in Figures 5b–e
(i.e., $k=21$ and $r\in\\{6,16,25,36,49,64,81,100\\}$). The lower bound for
the self-reconstructions is not too far from the actual reconstruction
accuracy for unweighted reconstructed networks in Figures 5b–e (it is within
$20\%$ for Caltech and UCLA and within $10\%$ for MIT and Harvard for all
$r$), but there is a much larger accuracy gap for the cross-reconstruction
experiments (e.g., there is at least a $50\%$ difference for
UCLA$\leftarrow$Caltech). This indicates that, even if one uses latent motifs
that are not very efficient at mesoscale approximation, one can obtain
unweighted global reconstructions that are significantly more accurate than
what is guaranteed by the theoretically proven bounds.

### Network denoising using latent motifs

In Figure 9, we use our NDL and NDR algorithms to perform two types of network
denoising. We can think of these two types of network denoising as distinct
binary-classification problems. One of them is network denoising with
subtractive noise (which is often called edge ‘prediction’ [70, 20, 33, 58,
62], and the other is network denoising with additive noise [57]. In each
scenario, we suppose that we are given an observed network $G=(V,E)$ with a
node set $V$ and an edge set $E$ and that we are asked to find an unknown
network $G^{\prime}=(V,E^{\prime})$ with the same node set $V$ but a possibly
different edge set $E^{\prime}$. We interpret $G$ as a corrupted version of a
‘true network’ $G^{\prime}$ that we observe with some uncertainty. We
interpret edges and non-edges in $G$ as ‘false relations’ and ‘false non-
relations’, respectively. In the subtractive-noise setting, we assume that $G$
is a partially observed version of $G^{\prime}$ (i.e., $E\subsetneq
E^{\prime}$), and we seek to classify all non-edges in $G$ as ‘positives’
(i.e., edges in $G^{\prime}$) or ‘negatives’ (i.e., non-edges in
$G^{\prime}$). In the additive-noise setting, we suppose that $G$ is a
corrupted version of $G^{\prime}$ that includes false edges (i.e., $E\supseteq
E^{\prime}$), and we seek to classify all edges in $G$ into ‘positives’ (i.e.,
edges in $G^{\prime}$) and ‘negatives’ (i.e., non-edges in $G^{\prime}$).

In our study of these two binary-classification problems, we use the following
five real-world networks: Caltech, SNAP FB, arXiv, Coronavirus, and H.
sapiens. Given a network $G=(V,E)$, our experiments proceed as follows. In the
subtractive-noise setting, we create a smaller network by removing half of the
existing edges, which we choose uniformly at random, such that the remaining
network is connected. We refer to this noise type as ‘$-\textup{ER}$’.  In the
additive-noise setting, we create two types of corrupted networks. We create
the first type of corrupted network by adding edges between node pairs that we
choose independently with a fixed probability so that $50\%$ of the edges in a
corrupted network are new. We refer to this noise type as ‘${+}\textup{ER}$’.
We create the second type of corrupted network by adding false edges in a
structured way by generating them using the WS model. We select 100 nodes for
four of the networks (the exception is that we select 500 for H. sapiens)
uniformly at random and generate 1,000 new edges (we generate 30,000 new edges
for H. sapiens) according to the WS model. In this corrupting WS network, each
node in a ring of 100 nodes is adjacent to its 20 nearest neighbors and we
uniformly randomly choose 30$\%$ of the edges to rewire. When rewiring an
edge, we choose between its two ends with equal probability of each, and we
attach this end to a node in the network that we choose uniformly at random.
We add those newly generated edges to the original network. We refer to this
noise type as ‘${+}\textup{WS}$’.

For each corrupted network, we apply NDL to learn a network dictionary We then
use each dictionary to reconstruct the corrupted network. We anticipate that
the reconstructed network is similar to its corresponding original (i.e.,
uncorrupted) network. The reconstruction algorithms output a weighted network
$G_{\textup{recons}}$, where the weight of each edge is our confidence that
the edge is a true edge of that network. For denoising subtractive
(respectively, additive) noise, we classify each non-edge (respectively, each
edge) in a corrupted network as ‘positive’ if its weight in
$G_{\textup{recons}}$ is strictly larger than some threshold $\theta$ and
‘negative’ otherwise. By varying $\theta$, we construct a receiver-operating
characteristic (ROC) curve that consists of points whose horizontal and
vertical coordinates are the false-positive rates and true-positive rates,
respectively. For denoising the ${-}\textup{ER}$ (respectively,
${+}\textup{ER}$ or ${+}\textup{WS}$) noise, one can also infer an optimal
value of $\theta$ for a 50% training set of non-edges (respectively, edges) in
$G$ with known labels and then use this value of $\theta$ to compute
classification measures such as accuracy and precision.

Figure 8.  Denoising Caltech corrupted by $50\%$ additive noise of type ER.
There are 16,656 true (i.e., original) edges and 7,854 false (i.e., added)
edges to correctly classify. As we illustrate in panel a, when a $k$-walk
connects two nodes $x$ and $y$, we define the distance between $x$ and $y$
along the $k$-walk to be the shortest-path distance between $x$ and $y$.
During reconstruction, we sample a sequence of $k$-walks in a network using a
Markov-chain Monte Carlo (MCMC) algorithm (see Algorithm MP in the SI).
Suppose that this sequence is $\mathbf{x}_{1},\ldots,\mathbf{x}_{T}$. We
compute the mean of the distances between $x$ and $y$ along the $k$-walk
$\mathbf{x}_{t}$ for all $t\in\\{1,\ldots,T\\}$ such that $\mathbf{x}_{t}$
connects $x$ and $y$. In panel g, we show the histogram of the mean distances
along the $k$-walks $\mathbf{x}_{1},\ldots,\mathbf{x}_{T}$ between the two
ends of true edges and the two ends of false edges. In panels b–f, we show
histograms of edge weights for various reconstructions using $21$-node latent
motifs (as weighted adjacency matrices) in panels h–l. We learn the 25 latent
motifs in panel h from the corrupted network. We multiply the on-chain edge
weights of these latent motifs by a thinning parameter of $\xi=0.5$ and
$\xi=0$ to obtain the matrices in panels i and j, respectively. We randomly
choose the latent motifs in panel k by drawing each entry of its $k\times k$
weighted adjacency matrix independently and uniformly from $[0,1]$. Setting
the on-chain entries to $0$ of such random latent motifs gives the matrices in
panel l. The classification accuracy (ACC) and precision (PRE) use the optimal
threshold for truncating weighted edges in the reconstruction that we compute
from a uniformly randomly chosen training set of edges (with $50\%$ of the
edges in the observed network). The area under the curve (AUC) refers to the
area under the receiver-operating characteristic (ROC) curve.

Based on our experiments, the network-reconstruction algorithm that we
described earlier seems to mostly be effective when denoising subtractive
noise. However, when denoising additive noise (for both ${+}\textup{ER}$ and
${+}\textup{WS}$), the edge weights in a reconstructed network can result it
ROC curves with an area under the curve (AUC) that is as small as 0.5. For
instance, as we show in the histogram in Figure 8b, when we denoise Caltech
with ${+}\textup{ER}$ noise, the false edges are assigned weights that are
significantly larger than those of the true edges when we reconstruct the
observed network using latent motifs that we learn from the corrupted network
(see Figure 8h). Consequently, we obtain an AUC of 0.5 for our classification.
There is a simple explanation of this result. The Markov-chain Monte Carlo
(MCMC) sampling algorithm, which incrementally evolves a specified $k$-walk,
that we use throughout the denoising process tends to connect false edges
using a smaller number of edges than it uses to connect true edges. In other
words, if there is an edge between nodes $x$ and $y$ in an additively
corrupted network and we uniformly randomly sample a $k$-walk that uses both
$x$ and $y$, then the number of edges between these two nodes along the
sampled $k$-walk tends to be small if the edge between $x$ and $y$ is false
and tends to be large if it is true (see Figure 8a). This indicates that there
are not many ways to connect the two ends of a false edge using a $k$-walk
that avoids using that false edge. Consequently, of the edges between nodes in
a uniformly sampled $k$-walk in an additively corrupted network, false edges
are more likely to appear as on-chain edges than as off-chain edges. When we
consider subtractive noise, an analogous observation holds for true non-edges
and false non-edges.

To address this issue, we modify our network-reconstruction algorithm for the
specific purpose of network denoising by thinning out the on-chain edge
weights of latent motifs prior to the reconstruction. To do this, we multiply
the weights of on-chain edges in the latent motifs and in all sampled
mesoscale patches by a scalar chain-edge thinning parameter $\xi\in[0,1]$. For
instance, the latent motifs in panels i and j in Figure 8 use $\xi=0.5$ and
$\xi=0$, respectively. As we see in the histograms in panels c and d, the
negative edges in the resulting reconstruction have significantly smaller
weights than the positive edges. For example, with $\xi=0$, this results in a
classification AUC of 0.91.

In Figure 9, we compare the performance of our network-denoising method to the
performance of several existing approaches. We use four classical approaches
(the Jaccard index, preferential attachment, the Adamic–Adar index, and a
spectral embedding) [20, 32] and two more recent methods (DeepWalk [46] and
node2vec [50]) that are based on network embeddings. Let $N(x)$ denote the set
of neighbors of node $x$ in a network. For the Jaccard index, preferential
attachment, and the Adamic–Adar index, the confidence score (which plays the
same role as an edge weight in a reconstructed network) that the nodes $x$ and
$y$ are adjacent via a true edge is $|N(x)\cap N(y)|/|N(x)\cup N(y)|$,
$|N(x)|\cdot|N(y)|$, and $\sum_{z\in N(x)\cap N(y)}1/\ln|N(z)|$, respectively.

Figure 9.  Application of our NDL and NDR algorithms to network denoising with
additive and subtractive noise on a variety of networks from empirical data
sets. In our experiments with subtractive noise, we corrupt a network by
removing $50\%$ of its edges uniformly at random. We seek to classify the
removed edges and the non-edges as true edges and false edges, respectively.
In our experiments with additive noise, we corrupt a network by uniformly
randomly adding $50\%$ of the number of its edges or 1,000 random edges for
all but one network (we add 30,000 random edges for H. sapiens) that we
generated using the WS model. We seek to classify the edges and non-edges in
the resulting corrupted network as ‘negative’ (i.e., false edges) or
‘positive’ (i.e., true edges). To perform classification in a network, we
first use NDL to learn latent motifs from a corrupted network and then
reconstruct the networks using NDR to assign a confidence value to each
potential edge. We then use these confidence values to infer the correct
labeling of potential edges in the uncorrupted network. Importantly, we never
use information from the original networks. For each network, we report AUC
scores of the ROC curves. See Figure 13 in the SI for accuracy, precision, and
recall scores.

We show results in the forms of means of the AUC of the ROC curves for five
independent runs of each method. The Adamic–Adar index is not defined for
nodes with self-edges, so we do not use it for networks with self-edges (such
as arXiv, Coronavirus, and H. sapiens). For our method, we apply NDL at scale
$k=21$ with $r\in\\{2,25\\}$ to learn latent motifs from the corrupted
networks and apply NDR with the chain-edge thinning parameter
${\xi}\in\\{0,1\\}$. A thinning parameter of ${\xi}=1$ corresponds to network
reconstruction with unchanged latent motifs and mesoscale patches (see Figure
8h), whereas ${\xi}=0$ removes all of the on-chain edges (see Figure 8j).

The existing methods Jaccard index and preferential attachment seem to work
well for noise types $+\textup{ER}$ and $-\textup{ER}$. For instance, the high
performance of preferential attachment for Coronavirus with noise type
$+\textup{ER}$ indicates that the product of the degrees of two nodes that are
attached via a true edge is larger than the product of the degrees of two
nodes that are attached via a false else. To illustrate why this is the case,
consider a star network with five leaves (i.e., degree-1 nodes). In this
network, a uniformly randomly chosen non-edge is always attached to two
degree-1 nodes, but a uniformly randomly chosen edge is always attached to one
degree-5 node (i.e., the center node) and one degree-1 node. As in this star
example, all real-world networks that we consider are sparse (with edge
densities of at most $0.057$, which is the density for Caltech). One can
approximately uniformly randomly choose a non-edge from such sparse networks
by choosing two nodes uniformly at random. However, as in the star example,
uniformly choosing an edge may induce a bias towards choosing at least one
large-degree node.

Our method outperforms all of the other examined methods for denoising
$+\textup{WS}$ structured noise. For example, when we add 1,000 false edges
that we generate from the WS model to Coronavirus (which has 2,463 true
edges), our method yields an AUC of 0.94. We obtain the second best AUC (of
only 0.61) using preferential attachment. There are at least two confounding
factors to consider in the high accuracy of our network denoising method.
First, as we noted before, false edges are more likely to appear as on-chain
edges than as off-chain edges in a subgraph that one samples by a uniformly
random $k$-walk. Consequently, the mean distance between two nodes along
$k$-walks tends to be larger when they form a true edge than when they form a
false edge (as shown in Figure 8a). This distances already somewhat
distinguishes the true edges and the false edges. Consequently, we see in
Figure 8k,l that we can reasonably successfully denoise the network Caltech
with additive noise of type ${+}\textup{ER}$ using randomized latent motifs
(in which each we draw each entry of their weighted adjacency matrices
independently and uniformly from $[0,1]$), regardless of whether the on-chain
thinning parameter is $\xi=1$ (see Figure 8k) or $\xi=0$ (see Figure 8l).
Second, knowing the exact subgraph mesoscale structure (i.e., the latent
motifs of the original network, which by definition does not have false edges)
improves the denoising accuracy, which one can infer directly from the
observed network (which has false edges) because of the low-rank nature of the
examined mesoscale structures of the social and PPI networks (see Figure 5a).
Specifically, because NDL learns a small number of latent motifs that are able
to successfully give an approximate basis for all mesoscale patches, they
should not be affected significantly by noise.

### Conclusions

We introduced a mesoscale network structure, which we call latent motifs, that
consists of $k$-node subgraphs that are building blocks of all connected
$k$-node subgraphs of a network. In contrast to motifs [9], which are used
frequently to capture overrepresented $k$-node subgraphs in a network,
nonnegative linear combinations of our latent motifs approximate $k$-node
subgraphs that are induced by uniformly random $k$-paths in a network. We also
established theoretically that one can approximate an entire network
accurately if one has a dictionary of latent motifs that can accurately
approximate the network’s mesoscale structures.

Our computational experiments in Figures 2, 3, and 4 demonstrated that latent
motifs can have distinctive network structures at various mesoscales. Our
computational experiments in Figures 5, 8, and 9 illustrated that various
social, collaboration and PPI networks have low-rank [39] mesoscale
structures, in the sense that a few (e.g., $r=25$, but see Figure 5 for other
choices of $r$) latent motifs that we learn using NDL are able to reconstruct,
infer, and denoise the edges in an entire network by employing our NDR
algorithm. We hypothesize that such low-rank mesoscale structures are a common
feature of networks beyond the social, collaboration, and PPI networks that we
examined. As we have illustrated in this paper, one can leverage mesoscale
structures to perform important tasks like network denoising, so it is
important in future studies to explore the level of generality of our
insights.

Notions of motifs have been extended to more complicated network structures,
such as temporal networks (in which nodes, edges, and edge weights can change
with time) [53] and multilayer networks (in which, e.g., nodes can be adjacent
via multiple types of relationships) [51]. We did not consider latent motifs
in these generalized network structures in our work, but it is clearly of
interest to extend our approach and algorithms to these situations.

## Acknowledgements

HL was supported by the National Science Foundation through grants 2206296 and
2010035. JV was supported by the National Science Foundation through grant
1740325.

LEARNING LOW-RANK LATENT MESOSCALE STRUCTURES IN NETWORKS

## References

  * [1] Paul Erdős and Alfréd Rényi  “On random graphs. I”  In _Publicationes Mathematicae_ 6.18, 1959, pp. 290–297
  * [2] William M. Rand  “Objective criteria for the evaluation of clustering methods”  In _Journal of the American Statistical Association_ 66.336 Taylor & Francis Group, 1971, pp. 846–850
  * [3] Paul W. Holland, Kathryn Blackmond Laskey and Samuel Leinhardt  “Stochastic blockmodels: First steps”  In _Social Networks_ 5.2 Elsevier, 1983, pp. 109–137
  * [4] Robert Tibshirani  “Regression shrinkage and selection via the lasso”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 58.1 Wiley Online Library, 1996, pp. 267–288
  * [5] Duncan J. Watts and Steven H. Strogatz  “Collective dynamics of ‘small-world’ networks”  In _Nature_ 393.6684 Nature Publishing Group, 1998, pp. 440–442
  * [6] Albert-László Barabási and Réka Albert  “Emergence of scaling in random networks”  In _Science_ 286.5439 American Association for the Advancement of Science, 1999, pp. 509–512
  * [7] Daniel D. Lee and H. Sebastian Seung  “Learning the parts of objects by non-negative matrix factorization”  In _Nature_ 401.6755 Nature Publishing Group, 1999, pp. 788–791
  * [8] Daniel D. Lee and H. Sebastian Seung  “Algorithms for non-negative matrix factorization”  In _Advances in Neural Information Processing Systems_ , 2001, pp. 556–562
  * [9] Ron Milo et al.  “Network motifs: Simple building blocks of complex networks”  In _Science_ 298.5594 American Association for the Advancement of Science, 2002, pp. 824–827
  * [10] Gavin C. Conant and Andreas Wagner  “Convergent evolution of gene circuits”  In _Nature Genetics_ 34.3 Nature Publishing Group, 2003, pp. 264–266
  * [11] Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani  “Least angle regression”  In _The Annals of Statistics_ 32.2 Institute of Mathematical Statistics, 2004, pp. 407–499
  * [12] Paul Glasserman  “Monte Carlo methods in financial engineering”  Springer, 2004
  * [13] Nadav Kashtan, Shalev Itzkovitz, Ron Milo and Uri Alon  “Efficient sampling algorithm for estimating subgraph concentrations and detecting network motifs”  In _Bioinformatics_ 20.11 Oxford University Press, 2004, pp. 1746–1758
  * [14] Olaf Sporns, Rolf Kötter and Karl J. Friston  “Motifs in brain networks”  In _PLoS Biology_ 2.11 Public Library of Science San Francisco, USA, 2004, pp. e369
  * [15] Michael Elad and Michal Aharon  “Image denoising via sparse and redundant representations over learned dictionaries”  In _IEEE Transactions on Image Processing_ 15.12 IEEE, 2006, pp. 3736–3745
  * [16] Jure Leskovec and Christos Faloutsos  “Sampling from large graphs”  In _Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2006, pp. 631–636
  * [17] Sebastian Wernicke  “Efficient detection of network motifs”  In _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ 3.4 IEEE, 2006, pp. 347–359
  * [18] Uri Alon  “Network motifs: Theory and experimental approaches”  In _Nature Reviews Genetics_ 8.6 Nature Publishing Group, 2007, pp. 450–461
  * [19] Honglak Lee, Alexis Battle, Rajat Raina and Andrew Y. Ng  “Efficient sparse coding algorithms”  In _Advances in Neural Information Processing Systems_ , 2007, pp. 801–808
  * [20] David Liben-Nowell and Jon Kleinberg  “The link-prediction problem for social networks”  In _Journal of the American Society for Information Science and Technology_ 58.7 Wiley Online Library, 2007, pp. 1019–1031
  * [21] Julien Mairal, Michael Elad and Guillermo Sapiro  “Sparse representation for color image restoration”  In _IEEE Transactions on Image Processing_ 17.1 IEEE, 2007, pp. 53–69
  * [22] Vincent D. Blondel, Jean-Loup Guillaume, Renaud Lambiotte and Etienne Lefebvre  “Fast unfolding of communities in large networks”  In _Journal of Statistical Mechanics: Theory and Experiment_ 2008.10 IOP Publishing, 2008, pp. P10008
  * [23] Julien Mairal, Michael Elad and Guillermo Sapiro  “Sparse learned representations for image restoration”  In _Proceedings of the 4th World Conference of the International Association for Statistical Computing_ , 2008, pp. 118  IASC
  * [24] Krzysztof Juszczyszyn, Przemyslaw Kazienko and Bogdan Gabrys  “Temporal changes in local topology of an email-based social network”  In _Computing and Informatics_ 28.6, 2009, pp. 763–779
  * [25] Julien Mairal et al.  “Non-local sparse models for image restoration”  In _IEEE 12th International Conference on Computer Vision_ , 2009, pp. 2272–2279  IEEE
  * [26] Gabriel Peyré  “Sparse modeling of textures”  In _Journal of Mathematical Imaging and Vision_ 34.1 Springer, 2009, pp. 17–31
  * [27] Mason A. Porter, Jukka-Pekka Onnela and Peter J. Mucha  “Communities in networks”  In _Notices of the American Mathematical Society_ 56.9, 2009, pp. 1082–1097, 1164–1166
  * [28] Richard T. Durrett  “Probability: Theory and Examples”, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge, UK: Cambridge University Press, 2010, pp. 428
  * [29] Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro  “Online learning for matrix factorization and sparse coding”  In _Journal of Machine Learning Research_ 11, 2010, pp. 19–60
  * [30] Takaaki Ohnishi, Hideki Takayasu and Misako Takayasu  “Network motifs in an inter-firm network”  In _Journal of Economic Interaction and Coordination_ 5.2 Springer, 2010, pp. 171–180
  * [31] Jason M. K. Rip, Kevin S. McCann, Denis H. Lynn and Sonia Fawcett  “An experimental test of a fundamental food web motif”  In _Proceedings of the Royal Society B: Biological Sciences_ 277.1688 The Royal Society, 2010, pp. 1743–1749
  * [32] Mohammad Al Hasan and Mohammed J. Zaki  “A survey of link prediction in social networks”  In _Social Network Data Analytics_ Springer, 2011, pp. 243–275
  * [33] Aditya Krishna Menon and Charles Elkan  “Link prediction via matrix factorization”  In _Machine Learning and Knowledge Discovery in Databases_ Heidelberg, Germany: Springer-Verlag, 2011, pp. 437–452
  * [34] Veronica Red, Eric D. Kelsic, Peter J. Mucha and Mason A. Porter  “Comparing community structure to characteristics in online collegiate social networks”  In _SIAM Review_ 53, 2011, pp. 526–543
  * [35] Lei Tang and Huan Liu  “Leveraging social media networks for classification”  In _Data Mining and Knowledge Discovery_ 23.3 Springer, 2011, pp. 447–478
  * [36] Roger A. Horn and Charles R. Johnson  “Matrix Analysis”  Cambridge, UK: Cambridge University Press, 2012
  * [37] Jure Leskovec and Julian J. Mcauley  “Learning to discover social circles in ego networks”  In _Advances in Neural Information Processing Systems_ , 2012, pp. 539–547
  * [38] László Lovász  “Large Networks and Graph Limits” 60, Colloquium Publications  Providence, RI, USA: American Mathematical Society, 2012, pp. 475
  * [39] Ivan Markovsky and Konstantin Usevich  “Low Rank Approximation”  Heidelberg, Germany: Springer-Verlag, 2012
  * [40] Sean P Meyn and Richard L Tweedie  “Markov Chains and Stochastic Stability”  Heidelberg, Germany: Springer-Verlag, 2012
  * [41] Jukka-Pekka Onnela et al.  “Taxonomies of networks from community structure”  In _Physical Review E_ 86.3 APS, 2012, pp. 036104
  * [42] Amanda L. Traud, Peter J. Mucha and Mason A. Porter  “Social structure of Facebook networks”  In _Physica A_ 391.16 Elsevier, 2012, pp. 4165–4180
  * [43] Julien Mairal  “Stochastic majorization-minimization algorithms for large-scale optimization”  In _Advances in Neural Information Processing Systems_ , 2013, pp. 2283–2291
  * [44] Tomas Mikolov et al.  “Distributed Representations of Words and Phrases and their Compositionality”  In _Advances in Neural Information Processing Systems_ 26 Curran Associates, Inc., 2013
  * [45] Xu Hong-lin, Yan Han-bing, Gao Cui-fang and Zhu Ping  “Social network analysis based on network motifs”  In _Journal of Applied Mathematics_ 2014 Hindawi, 2014, pp. 874708
  * [46] Bryan Perozzi, Rami Al-Rfou and Steven Skiena  “DeepWalk: Online learning of social representations”  In _Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2014, pp. 701–710
  * [47] Konstantin Ristl, Sebastian J. Plitzko and Barbara Drossel  “Complex response of a food-web module to symmetric and asymmetric migration between several patches”  In _Journal of Theoretical Biology_ 354 Elsevier, 2014, pp. 54–59
  * [48] Lucas G. S. Jeub et al.  “Think locally, act locally: Detection of small, medium-sized, and large communities in large networks”  In _Physical Review E_ 91 American Physical Society, 2015, pp. 012821
  * [49] Santo Fortunato and Darko Hric  “Community detection in networks: A user guide”  In _Physics Reports_ 659 Elsevier, 2016, pp. 1–44
  * [50] Aditya Grover and Jure Leskovec  “node2vec: Scalable feature learning for networks”  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2016, pp. 855–864
  * [51] Federico Battiston, Vincenzo Nicosia, Mario Chavez and Vito Latora  “Multilayer motif analysis of brain networks”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 27.4 AIP Publishing LLC, 2017, pp. 047404
  * [52] David A Levin and Yuval Peres  “Markov Chains and Mixing Times”  Providence, RI, USA: American Mathematical Society, 2017
  * [53] Ashwin Paranjape, Austin R Benson and Jure Leskovec  “Motifs in temporal networks”  In _Proceedings of the tenth ACM international conference on web search and data mining_ , 2017, pp. 601–610
  * [54] Ankit N. Khambhati, Ann E. Sizemore, Richard F. Betzel and Danielle S. Bassett  “Modeling and interpreting mesoscale network dynamics”  In _NeuroImage_ 180 Elsevier, 2018, pp. 337–349
  * [55] Mark E. J. Newman  “Networks”  Oxford, UK: Oxford University Press, 2018
  * [56] Frank W. Takes, Walter A. Kosters, Boyd Witte and Eelke M. Heemskerk  “Multiplex network motifs as building blocks of corporate networks”  In _Applied Network Science_ 3.1 Springer, 2018, pp. 1–22
  * [57] Fernanda B. Correia, Edgar D. Coelho, José L. Oliveira and Joel P. Arrais  “Handling noise in protein interaction networks”  In _BioMed Research International_ 2019 Hindawi, 2019, pp. 1–13
  * [58] Katja Kovács et al.  “Network-based prediction of protein interactions”  In _Nature Communications_ 10.1 Nature Publishing Group, 2019, pp. 1240
  * [59] Hanbaek Lyu, Facundo Memoli and David Sivakoff  “Sampling random graph homomorphisms and applications to network data analysis”  In _arXiv:1910.09483_ , 2019
  * [60] Rose Oughtred et al.  “The BioGRID interaction database: 2019 update”  In _Nucleic Acids Research_ 47.D1 Oxford University Press, 2019, pp. D529–D541
  * [61] David E. Gordon et al.  “A SARS-CoV-2 protein interaction map reveals targets for drug repurposing”  In _Nature_ 583 Nature Publishing Group, 2020, pp. 1–13
  * [62] Roger Guimerà  “One model to rule them all in network science?”  In _Proceedings of the National Academy of Sciences of the United States of America_ 117.41 National Acad Sciences, 2020, pp. 25195–25197
  * [63] Jure Leskovec and Andrej Krevl  “SNAP Datasets: Stanford Large Network Dataset Collection” Available at http://snap.stanford.edu/data, 2020
  * [64] Hanbaek Lyu, Deanna Needell and Laura Balzano  “Online matrix factorization for Markovian data and applications to network dictionary learning”  In _Journal of Machine Learning Research_ 21, 2020, pp. 1–49
  * [65] C. Seshadhri, Aneesh Sharma, Andrew Stolman and Ashish Goel  “The impossibility of low-rank representations for triangle-rich complex networks”  In _Proceedings of the National Academy of Sciences of the United States of America_ 117.11 National Academy of Sciences, 2020, pp. 5631–5637
  * [66] theBiogrid.org  “Coronavirus PPI network” Retrieved from https://wiki.thebiogrid.org/doku.php/covid (downloaded 24 July 2020, Ver. 3.5.187.tab3), 2020
  * [67] theBiogrid.org  “Homo sapiens PPI network” Retrieved from https://wiki.thebiogrid.org/doku.php/covid (downloaded 24 July 2020, Ver. 3.5.180.tab2), 2020
  * [68] Marco Bressan  “Efficient and near-optimal algorithms for sampling connected subgraphs”  In _Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing_ , STOC 2021  Virtual, Italy: Association for Computing Machinery, 2021, pp. 1132–1143
  * [69] Alice C. Schwarze and Mason A. Porter  “Motifs for processes on networks”  In _SIAM Journal on Applied Dynamical Systems_ 20.4, 2021, pp. 2516–2557
  * [70] Tao Zhou  “Progresses and challenges in link prediction”  In _iScience_ 24.11, 2021, pp. 103217

## Methods

We briefly discuss our algorithms for network dictionary learning (NDL) and
network denoising and reconstruction (NDR). We also provide a detailed
description of the employed real-world and synthetic networks.

In this section, we restrict our discussion to networks that one can represent
as a graph $G=(V,E)$ with a node-set $V$ and an edge set $E$ without directed
edges or multi-edges (but possibly with self-edges). In our SI, we give an
extended discussion that applies to more general types of networks.
Specifically, in that discussion, we no longer restrict edges to have binary
weights; instead, the weights may be continuous-valued. See Section B.1 of the
SI.

### Motif sampling and mesoscale patches of networks

Natural candidates for the mesoscale patches of a network are its connected
$k$-node subgraphs. These subgraphs have $k$ nodes that inherit their
adjacency structures from the original networks from which we draw them. It is
convenient to consider the $k\times k$ adjacency matrices of these subgraphs,
as it allows us to perform computations on the space of subgraphs. However, we
need to address two issues to do this. First, because the same $k$-node
subgraph can have multiple (specifically, $k!$) different representations as
an adjacency matrix (depending on the ordering of its nodes), we need an
unambiguous way to choose an ordering of its nodes. Second, because most real-
world networks are sparse [55], independently choosing a set of $k$ nodes from
a network may yield only a few edges and may thus often result in a
disconnected subgraph. Therefore, we need an efficient sampling algorithm to
guarantee that we obtain connected $k$-node subgraphs when we sample from
sparse networks.

We employ an approach that is based on motif sampling [59] both to choose an
ordering of the nodes in a $k$-node subgraph and to ensure that we sample
connected subgraphs from sparse networks. The key idea is to consider the
random $k$-node subgraph that we obtain by first sampling a copy of a
‘template’ subgraph uniformly at random and including all edges between the
sampled nodes of the network. We take such a template to be a ‘$k$-path’. A
sequence $\mathbf{x}=(x_{1},\ldots,x_{k})$ of $k$ (not necessarily distinct)
nodes is a $k$-walk if $x_{i}$ and $x_{i+1}$ are adjacent for all
$i\in\\{1,\ldots,k-1\\}$. A $k$-walk is a $k$-path if all nodes in the walk
are distinct (see Figure 1). For each $k$-path
$\mathbf{x}=(x_{1},\ldots,x_{k})$, we define the corresponding mesoscale patch
of a network to be the $k\times k$ matrix $A_{\mathbf{x}}$ such that
$A_{\mathbf{x}}(i,j)=1$ if nodes $x_{i}$ and $x_{j}$ are adjacent and
$A_{\mathbf{x}}(i,j)=0$ if they are not adjacent. This is the adjacency matrix
of the $k$-node subgraph of the network with nodes $x_{1},\ldots,x_{k}$.  One
can use one of the Markov-chain Monte Carlo (MCMC) algorithms for motif
sampling from [59] to efficiently and uniformly randomly sample a $k$-walk
from a sparse network. By only accepting samples in which the $k$-walk has $k$
distinct nodes (i.e., so that it is $k$-path), we efficiently sample a
uniformly random $k$-path from a network, as long as $k$ is not too large. If
$k$ is too large, one has to ‘reject’ too many samples of $k$-walks that are
not $k$-paths. The expected number of rejected samples is approximately the
number of $k$-walks divided by the number of $k$-paths. The number of
$k$-walks in a network grows monotonically with $k$, but the number of
$k$-paths can decrease with $k$ and it becomes identically $0$ if $k$ exceeds
the shortest-path diameter of the network. (See Section C of the SI for a
detailed discussion.) Consequently, by repeatedly sampling $k$-paths
$\mathbf{x}$, we obtain a data set of mesoscale patches $A_{\mathbf{x}}$ of a
network.

### Algorithm for network dictionary learning

We now present the basic structure of the algorithms that we employ for
network dictionary learning (NDL) [64]. Suppose that we compute all possible
$k$-paths $\mathbf{x}_{1},\dots,\mathbf{x}_{M}$ and their corresponding
mesoscale patches $A_{\mathbf{x}_{t}}$ (which are $k\times k$ binary matrices)
for $t\in\\{1,\ldots,M\\}$ of a network. We column-wise vectorize (i.e., we
place the second column underneath the first column and so on; see Algorithm
A4 in the SI) each of these $k\times k$ mesoscale patches to obtain a
$k^{2}\times M$ data matrix $X$. We then apply nonnegative matrix
factorization (NMF) [7] to obtain a $k^{2}\times r$ nonnegative matrix $W$ for
some fixed integer $r\geq 1$ to yield an approximate factorization $X\approx
WH$ for some nonnegative matrix $H$. From this procedure, we approximate each
column of $X$ by the nonnegative linear combination of the $r$ columns of $W$;
its coefficients are the entries of the corresponding column of $H$. If we let
$\mathcal{L}_{i}$ be the $k\times k$ matrix that we obtain by reshaping the
$i^{\textup{th}}$ column of $W$ (using Algorithm A5 in SI), then
$\mathcal{L}_{1},\ldots,\mathcal{L}_{r}$ are the learned latent motifs; they
form a network dictionary. The set of these latent motifs is an approximate
basis (but not a subset) of the set
$A_{\mathbf{x}_{1}},\dots,A_{\mathbf{x}_{M}}$ of mesoscale patches. For
instance, latent motifs have entries that take continuous values between 0 and
1, but mesoscale patches have binary entries. We can regard each
$\mathcal{L}_{i}$ as the $k$-node weighted network with node set
$\\{1,\ldots,k\\}$ and weighted adjacency matrix $\mathcal{L}_{i}$. See Figure
2 for an illustration of latent motifs as weighted networks.

The scheme in the paragraph above requires us to store all possible mesoscale
patches of a network, entailing a memory requirement that is at least of order
$k^{2}M$, where $M$ denotes the total number of mesoscale patches in a
network. Because $M$ scales with the size of the network from which we sample
subgraphs, we need unbounded memory to handle arbitrarily large networks. To
address this issue, Algorithm NDL implements the above scheme in the setting
of ‘online learning’, where subsets (so-called ‘minibatches’) of data arrive
in a sequential manner and one does not store previous subsets of the data
before processing new subsets. Specifically, at each iteration
$t=1,2,\ldots,T$, we process a sample matrix $X_{t}$ that is smaller than the
full matrix $X$ and includes only $N\ll M$ mesoscale patches, where one can
take $N$ to be independent of the network size. Instead of a standard NMF
algorithm for a fixed matrix [8], we use an ‘online’ NMF algorithm [29, 64]
that one can use on sequences of matrices, where the intermediate dictionary
matrices $W_{t}$ that we obtain by factoring the sample matrix $X_{t}$
typically improve as we iterate (see [29, 64]). In Algorithm NDL in the SI, we
give a complete implementation of the NDL algorithm.

### Algorithm for network denoising and reconstruction

Suppose that we have an image $\gamma$ of size $k\times k$ pixels and a set of
‘basis images’ $\beta_{1},\ldots,\beta_{r}$ of the same size. We can
‘reconstruct’ the image $\gamma$ using the basis images
$\beta_{1},\ldots,\beta_{r}$ by finding nonnegative coefficients
$a_{1},\ldots,a_{r}$ such that the linear combination
$\hat{\gamma}=a_{1}\beta_{1}+\cdots+a_{r}\beta_{r}$ is as close as possible to
$\gamma$. The basis images determine what shapes and colors of the original
image to capture in the reconstruction $\hat{\gamma}$. In the standard
pipeline for image denoising and reconstruction [15, 23, 25], one assumes that
the size $k\times k$ of the square ‘patches’ is much smaller than the full
image $\gamma$. One can then sample a large number of $k\times k$ overlapping
patches $\gamma_{1},\ldots,\gamma_{M}$ of the image $\gamma$ and obtain their
best linear approximations $\hat{\gamma}_{1},\dots,\hat{\gamma}_{M}$ of them
using the basis images $\beta_{1},\ldots,\beta_{r}$. Because the $k\times k$
patches $\gamma_{1},\ldots,\gamma_{M}$ overlap, each pixel $(I,J)$ of $\gamma$
can occur in multiple instances of $\gamma_{1},\ldots,\gamma_{M}$. Therefore,
we take the mean of the corresponding values in the mesoscale reconstruction
$\hat{\gamma}_{0}$ as the value of the pixel $(I,J)$ in the global
reconstruction $\hat{\gamma}$.

Figure 10.  a image Woman with a Parasol - Madame Monet and Her Son (Claude
Monet, 1875). We show reconstructions of it in b and c. The image in panel b
is a reconstruction of the image in panel a using the dictionary with 25 basis
images of size $21\times 21$ pixels in panel d. We uniformly randomly choose
the color of each pixel from all possible colors (which we represent as
vectors in $[0,256]^{3}$ for red–green–blue (RGB) weights). The image in panel
d is a reconstruction of the image in panel a using the dictionary with 25
basis images of size $21\times 21$ pixels in panel e. We learn this basis from
the image in panel f using NMF. The image in panel f is from the collection
Die Graphik Ernst Ludwig Kirchners bis 1924, von Gustav Schiefler Band I bis
1916 (Accession Number 2007.141.9, Ernst Ludwig Kirchner, 1926). [We use the
images in panels a and f with permission from the National Gallery of Art in
Washington, DC, USA.]

As an illustration, consider reconstructing the color image in Figure 10a in
two ways, which yield the images in panels b and c. In panel d, we show the
dictionary with 25 basis images of size $21\times 21$ pixels. We uniformly
randomly choose the color of each pixel from all possible colors (which we
represent as vectors in $[0,256]^{3}$ for red–green–blue (RGB) weights). The
basis images do not include any information about the original image in panel
a, so the linear approximation of the $21\times 21$ mesoscale patches of the
image in panel a using the basis images in panel d may be inaccurate. However,
when we reconstruct the entire image from panel a using the basis images in
panel d, we do observe some basic geometric information from the original
image. We show the image that results from this reconstruction in panel b.
Importantly, the image reconstruction in b uses both the basis images and the
original image that one seeks to reconstruct. Unfortunately, the colors have
averaged out to be neutral, so the reconstructed image is monochrome. Using a
smaller (e.g., $5\times 5$) randomly generated (with each pixel again taking
an independently and uniformly chosen color) basis-image set for
reconstruction results in a similar monochrome (but sharper) reconstructed
image. Notably, one can learn the basis images in panel e from the image in
panel f using nonnegative matrix factorization (NMF) [7]. The image in f is in
black and white, and thus so are the images in e. The corresponding
reconstruction in panel c has nicely captured shapes from the original image,
although we have lost the color information in the original image in panel a
and the reconstruction in panel d is thus in black and white.

A network analog of the above patch-based image reconstruction proceeds as
follows. Given a network $G=(V,E)$ and latent motifs
$(\mathcal{L}_{1},\ldots,\mathcal{L}_{r})$ (which we do not necessarily
compute from $G$; see Figure 10), we obtain a weighted network
$G_{\textup{recons}}$ using the same node set $V$ and a weighted adjacency
matrix $A_{\textup{recons}}:V^{2}\rightarrow\mathbb{R}$. To do this, we first
uniformly randomly sample a large number $T$ of $k$-walks
$\mathbf{x}_{1},\ldots,\mathbf{x}_{T}:\\{1,\ldots,k\\}\rightarrow V$ in $G$.
We then determine the corresponding mesoscale patches
$A_{\mathbf{x}_{1}},\ldots,A_{\mathbf{x}_{T}}$ of $G$, which are $k\times k$
unweighted matrices. We then approximate each mesoscale patch
$A_{\mathbf{x}_{t}}$ by a nonnegative linear combination
$\hat{A}_{\mathbf{x}_{t}}$ of the latent motifs $\mathcal{L}_{i}$. We seek to
‘replace’ each $A_{\mathbf{x}_{t}}$ with $\hat{A}_{\mathbf{x}_{t}}$ to
construct the weighted adjacency matrix $A_{\textup{recons}}$. To do this, we
define $A_{\textup{recons}}(x,y)$ for each $x,y\in V$ as the mean of
$\hat{A}_{\mathbf{x}_{t}}(a,b)$ over all $t\in\\{1,\ldots,T\\}$ and all
$a,b\in\\{1,\ldots,k\\}$ such that $\mathbf{x}_{t}(a)={x}$ and
$\mathbf{x}_{t}(b)={y}$. Our network denoising and reconstruction (NDR)
algorithm (see Algorithm NDR) builds on this idea. See Section E of the SI for
more details.

### Comparison of our work with the prior research in [64]

Recently, Lyu et al. [64] proposed a preliminary approach for the algorithms
that we propose in the present work — the NDL algorithm with $k$-walk sampling
and the network-reconstruction — as an application to showcase a theoretical
result about the convergence of online NMF for data samples that are not
independently and identically distributed (IID) [64, Thm. 1]. A notable
limitation of the NDL algorithm in [64] is that one cannot interpret the
elements of a network dictionary as latent motifs and one thus cannot
associate them directly with mesoscale structures in networks. Additionally,
Lyu et al. [64] does not include any theoretical analysis of either the
convergence or the correctness of network reconstruction, so it was unclear
from that work whether or not one can reconstruct a network using the low-rank
mesoscale structure that is encoded in a network dictionary. Moreover, one
cannot use the network-reconstruction algorithm that was proposed in [64] to
denoise additive noise unless one knows in advance that the noise is additive
(rather than subtractive) before denoising (see [64, Rmk. 4.]). In the present
paper, we build substantially on the research in [64] and provide a much more
complete computational and theoretical framework to analyze low-rank mesoscale
structure of networks. In particular, we overcome all of the aforementioned
limitations. In Table 1, we summarize the key differences between the present
work and the recent paper by Lyu et al. [64].

NDL | Sampling | Latent motifs | Convergence | Efficient MCMC
---|---|---|---|---
Lyu et al. [64] |  $k$-walks | ✗ | Non-bipartite networks | ✗
The present work | $\begin{matrix}\textup{$k$-paths}\\\ \textup{(Alg.\leavevmode\nobreak\ \ref{alg:motif_inj}) }\end{matrix}$ | ✓ | $\begin{matrix}\textup{Non-bipartite and}\\\ \textup{{bipartite networks}}\\\ \textup{(Thm.\leavevmode\nobreak\ \ref{thm:NDL}, Thm.\leavevmode\nobreak\ \ref{thm:NDL2}) }\end{matrix}$ | $\begin{matrix}\checkmark\\\ \textup{(Prop.\leavevmode\nobreak\ \ref{prop:approximate_pivot}) }\end{matrix}$
|  |  |  |

NDR | Sampling | Reconstruction | Denoising | Convergence | Error bound
---|---|---|---|---|---
Lyu et al. [64] |  $k$-walks | ✓ | ✗ | ✗ | ✗
The present work | $\begin{matrix}\textup{$k$-walks}\\\ \textup{$k$-paths}\end{matrix}$ | ✓ | ✓ | $\begin{matrix}\checkmark\\\ \textup{ (Thm.\leavevmode\nobreak\ \ref{thm:NR}) }\end{matrix}$ | $\begin{matrix}\checkmark\\\ \textup{ (Thm.\leavevmode\nobreak\ \ref{thm:NR}) }\end{matrix}$
|  |  |  |  |

Table 1.  A comparison of the contributions in the present work with those in
Lyu et al. [64]. In the SI, we give the statements and proofs of the
proposition and theorems in this table.

The most significant theoretical advance that we make in the present paper
concerns the relationship between global reconstruction error and an
approximation error of mesoscale patches by latent motifs, with an explicit
dependence on the number of nodes in subgraphs at the mesoscale that we state
in Theorem G.10(iii) of our SI. Informally, Theorem G.10(iii) states that one
can accurately reconstruct a network if one has a dictionary of latent motifs
that can accurately approximate mesoscale patches of a network. We illustrate
this theoretical result with supporting experiments in Figure 7. A crucial
part of our proof of Theorem G.10 (iii) is our use of an explicit formula for
the limiting reconstructed network as the number of iterations that one uses
for network reconstruction tends to infinity. We state this explicit formula
for the limiting reconstructed network, which is also a novel contribution of
the present paper, in Theorem G.10(ii).

We now elaborate on the use of $k$-path sampling in our new NDL algorithm to
ensure the interpretability of the network dictionary elements as latent
motifs. The NDL algorithm in [64] used the $k$-walk motif-sampling algorithm
of [59]. That algorithm samples a sequence of $k$ nodes (which are not
necessarily distinct) in which the $i$th node is adjacent to the $(i+1)$th
node for all $i\in\\{1,\ldots,k-1\\}$. The $k$-walks that sample $k\times k$
subgraph adjacency matrices can have overlapping nodes, so some of the
$k\times k$ adjacency matrices can correspond to subgraphs with fewer than $k$
nodes. If there are a large number of such subgraphs of a network, then the
$k$-node latent motifs that one learns from the set of subgraph adjacency
matrices can have misleading patterns that may not exist in any $k$-node
subgraph of the network. This situation occurs in the network CORONAVIRUS PPI,
where one obtains clusters of large-degree nodes from the learned latent
motifs if one uses $k$-walk sampling. However, this misleading result arises
from the $k$-walk visiting the same large-degree node many times, rather than
because $k$ distinct nodes of the network actually exhibit this type of
connection (see Figure 11). To resolve this issue in the present paper, during
the dictionary-learning phase, we combine MCMC $k$-walk sampling with
rejection sampling so that we use only $k$-walks with $k$ distinct nodes
(i.e., we use $k$-paths). Consequently, we now learn the $k$-node latent
motifs only from $k\times k$ adjacency matrices that correspond to $k$-node
subgraphs of the network. This guarantees that any network structure (e.g.,
large-degree nodes, communities, and so on) in the latent motifs must also
exist in the network at scale $k$.

Figure 11.  comparison of subgraphs of Coronavirus PPI that are induced by
node sets that we sample using a uniformly random $k$-paths and b uniformly
random $k$-walks with $k=10$. We also compare the network dictionary with
$r=9$ latent motifs of Coronavirus PPI that we determine using c the NDL
algorithm (see Algorithm NDL) in the present work to d the network dictionary
that we determine using the NDL algorithm from [64]. We also show the weighted
adjacency matrices of the latent motifs. The $10$-walks in the network tend to
visit the same nodes many times. Consequently, one cannot regard the $10\times
10$ mesoscale patches that correspond to those walks as the adjacency matrices
of $k$-node subgraphs of the network. Additionally, the networks in the
network dictionary in panel d have clusters of several large-degree nodes,
even though the original network does not possess such mesoscale structures.

### Data sets

We examine the following eight real-world networks:

  1. (1)

Caltech: This connected network, which is part of the Facebook100 data set
[42] (and which was studied previously as part of the Facebook5 data set
[34]), has 762 nodes and 16,651 edges. The nodes represent users in the
Facebook network of Caltech on one day in fall 2005, and the edges encode
Facebook ‘friendships’ between these accounts.

  2. (2)

MIT: This connected network, which is part of the Facebook100 data set [42],
has 6,402 nodes and 251,230 edges. The nodes represent user accounts in the
Facebook network of MIT on one day in fall 2005, and the edges encode Facebook
‘friendships’ between these accounts.

  3. (3)

UCLA: This connected network, which is part of the Facebook100 data set [42],
has 20,453 nodes and 747,604 edges. The nodes represent user accounts in the
Facebook network of UCLA on one day in fall 2005, and the edges encode
Facebook ‘friendships’ between these accounts.

  4. (4)

Harvard: This connected network, which is part of the Facebook100 data set
[42], has 15,086 nodes and 824,595 edges. The nodes represent user accounts in
the Facebook network of Harvard on one day in fall 2005, and the edges
represent Facebook ‘friendships’ between these accounts.

  5. (5)

SNAP Facebook (with the shorthand SNAP FB) [37]: This connected network has
4,039 nodes and 88,234 edges. This network is a Facebook network that has been
used as an example for edge inference [50]. The nodes represent user accounts
in the Facebook network on one day in 2012, and the edges represent Facebook
‘friendships’ between these accounts.

  6. (6)

arXiv ASTRO-PH (with the shorthand arXiv) [63, 50]: This network has 18,722
nodes and 198,110 edges. Its largest connected component has 17,903 nodes and
197,031 edges. We use the complete network in our experiments. This network is
a collaboration network between authors of astrophysics papers that were
posted to the arXiv preprint server. The nodes represent scientists and the
edges indicate coauthorship relationships. This network has 60 self-edges;
these edges encode single-author papers.

  7. (7)

Coronavirus PPI (with the shorthand Coronavirus): This connected network is
curated by theBiogrid.org [60, 66, 61] from 142 publications and preprints. It
has 1,536 proteins that are related to coronaviruses and 2,463 protein–protein
interactions (in the form of physical contacts) between them. This network is
the largest connected component of the Coronavirus PPI network that we
downloaded on 24 July 2020; in total, there are 1,555 proteins and 2,481
interactions. Of the 2,481 interactions, 1,536 of them are for SARS-CoV-2 and
were reported by 44 publications and preprints; the rest are related to
coronaviruses that cause Severe Acute Respiratory Syndrome (SARS) or Middle
Eastern Respiratory Syndrome (MERS).

  8. (8)

Homo sapiens PPI (with the shorthand H. sapiens) [60, 67, 50]: This network
has 24,407 nodes and 390,420 edges. Its largest connected component has 24,379
nodes and 390,397 edges. We use the complete network in our experiments. The
nodes represent proteins in the organism _Homo sapiens_ , and the edges encode
physical interactions between these proteins.

We now describe the eight synthetic networks that we examine:

  1. (9)

ER1 and ER2: An Erdős–Rényi (ER) network [1, 55], which we denote by
$\textup{ER}(n,p)$, is a random-graph model. The parameter $n$ is the number
of nodes and the parameter $p$ is the independent, homogeneous probability
that each pair of distinct nodes has an edge between them. The network ER1 is
an individual graph that we draw from $\textup{ER}({5000},0.01)$, and ER2 is
an individual graph that we draw from $\textup{ER}({5000},0.02)$.

  2. (10)

WS1 and WS2: A Watts–Strogatz (WS) network, which we denote by
$\textup{WS}(n,k,p)$, is a random-graph model to study the small-world
phenomenon [5, 55]. In the version of WS networks that we use, we start with
an $n$-node ring network in which each node is adjacent to its $k$ nearest
neighbors. With independent probability $p$, we then remove and rewire each
edge so that it connects a pair of distinct nodes that we choose uniformly at
random. The network WS1 is an individual graph that we draw from
$\textup{WS}({5000},50,0.05)$, and WS2 is an individual graph that we draw
from $\textup{WS}({5000},50,0.10)$.

  3. (11)

BA1 and BA2: A Barabási–Albert (BA) network, which we denote by
$\textup{BA}(n,{n_{0}})$, is a random-graph model with a linear preferential-
attachment mechanism [6, 55]. In the version of BA networks that we use, we
start with ${n_{0}}$ isolated nodes and we introduce new nodes with ${n_{0}}$
new edges each that attach preferentially (with a probability that is
proportional to node degree) to existing nodes until we obtain a network with
$n$ nodes. The network BA1 is an individual graph that we draw from
$\textup{BA}({5000},25)$, and BA2 is an individual graph that we draw from
$\textup{BA}({5000},50)$.

  4. (12)

SBM1 and SBM2: We use stochastic-block-model (SBM) networks in which each
block is an ER network [3]. Fix disjoint finite sets $C_{1}\cup\cdots\cup
C_{{k_{0}}}$ and a $k_{0}\times k_{0}$ matrix $B$ whose entries are real
numbers between $0$ and $1$. An SBM network, which we denote by
$\textup{SBM}(C_{1},\ldots,C_{{k_{0}}},B)$, has the node set
$V=C_{1}\cup\cdots\cup C_{{k_{0}}}$. For each node pair $(x,y)$, there is an
edge between $x$ and $y$ with independent probabilities $B[i_{0},j_{0}]$, with
indices $i_{0},j_{0}\in\\{1,\ldots,k_{0}\\}$ such that $x\in C_{i_{0}}$ and
$y\in C_{j_{0}}$. If $k_{0}=1$ and $B$ has a constant $p$ in all entries, this
SBM specializes to the Erdős–Rényi (ER) random-graph model $\textup{ER}(n,p)$
with $n=|C_{1}|$. The networks SBM1 and SBM2 are individual graphs that we
draw from $\textup{SBM}(C_{1},\ldots,C_{k_{0}},B)$ with
$|C_{1}|=|C_{2}|=|C_{3}|=\text{1,000}$, where $B$ is the $3\times 3$ matrix
whose diagonal entries are $0.5$ in both cases and whose off-diagonal entries
are $0.001$ for SBM1 and $0.1$ for SBM2. Both networks have 3,000 nodes; SBM1
has 752,450 edges and SBM2 has 1,049,365 edges.

### Types of noise

We now describe the three noise types that we consider in our network-
denoising experiments (see Figures 8 and 9). These noise types are as follows:

  1. (1)

(Noise type: ${-}\textup{ER}$) Given a network $G=(V,E)$, we choose a spanning
tree $T$ of $G$ (such a tree includes all nodes of $G$) uniformly at random
from all possible spanning trees. Let $E_{0}$ denote the set of all edges in
$G$ that are not in the edge set of $T$. We then obtain a corrupted network
$G^{\prime}$ by uniformly randomly removing half of the edges in $E_{0}$ from
$G$. Note that $G^{\prime}$ is guaranteed to be connected.

  2. (2)

(Noise type: ${+}\textup{ER}$) Given a network $G=(V,E)$, we uniformly
randomly choose a set $E_{2}$ of pairs of non-adjacent nodes in $G$ of size
$|E_{2}|=\lfloor|E|/2\rfloor$. The corrupted network is $G^{\prime}=(V,E\cup
E_{2})$, and we note that $50\%$ of the edges in $G^{\prime}$ are new.

  3. (3)

(Noise type: ${+}\textup{WS}$) Given a network $G=(V,E)$, fix integers
$n_{0}\in\\{1,\dots,|V|\\}$ and $k\in\\{1,\dots,n_{0}\\}$ and fix a real
number $p\in[0,1]$. We uniformly randomly choose a subset $V_{0}\subseteq V$
(with $|V_{0}|=n_{0}$) of the nodes of $G$. We generate a network
$H=(V_{0},E_{3})$ from the Watts–Strogatz model $\textup{WS}(n_{0},k_{0},p)$
using the node set $V_{0}$. We then obtain the corrupted network
$G^{\prime}=(V,E\cup E_{3})$, which has $|E_{3}|=n_{0}\lfloor k/2\rfloor$ new
edges. If $G$ is

Caltech, SNAP FB, arXiv, or Coronavirus, we use the parameters $n_{0}=100$,
$k=20$, and $p=0.3$. In this case, $G^{\prime}$ has 1,000 new edges. If $G$ is
H. sapiens, we use the parameters $n_{0}=500$, $k=120$, and $p=0.3$. In this
case, $G^{\prime}$ has 30,000 new edges.

LEARNING LOW-RANK LATENT MESOSCALE STRUCTURES IN NETWORKS

### Data availability

The datasets generated during the current study are available in the
repository https://github.com/HanbaekLyu/NDL_paper. For datasets of real-world
networks analysed during the current study, see the ‘Data sets’ section in our
Methods.

### Code availability

Our code for the main algorithms and simulations is publicly available in the
repository https://github.com/HanbaekLyu/NDL_paper. We also provide a user-
friendly version as a Python package ndlearn at
https://github.com/jvendrow/Network-Dictionary-Learning.

LEARNING LOW-RANK LATENT MESOSCALE STRUCTURES IN NETWORKS

Supplementary Information:

LEARNING LOW-RANK LATENT MESOSCALE STRUCTURES IN NETWORKS

##  Appendix A Overview

In this supplement, we present our algorithms for network dictionary learning
(NDL) and network denoising and reconstruction (NDR), and we prove theoretical
results about their convergence and error bounds. We give the complete NDL
algorithm (see Algorithm NDL) in Section D, and we give the complete NDR
algorithm (see Algorithm NDR) in Section E. We introduce the notion of ‘latent
motif dominance’ in Section D.2 to measure the significance of each latent
motif that we learn from a network. In Section G, we give a rigorous analysis
of the NDL and NDR algorithms.

##  Appendix B Problem Formulation for Network Dictionary Learning (NDL)

###  B.1. Definitions and notation

To facilitate our discussions, we use terminology and notation from [38, Ch.
3]. In the main manuscript, we described a network as a graph $G=(V,E)$ with a
node set $V$ and an edge set $E$ without directed or multi-edges, but possibly
with self-edges. One can characterize the edge set $E$ of $\mathcal{G}$ using
an adjacency matrix $A_{G}:V^{2}\rightarrow\\{0,1\\}$, where
$A(x,y)=\mathbb{1}(\\{x,y\\}\in E)$ for each $x,y\in V$. The function
$\mathbb{1}(S)$ denotes the indicator of the event $S$; it takes the value $1$
if $S$ occurs and takes the value $0$ if it does not occur. In this
supplementary information, we formulate our NDL framework in the more general
setting in which edges in a network may have weights. Although one can extend
the above definition of networks to include weighted edges by adjoining an
additional entry to $G=(V,E)$ for edge weights, it is convenient to extend the
range of adjacency matrices from $\\{0,1\\}$ to the interval $[0,\infty)$.

We define a network as a pair $\mathcal{G}=(V,A_{\mathcal{G}})$ with node set
$V$ and a weight matrix (which is also often called a ‘weighted adjacency
matrix’) $A_{\mathcal{G}}:V^{2}\rightarrow[0,\infty)$ that encodes the weights
of the contacts between nodes. For simplicity, we often drop the subscript
$\mathcal{G}$ in $A_{\mathcal{G}}$ and denote it as $A$. A graph $G=(V,E)$
determines a unique network $\mathcal{G}=(V,A_{G})$, where $A_{G}$ is the
adjacency matrix of $G$. The set $V(\mathcal{G})$ is the node set of the
network $\mathcal{G}$, which has size $|V(\mathcal{G})|$, where $|S|$ is the
number of elements in the set $S$. A pair $(x,y)$ of nodes in $\mathcal{G}$ is
called a directed edge if $A(x,y)>0$ and it is called a self-edge if $x=y$ and
$A(x,x)>0$. We say that a network $\mathcal{G}=(V,A)$ is symmetric if its
weight matrix is symmetric (i.e., $A(x,y)=A(y,x)$ for all $x,y\in V$) and we
say that it is binary (i.e., unweighted) if $A(x,y)\in\\{0,1\\}$ for all
$x,y\in V$. Given a symmetric network $\mathcal{G}=(V,A)$, an unordered pair
$\\{x,y\\}$ of nodes in $\mathcal{G}$ is an edge if $A(x,y)=A(y,x)>0$. The
network $\mathcal{G}$ is bipartite if it admits a ‘bipartition’, which is a
partition $V=V_{1}\cup V_{2}$ of the node set $V$ such that $V=V_{1}\cup
V_{2}$ and $A(x,y)=0$ if $x,y\in V_{1}$ or $x,y\in V_{2}$ for each $x,y\in V$.
If two networks $\mathcal{G}=(V,A)$ and
$\mathcal{G}^{\prime}=(V^{\prime},A^{\prime})$ satisfy $V^{\prime}\subseteq V$
and $A^{\prime}(x,y)\leq A(x,y)$ for all $x,y\in V^{\prime}$, then we say that
$\mathcal{G}^{\prime}$ is a subgraph of $\mathcal{G}$ and write
$\mathcal{G}^{\prime}\subseteq\mathcal{G}$. If $A^{\prime}(x,y)=A(x,y)$ for
all $x,y\in V^{\prime}$, then we say that $\mathcal{G}^{\prime}$ an induced
subgraph of $\mathcal{G}$ that is induced by the node set $V^{\prime}$.

For an integer $k\geq 2$ and nodes $x,y\in V$, we refer to a sequence
$(x_{1},\ldots,x_{k})$ of (not necessarily distinct) nodes in $\mathcal{G}$ as
a $k$-walk from $x$ to $y$ if $A(x_{i},x_{i+1})>0$ for all
$i\in\\{1,\ldots,k-1\\}$ and $(x_{1},x_{k})=(x,y)$. A $k$-walk
$(x_{1},\cdots,x_{k})$ is a $k$-path if all nodes $x_{1},\dots,x_{k}$ are
distinct. We say that the network $\mathcal{G}$ is connected if for any nodes
$x,y\in V$, there exists a $k$-path from $x$ to $y$ for some $k\geq 1$. If
$\mathcal{G}$ is connected, then for any two distinct nodes $x,y\in V$, we
define $d_{\mathcal{G}}(x,y)$ to be the smallest possible integer $k\geq 1$
such that there exists a $k$-walk from $x$ to $y$. The quantity
$d_{\mathcal{G}}(x,y)$ is called the shortest-path distance between $x$ and
$y$. The maximum of $d_{\mathcal{G}}(x,y)$ over all node pairs $(x,y)$ is the
shortest-path diameter $\textup{diam}(\mathcal{G})$ of $\mathcal{G}$. It
equals the minimum number of edges in a walk in $\mathcal{G}$ that connects
$x$ and $y$.

Suppose that we are given $N$ elements $\mathbf{v}_{1},\ldots,\mathbf{v}_{m}$
in some vector space. When we say that we take their mean, we refer to their
sample mean $\bar{\mathbf{v}}=N^{-1}\sum_{i=1}^{N}\mathbf{v}_{i}$. When we say
that we take a weighted average of them, we refer to the expectation
$\sum_{i=1}^{N}\mathbf{v}_{i}p_{i}$, where $(p_{1},\ldots,p_{N})$ is a
probability distribution on the set of $N$ elements.

###  B.2. Homomorphisms between networks and motif sampling

The ability to sample from a complex data set according to a known probability
distribution (e.g., a uniform one) is a crucial ingredient in dictionary-
learning problems. For instance, in image-processing applications [15, 23,
25], it is straightforward to uniformly randomly sample a $k\times k$ patch
from an image. However, it is not straightforward to uniformly randomly sample
a connected $k$-node subgraph of a network [68, 13, 17, 16]. For our purpose
of developing dictionary learning for networks, we use motif sampling, which
was introduced recently in [59]. In motif sampling, instead of directly
sampling a connected subgraph, one samples a random function from the node set
of a smaller network (i.e., a motif) into the node set of a target network
that preserves adjacency relations, and one then uses the subgraph that is
induced by the nodes in the image of the function. As we discuss below, such a
function between networks is a homomorphism.

Fix an integer $k\geq 1$ and a weight matrix
$A_{F}:[k]^{2}\rightarrow[0,\infty)$, where we use the shorthand notation
$[k]=\\{1,\ldots,k\\}$. We use the term motif for the corresponding network
$F=([k],A_{F})$. A motif is a network, and we use motifs to sample from a
given (and much larger) network. The type of motif that particularly interests
us is a $k$-chain, for which
$A_{F}=\mathbb{1}(\\{(1,2),(2,3),\ldots,(k-1,k)\\})$. A $k$-chain is a
directed path with node set $[k]$. For a general $k$-node motif
$F=([k],A_{F})$ and a network $\mathcal{G}=(V,A)$, we define the probability
distribution $\pi_{F\rightarrow\mathcal{G}}$ on the set $V^{[k]}$ of all node
maps (i.e., functions between node sets) $\mathbf{x}:[k]\rightarrow V$ by

$\pi_{F\rightarrow\mathcal{G}}(\mathbf{x}):=\frac{1}{\mathtt{Z}}\left(\prod_{i,j\in\\{1,\ldots,k\\}}A(\mathbf{x}(i),\mathbf{x}(j))^{A_{F}(i,j)}\right)\,,$
(1)

where $\mathtt{Z}=\mathtt{Z}(F,\mathcal{G})$ is a normalization constant that
we call the homomorphism density of $F$ in $\mathcal{G}$ [38]. A node map
$\mathbf{x}:[k]\rightarrow V$ a homomorphism $F\rightarrow\mathcal{G}$ if
$\pi_{F\rightarrow\mathcal{G}}(\mathbf{x})>0$, which is the case if and only
if $A(\mathbf{x}(a),\mathbf{x}(b))>0$ for all $a,b\in[k]$ with $A_{F}(a,b)>0$
(with the convention that ${\zeta}^{0}=1$ for all ${\zeta}\in\mathbb{R}$).
Informally, this means that if $(i,j)$ is a directed edge in $F$, then
$(\mathbf{x}(i),\mathbf{x}(j))$ is a directed edge in $\mathcal{G}$. Motif
sampling refers to the problem of sampling a random homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$ according to the distribution (1).

For our goal of learning interpretable latent motifs, it is important to
sample a homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$ such that
$\mathbf{x}$ is injective. This ensures that the nodes
$\mathbf{x}(1),\ldots,\mathbf{x}(k)$ in $V$ that correspond to nodes
$1,\ldots,k$ in the motif $F$ through the homomorphism $\mathbf{x}$ are all
distinct. When $\mathbf{x}:F\rightarrow\mathcal{G}$ is an injective
homomorphism, we denote $\mathbf{x}:F\hookrightarrow\mathcal{G}$. Using the
subgraph of $\mathcal{G}$ that is induced by the node set
$\\{\mathbf{x}(1),\ldots,\mathbf{x}(k)\\}$ when $\mathbf{x}$ is injective
returns a $k$-node subgraph of $\mathcal{G}$. For convenience, we define the
probability distribution

$\pi_{F\hookrightarrow\mathcal{G}}(\mathbf{x}):=C\,\pi_{F\rightarrow\mathcal{G}}(\mathbf{x})\cdot\mathbb{1}(\textup{$\mathbf{x}(1),\ldots,\mathbf{x}(k)$
are distinct})\,,$ (2)

where $C>0$ is a normalization constant. Injective motif sampling refers to
the problem of sampling a random injective homomorphism
$\mathbf{x}:F\hookrightarrow\mathcal{G}$ according to the distribution (2).
The probability distribution (2) is well-defined as long as there exists an
injective homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$. For instance, if
$F$ is a $k$-chain motif and $k$ is larger than the shortest-path diameter of
$\mathcal{G}$, then there is no injective homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$ and the probability distribution (2) is
not well-defined. As a special case of interest, for a symmetric and binary
motif $F$ and a network $\mathcal{G}$, the distributions
$\pi_{F\rightarrow\mathcal{G}}$ and $\pi_{F\hookrightarrow\mathcal{G}}$ are
the uniform distribution among all homomorphisms $F\rightarrow\mathcal{G}$ and
among all injective homomorphisms $F\hookrightarrow\mathcal{G}$, respectively.
That is,

$\displaystyle\textup{$F$, $\mathcal{G}$ are symmetric and
binary}\quad\Longrightarrow\quad\begin{matrix}\pi_{F\rightarrow\mathcal{G}}=\textup{Uniform}(\\{\mathbf{x}:F\rightarrow\mathcal{G}\\})\\\
\pi_{F\hookrightarrow\mathcal{G}}=\textup{Uniform}(\\{\mathbf{x}:F\hookrightarrow\mathcal{G}\\})\end{matrix}\,,$
(3)

which is the case for all of our examples in the main manuscript. In Section
C, we discuss three Markov-chain Monte Carlo (MCMC) algorithms for motif
sampling and propose corresponding algorithms for injective motif sampling by
combining them with rejection sampling.

###  B.3. Mesoscale patches of networks

A homomorphism $F\rightarrow\mathcal{G}$ is a node map $V(F)\rightarrow
V(\mathcal{G})$ that maps the edges of $F$ to edges of $\mathcal{G}$, so it
maps $F$ onto a subgraph of $\mathcal{G}$. It thereby maps $F$ ‘into’
$\mathcal{G}$.  For each homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$
from a motif $F=([k],A_{F})$ into a network $\mathcal{G}=(V,A)$, we define a
$k\times k$ matrix

$\displaystyle
A_{\mathbf{x}}(a,b):=A\big{(}\mathbf{x}(a),\mathbf{x}(b)\big{)}\quad\text{for
all\,\, $a,b\in\\{1,\ldots,k\\}$}\,.$ (4)

We say that $A_{\mathbf{x}}$ in (4) is the mesoscale patch of $\mathcal{G}$
that is induced by the homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$,
which is specified uniquely by the homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$ and the weight matrix $A$. Given a
$k\times k$ matrix $B$ and a homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$, we say that the $(a,b)$ entries of $B$
are on-chain if $A_{F}(a,b)>0$ and off-chain otherwise. The condition
$A_{F}(a,b)>0$ implies that $A(\mathbf{x}(a),\mathbf{x}(b))>0$ by the
definition of the homomorphism $\mathbf{x}$, so the on-chain entries of
$A_{\mathbf{x}}$ are always positive (and are always $1$ if $\mathcal{G}$ is
unweighted). However, the off-chain entries of $A_{\mathbf{x}}$ are not
necessarily positive, so they encode meaningful information about a network
that one ‘detects’ with the homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$.
As an illustration, suppose that $F$ is the $6$-chain motif and that
$\mathcal{G}=(V,A)$ is an undirected and binary graph. For any homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$, we have

$\displaystyle A_{\mathbf{x}}=\begin{bmatrix}0&1&*&*&*&*\\\ 1&0&1&*&*&*\\\
*&1&0&1&*&*\\\ *&*&1&0&1&*\\\ *&*&*&1&0&1\\\ *&*&*&*&1&0\end{bmatrix}\,,$ (5)

where each entry $*$ of $A_{\mathbf{x}}$ is either $0$ or $1$. In this
example, the entries that we mark as $1$ are the on-chain entries of
$A_{\mathbf{x}}$ and the other entries are off-chain entries.

Let $\mathcal{G}_{\mathbf{x}}$ denote the induced subgraph of $\mathcal{G}$
whose node set is the image
$\textup{Im}(\mathbf{x}_{t})=\\{\mathbf{x}(a)\,|\,a\in\\{1,\ldots,k\\}\\}$ of
the homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$. If $\mathbf{x}$ has $k$
distinct nodes in its image, then the weight matrix of
$\mathcal{G}_{\mathbf{x}}$ is exactly the mesoscale patch $A_{\mathbf{x}}$.
However, this is not the case when $\mathbf{x}$ has fewer than $k$ distinct
nodes. In that situation, we cannot interpret the mesoscale patch
$A_{\mathbf{x}}$ as the weight matrix of the induced subgraph
$\mathcal{G}_{\mathbf{x}}$ of $\mathcal{G}$. (For example, see Figure 11 in
the main manuscript.) This motivates us to sample an injective homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$ according to the distribution (2),
instead of according to the distribution (1).

###  B.4. Problem formulation for network dictionary learning (NDL)

The goal of the NDL problem is to learn, for a fixed integer $r\geq 1$, a set
of $r$ nonnegative matrices $\mathcal{L}_{1},\ldots,\mathcal{L}_{r}$, of size
$k\times k$ and Frobenius norms of at most $1$, such that

$\displaystyle A_{\mathbf{x}}\approx
a_{1}(\mathbf{x})\mathcal{L}_{1}+\cdots+a_{r}(\mathbf{x})\mathcal{L}_{r}$ (6)

for each injective homomorphism $\mathbf{x}:F\hookrightarrow\mathcal{G}$ for
some coefficients $a_{1}(\mathbf{x}),\ldots,a_{r}(\mathbf{x})\geq 0$. For each
homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$, this implies that one can
approximate the mesoscale patch $A_{\mathbf{x}}$ of $\mathcal{G}$ that is
induced by $\mathbf{x}$ as a linear combination of the $r$ matrices
$\mathcal{L}_{1},\ldots,\mathcal{L}_{r}$. We say that the tuple
$(\mathcal{L}_{1},\ldots,\mathcal{L}_{r})$ is a network dictionary for
$\mathcal{G}$, and we say that each $\mathcal{L}_{i}$ is a latent motif of
$\mathcal{G}$. We identify a network dictionary
$(\mathcal{L}_{1},\ldots,\mathcal{L}_{r})$ with the nonnegative matrix
$W\in\mathbb{R}_{\geq 0}^{k^{2}\times r}$ whose $j^{\textup{th}}$ column is
the vectorization of the $j^{\textup{th}}$ latent motif $\mathcal{L}_{j}$ for
$j\in\\{1,\ldots,r\\}$. The choice of vectorization $\mathbb{R}^{k\times
k}\rightarrow\mathbb{R}^{k^{2}}$ is arbitrary, but we use a column-wise
vectorization in Algorithm A4. One can interpret each $\mathcal{L}_{i}$ as the
$k$-node weighted network with node set $\\{1,\ldots,k\\}$ and weight matrix
$\mathcal{L}_{i}$. See Figure 2 in the main manuscript for plots of latent
motifs as weighted networks.

For the latent motifs $\mathcal{L}_{i}$ to be interpretable as subgraphs of
$\mathcal{G}$, it is crucial that we require both their entries and the
coefficients $a_{i}(\mathbf{x})$ to be nonnegative. The nonnegativity
constraint on each latent motif $\mathcal{L}_{i}$ allows one to interpret each
$\mathcal{L}_{i}$ as the weight matrix of a $k$-node network. Additionally,
because the coefficients $a_{j}(\mathbf{x})$ are also nonnegative, the
approximate decomposition (6) implies that
$a_{i}(\mathbf{x})\mathcal{L}_{i}\lessapprox A_{\mathbf{x}}$. Therefore, if
$a_{i}(\mathbf{x})>0$, any network structure (e.g., large-degree nodes,
communities, and so on) in the latent motif $\mathcal{L}_{i}$ must also exist
in $A_{\mathbf{x}}$. Therefore, one can consider the latent motifs as
approximate $k$-node subgraphs $\mathcal{G}$ that exhibit ‘typical’ network
structures of $\mathcal{G}$ at scale $k$. In the spirit of Lee and Seung [7],
one can view the latent motifs as ‘parts’111Lee and Seung [7] discussed a
similar nonnegative decomposition in which the $A_{\mathbf{x}}$ are images of
faces. In that scenario, the learned factors capture parts of human faces
(such as eyes, noses, and mouths). of a network $\mathcal{G}$.

As a more precise formulation of (6), consider the stochastic optimization
problem

$\displaystyle\operatorname*{arg\,min}_{\begin{subarray}{c}\mathcal{L}_{1},\ldots,\mathcal{L}_{r}\in\mathbb{R}_{\geq
0}^{k\times k}\\\
{\lVert\mathcal{L}_{1}\rVert_{F},\ldots,\lVert\mathcal{L}_{r}\rVert_{F}\leq
1}\end{subarray}}\mathbb{E}_{\mathbf{x}\sim\pi_{F\hookrightarrow\mathcal{G}}}\left[\inf_{a_{1}(\mathbf{x}),\ldots,a_{r}(\mathbf{x})\geq
0}\left\lVert
A_{\mathbf{x}}-\sum_{i=1}^{r}a_{i}(\mathbf{x})\mathcal{L}_{i}\right\rVert_{F}\right]\,,$
(7)

where $\pi_{F\hookrightarrow\mathcal{G}}$ is the probability distribution that
we defined in (2) and $\lVert\cdot\rVert_{F}$ denotes the matrix Frobenius
norm.  The choice of the probability distribution
$\pi_{F\hookrightarrow\mathcal{G}}$ for the injective homomorphisms
$\mathbf{x}:F\hookrightarrow\mathcal{G}$ is natural because it becomes the
uniform distribution on the set of all injective homomorphisms
$F\hookrightarrow\mathcal{G}$ when the adjacency matrices for $\mathcal{G}$
and $F$ are both unweighted. Exactly solving the NDL optimization problem (7)
is computationally difficult because the objective function in it is non-
convex and it is not obvious how to sample an injective homomorphism
$F\hookrightarrow\mathcal{G}$ according to the distribution
$\pi_{F\hookrightarrow\mathcal{G}}$ that we defined in (2). In Section D, we
provide an algorithm for NDL that approximately solves (7).

###  B.5. Overview of our algorithms and their theoretical guarantees

We overview our algorithms and their theoretical guarantees. Our main
theoretical results (which are all novel) are Theorems G.4 and G.7 for NDL and
Theorems G.10 and G.14 for NDR. We summarize our algorithms and main results
in Table 2, and we compare and contrast them with the results in [64].

NDL | Sampling | Convergence | Efficient MCMC
---|---|---|---
Lyu et al. [64] |  $k$-walks | Non-bipartite networks | ✗
$\begin{matrix}\textup{The present}\\\ \textup{work}\end{matrix}$ | $\begin{matrix}\textup{$k$-paths}\\\ \textup{(Alg.\leavevmode\nobreak\ \ref{alg:motif_inj}) }\end{matrix}$ |  $\begin{matrix}\textup{Non-bipartite networks}\\\ \textup{(Thm.\leavevmode\nobreak\ \ref{thm:NDL}) }\end{matrix}$ $\begin{matrix}\textup{Bipartite networks}\\\ \textup{(Thm.\leavevmode\nobreak\ \ref{thm:NDL2}) }\end{matrix}$ | $\begin{matrix}\checkmark\\\ \textup{(Prop.\leavevmode\nobreak\ \ref{prop:approximate_pivot}) }\end{matrix}$
|  |  |

NDR | Sampling | denoising | Convergence | Error bound
---|---|---|---|---
Lyu et al. [64] |  $k$-walks | F | ✗ | ✗
$\begin{matrix}\textup{The present work}\end{matrix}$ | $\begin{matrix}\textup{$k$-walks}\\\ \textup{$k$-paths}\end{matrix}$ |  F, T | $\begin{matrix}\checkmark\\\ \textup{ (Thm.\leavevmode\nobreak\ \ref{thm:NR}{(i)}{--}{(ii)}) }\\\ \textup{ (Thm.\leavevmode\nobreak\ \ref{thm:NR2}{(i)}{--}{(ii)}) }\end{matrix}$ | $\begin{matrix}\checkmark\\\ \textup{ (Thm.\leavevmode\nobreak\ \ref{thm:NR}{(iii)}) }\\\ \textup{ (Thm.\leavevmode\nobreak\ \ref{thm:NR2}{(iii)}) }\end{matrix}$
|  |  |  |

Table 2.  A comparison of the algorithms and theoretical contributions of the
present work with those in [64]. In the table for NDR, denoising refers to the
Boolean variable in our NDR algorithm (see Algorithm NDR). The special case of
the NDR algorithm with $k$-walk sampling (i.e., $\texttt{InjHom}=\texttt{F}$)
and $\texttt{denoising}=\texttt{F}$ corresponds to the network-reconstruction
algorithm in [64].

Algorithm NDL:

Given a network $\mathcal{G}$, the NDL algorithm (see Algorithm NDL) computes
a sequence $(W_{t})_{t\geq 0}$ of network dictionaries (which take the form of
$k^{2}\times r$ matrices) of latent motifs.

Algorithm NDR:

Given a network $\mathcal{G}$, a network dictionary $W$, and a threshold
parameter $\theta>0$, the NDR algorithm (see Algorithm NDR) computes a
sequence of weighted networks $\mathcal{G}_{\textup{recons}}$ and binary
(i.e., unweighted) networks $\mathcal{G}_{\textup{recons}}(\theta)$.

Theorem G.4:

Given a non-bipartite network $\mathcal{G}$ and a choice of the parameters in
Algorithm NDL, we prove that the sequence $(W_{t})_{t\geq 0}$ of network
dictionaries converges almost surely to the set of stationary points of the
associated objective function in (7).

Theorem G.7:

Given a bipartite network $\mathcal{G}$ and a choice of the parameters in
Algorithm NDL, we prove a convergence result that is analogous to the one in
Theorem G.4.

Theorem G.10:

Given a non-bipartite target network $\mathcal{G}$ and a network dictionary
$W$, we show that (i) the sequence of weighted reconstructed networks
$\mathcal{G}_{\textup{recons}}$ that we obtain using the NDR algorithm (see
Algorithm NDR) converges almost surely to some limiting network, and (ii) we
obtain a closed-form expression of this limiting network.  We also show that
(iii) a suitable Jaccard reconstruction error between the original network
$\mathcal{G}$ and the limiting reconstructed network satisfies

$\displaystyle\textup{Jaccard reconstruction error}\leq\frac{\textup{mesoscale
approximation error}}{2(k-1)}\,,$ (8)

where $k$ denotes the mesoscale parameter (i.e., the number of nodes in a
$k$-chain motif) and the mesoscale approximation error is the mean $L_{1}$
distance between the $k\times k$ mesoscale patches of $\mathcal{G}$ and their
nonnegative linear approximations from the latent motifs in $W$.

Theorem G.14:

We show a convergence result that is analogous to the one in Theorem G.14 for
a bipartite target network $\mathcal{G}$.

##  Appendix C Markov-Chain Monte Carlo (MCMC) Motif-Sampling Algorithms

We mentioned in Section B.4 that one of the main difficulties in solving the
optimization problem (7) is to directly sample an injective homomorphism
$\mathbf{x}:F\hookrightarrow\mathcal{G}$ from the distribution
$\pi_{F\hookrightarrow\mathcal{G}}$ (see (2)). To overcome this difficulty, we
use the Markov-chain Monte Carlo (MCMC) algorithms that were introduced in
[59]. Although the algorithms in [59] apply to networks with edge weights
and/or node weights, we only use the simplified forms of them that we give in
Algorithms MP and MG. Algorithm MP with the option
$\texttt{AcceptProb}=\texttt{Approximate}$ is a novel algorithm of the present
paper. Using these MCMC sampling algorithms, we generate a sequence
$(\mathbf{x}_{t})_{t\geq 0}$ of homomorphisms $F\rightarrow\mathcal{G}$ such
that the distribution of $\mathbf{x}_{t}$ converges to
$\pi_{F\rightarrow\mathcal{G}}$ under some mild conditions on $\mathcal{G}$
and $F$ [64, Thm. 5.7].

Once we have an iterative motif-sampling algorithm, we combine it with a
standard rejection-sampling algorithm to shrink the support of the probability
distribution $\pi_{F\rightarrow\mathcal{G}}$ to that of injective
homomorphisms $\mathbf{x}:F\hookrightarrow\mathcal{G}$. (See, e.g., [12] for
background information about rejection sampling.) The key idea is to ignore
(i.e., ‘reject’) the unwanted instances in the trajectory
$(\mathbf{x}_{t})_{t\geq 0}$. In our case, the instances that we reject are
the homomorphisms $\mathbf{x}_{t}$ that are not injective. That is, we reject
situations in which $\mathbf{x}_{t}(1),\ldots,\mathbf{x}_{t}(k)$ are not all
distinct. We state our algorithm for injective motif sampling in Algorithm RM.

1:Input: Network $\mathcal{G}=(V,A)$, motif $F=([k],A_{F})$, and homomorphism
$\mathbf{x}:F\rightarrow\mathcal{G}$

2:While: $\mathbf{x}:F\rightarrow\mathcal{G}$ is injective (i.e.,
$\mathbf{x}^{\prime}(1),\ldots,\mathbf{x}^{\prime}(k)$ are distinct)

3: Update $\mathbf{x}$ to a new homomorphism $F\rightarrow\mathcal{G}$ using
either Algorithm MP or Algorithm MG

4:Output: Injective homomorphism $\mathbf{x}:F\hookrightarrow\mathcal{G}$

Algorithm RM . Injective motif sampling

Algorithm RM restricts the state space of the MCMC motif-sampling algorithms
(see Algorithms MG and MP) onto the subset of injective homomorphisms
$F\hookrightarrow\mathcal{G}$. By the strong Markov property, this restriction
is a Markov chain. Therefore, Algorithm RM is an MCMC algorithm for the
injective motif-sampling problem. (See Proposition G.3 for details.) When
there are only a few injective homomorphisms $F\hookrightarrow\mathcal{G}$
relative to the number of homomorphisms $F\rightarrow\mathcal{G}$, the
rejection step (i.e., the while loop) in Algorithm RM may take a while to
terminate. (It is inversely proportional to the probability that a random
homomorphism under the probability distribution
$\pi_{F\rightarrow\mathcal{G}}$ is injective.) For example, this is the case
when $\mathcal{G}$ is the network Coronavirus PPI and $F$ is a $k$-chain motif
with $k\geq 21$.

We now give more details about the MCMC algorithms that we employ for non-
injective motif sampling. In the pivot chain (see Algorithm MP with
AcceptProb=Exact), for each update $\mathbf{x}_{t}\mapsto\mathbf{x}_{t+1}$,
the pivot $\mathbf{x}_{t}(1)$ first performs a random-walk move on
$\mathcal{G}$ (see (9)) to move to a new node $\mathbf{x}_{t+1}(1)\in V$. It
accepts this move with a suitable acceptance probability (see (10)) according
to the Metropolis–Hastings algorithm (see, e.g., [52, Sec. 3.2]), so the
stationary distribution is exactly the target distribution. After the move
$\mathbf{x}_{t}(1)\mapsto\mathbf{x}_{t+1}(1)$, we sample each
$\mathbf{x}_{t+1}(i)\in V$ for $i\in\\{2,3,\ldots,k\\}$ successively from the
conditional distribution (11). This ensures that the desired distribution
$\pi_{F\rightarrow\mathcal{G}}$ in (1) is a stationary distribution of the
resulting Markov chain. In the Glauber chain, we select one node $i\in[k]$ of
$F$ uniformly at random, and we resample its location $\mathbf{x}_{t}(i)\in
V(\mathcal{G})$ at time $t$ to $\ell=\mathbf{x}_{t+1}(i)\in V$ from the
conditional distribution (12) (see Figure 12a). See [52, Sec. 3.3] for
discussions of the Metropolis–Hastings algorithm and Glauber-chain MCMC
sampling.

1:Input: Symmetric network $\mathcal{G}=(V,A)$, a $k$-chain motif
$F=([k],A_{F})$, and homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$

2:Parameters:
$\mathtt{AcceptProb}\in\\{\mathtt{Exact},\mathtt{Approximate}\\}$

3:Do: $\mathbf{x}^{\prime}\leftarrow\mathbf{x}$

4: If $\sum_{c\in V}A(\mathbf{x}(1),c)=0$: Terminate

5: Else:

6: Sample $\ell\in V$ at random from the distribution  $\displaystyle
p_{1}(w)=\frac{A(\mathbf{x}(1),w)}{\sum_{c\in V}A(\mathbf{x}(1),c)}\,\,,\quad
w\in V$ (9)

7: Compute the acceptance probability $\alpha\in[0,1]$ by
$\displaystyle\alpha\leftarrow\begin{cases}\min\left(\frac{\sum_{c\in[n]}A^{k-1}(\ell,c)}{\sum_{c\in[n]}A^{k-1}(\mathbf{x}(1),c)}\frac{\sum_{c\in
V}A(c,\mathbf{x}(1))}{\sum_{c\in
V}A(\mathbf{x}(1),c)},\,1\right)\,,&\text{if\quad$\mathtt{AcceptProb}=\mathtt{Exact}$}\\\\[10.0pt]
\min\left(\frac{\sum_{c\in V}A(c,\mathbf{x}(1))}{\sum_{c\in
V}A(\mathbf{x}(1),c)},\,1\right)\,,&\text{if\quad$\mathtt{AcceptProb}=\mathtt{Approximate}$}\end{cases}$
(10)

8: Sample $U\in[0,1]$ uniformly at random, independently of everything else

9: $\ell\leftarrow\mathbf{x}(1)$ if $U>\lambda$ and
$\mathbf{x}^{\prime}(1)\leftarrow\ell$

10: For $i=2,3,\ldots,k$:

11: Sample $\mathbf{x}^{\prime}(i)\in V$ from the distribution  $\displaystyle
p_{i}(w)=\frac{A(\mathbf{x}(i-1),w)}{\sum_{c\in
V}A(\mathbf{x}(i-1),c)}\,\,,\quad w\in V$ (11)

12:Output: Homomorphism $\mathbf{x}^{\prime}:F\rightarrow\mathcal{G}$

Algorithm MP . Pivot-Chain Update

1:Input: Network $\mathcal{G}=(V,A)$, a $k$-chain motif $F=([k],A_{F})$, and
homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$

2:Do: Sample $v\in[k]$ uniformly at random

3: Sample $z\in V$ at random from the distribution  $\displaystyle
p(w)=\frac{1}{Z}\left(\prod_{u\in[k]}A(\mathbf{x}(u),w)^{A_{F}(u,v)}\right)\left(\prod_{u\in[k]}A(w,\mathbf{x}(u))^{A_{F}(v,u)}\right)\,,\quad
w\in V$ (12) where $Z=\sum_{c\in
V}\left(\prod_{u\in[k]}A(\mathbf{x}(u),c)^{A_{F}(u,v)}\right)\left(\prod_{u\in[k]}A(c,\mathbf{x}(u))^{A_{F}(v,u)}\right)$
is the normalization constant

4: Define a new homomorphism $\mathbf{x}^{\prime}:F\rightarrow\mathcal{G}$ by
$\mathbf{x}^{\prime}(w)=z$ if $w=v$ and $\mathbf{x}^{\prime}(w)=\mathbf{x}(w)$
otherwise

5:Output: Homomorphism $\mathbf{x}^{\prime}:F\rightarrow\mathcal{G}$

Algorithm MG . Glauber-Chain Update

Let $\Delta$ denote the maximum degree (i.e., number of neighbors) of the
nodes of the network $\mathcal{G}=(V,A)$. We also say that the network
$\mathcal{G}$ itself has a maximum degree of $\Delta$. The Glauber chain has
an efficient local update (with a computational complexity of $O(\Delta)$). It
converges quickly to the stationary distribution
$\pi_{F\rightarrow\mathcal{G}}$ only for networks that are dense enough so
that two homomorphisms that differ at one node have a probability of at least
$1/(2\Delta)$ to coincide after a single Glauber-chain update. See [59, Thm.
6.1] for a precise statement of this fact.

The pivot chain (see Algorithm MP with $\texttt{AcceptProb}=\texttt{Exact}$)
has more computationally expensive local updates than the Glauber chain. The
pivot chain has a computational complexity of $O(\Delta^{k-1})$ (as discussed
in [59, Remark 5.6]), but it converges as fast as a ‘lazy’ random walk on a
network. (In a lazy random walk, each move has a chance to be rejected; see
[59, Thm. 6.2].) In our computational experiments, we find that the Glauber
chain is slow, especially for sparse networks (e.g., for COVID PPI, which has
an edge density of $0.0010$, and UCLA, which has an edge density of $0.0037$)
and that the pivot chain is too expensive to compute for chain motifs with
$k\geq 21$. As a compromise, to simultaneously have low computational
complexity and fast convergence (which is as fast as the standard random
walk), we employ an approximate pivot chain, which is Algorithm MP with the
option $\texttt{AcceptProb}=\texttt{Approximate}$. Specifically, we compute
the acceptance probability $\alpha$ in (10) only approximately and thereby
reduce the computational cost to $O(\Delta)$. Our compromise, which we discuss
in the next paragraph, is that the stationary distribution of the approximate
pivot chain may be slightly different from our target distribution
$\pi_{F\rightarrow\mathcal{G}}$.

Define a probability distribution $\hat{\pi}_{F\rightarrow\mathcal{G}}$ on the
set of all node maps $\mathbf{x}:[k]\rightarrow V$ as

$\displaystyle\hat{\pi}_{F\rightarrow\mathcal{G}}(\mathbf{x}):=\frac{\prod_{i=1}^{k}A(\mathbf{x}(i-1),\mathbf{x}(i))}{|V|\sum_{y_{2},\ldots,y_{k}\in
V}A(\mathbf{x}(1),y_{2})\prod_{i=3}^{k}A(y_{i-1},y_{i})}\,.$ (13)

According to Proposition G.2, the stationary distribution of the approximate
pivot chain is (13). The distribution (13) is different from the desired
target distribution $\pi_{F\rightarrow\mathcal{G}}$. Specifically,
$\pi_{F\rightarrow\mathcal{G}}(\mathbf{x})$ is proportional only to the
numerator in (13) and the sum in the denominator in (13) is a weighted count
of the homomorphisms $\mathbf{y}:F\rightarrow\mathcal{G}$ for which
$\mathbf{y}(1)=\mathbf{x}(1)$. Therefore, under
$\hat{\pi}_{F\rightarrow\mathcal{G}}$, we penalize the probability of each
homomorphism $\mathbf{x}:F\rightarrow\mathcal{G}$ according to the number of
$k$-step walks in $\mathcal{G}$ that start from $\mathbf{x}(1)\in V$. (The
exact acceptance probability in (10) neutralizes this penalty.) It follows
that $\hat{\pi}_{F\rightarrow\mathcal{G}}$ is close to
$\pi_{F\rightarrow\mathcal{G}}$ when the $k$-step-walk counts that start from
each node in $\mathcal{G}$ do not differ too much for different nodes. For
example, on degree-regular networks like lattices, such counts do not depend
on the starting node, and it thus follows that
$\hat{\pi}_{F\rightarrow\mathcal{G}}=\pi_{F\rightarrow\mathcal{G}}$.
Nevertheless, despite the potential discrepancy between
$\pi_{F\rightarrow\mathcal{G}}$ and $\hat{\pi}_{F\rightarrow\mathcal{G}}$, the
approximate pivot chain gives good results for the reconstruction and
denoising experiments that we showed in Figures 5 and 9 in the main
manuscript.

##  Appendix D Algorithm for Network Dictionary Learning (NDL)

###  D.1. Algorithm overview and statement

The essential idea behind our algorithm for NDL (see Algorithm NDL) is as
follows. Suppose that we compute all possible injective homomorphisms
$\mathbf{x}_{1},\dots,\mathbf{x}_{M}:F\hookrightarrow\mathcal{G}$ and their
corresponding mesoscale patches $A_{\mathbf{x}_{t}}$ for
$t\in\\{1,\ldots,M\\}$. These $M$ mesoscale patches of $\mathcal{G}$ form the
data set (a so-called ‘batch’) in which we apply a dictionary-learning
algorithm. To do this, we column-wise vectorize each of these $k\times k$
matrices (using Algorithm A4) and obtain a $k^{2}\times M$ data matrix $X$,
and we then apply nonnegative matrix factorization (NMF) [7] to obtain a
$k^{2}\times r$ nonnegative matrix $W$ for some fixed integer $r\geq 1$ to
yield an approximate factorization $X\approx WH$ for some nonnegative matrix
$H$. From this procedure, we approximate each column of $X$ by the nonnegative
linear combination of the $r$ columns of $W$ with coefficients that are given