Blind identification of stochastic block models from dynamical observations

TL;DR

This paper introduces spectral algorithms for blind identification of stochastic block models from dynamical observations, enabling recovery of network structure and parameters without direct edge knowledge.

Contribution

It presents novel spectral methods with provable guarantees for inferring SBM partitions and parameters solely from nodal diffusive process observations.

Findings

Spectral algorithms achieve high-accuracy partition recovery.

The methods are supported by theoretical guarantees based on random matrix theory.

Numerical experiments validate the effectiveness of the proposed algorithms.

Abstract

We consider a blind identification problem in which we aim to recover a statistical model of a network without knowledge of the network's edges, but based solely on nodal observations of a certain process. More concretely, we focus on observations that consist of single snapshots taken from multiple trajectories of a diffusive process that evolves over the unknown network. We model the network as generated from an independent draw from a latent stochastic block model (SBM), and our goal is to infer both the partition of the nodes into blocks, as well as the parameters of this SBM. We discuss some non-identifiability issues related to this problem and present simple spectral algorithms that provably solve the partition recovery and parameter estimation problems with high accuracy. Our analysis relies on recent results in random matrix theory and covariance estimation, and associated concentration inequalities. We illustrate our results with several numerical experiments.

Equations

$${\mathbf{x}}_{t+1}^{(i)}=f(\mathbf{A}){\mathbf{x}}_t^{(i)},i=1,\dots ,s,$$

$$A_{ij}\sim \text{Ber}({\Omega }_{g_i,g_j}).$$

$$\mathcal{A}:=\mathbb{E}[\mathbf{A}|\mathbf{G}]=\mathbf{G}\mathbf{\Omega }{\mathbf{G}}^{\top }.$$

$$\mathcal{L}:=\mathbb{E}[\mathbf{D}{]}^{-1/2}\mathbb{E}[\mathbf{A}]\mathbb{E}[\mathbf{D}{]}^{-1/2}=\text{diag}(\mathcal{A}1{)}^{-1/2}\mathcal{A}\text{ diag}(\mathcal{A}1{)}^{-1/2}.$$

$$\mathbf{V}=\mathbf{G}({\mathbf{G}}^{\top }\mathbf{G}{)}^{-1/2}\mathcal{U}.$$

$$\tilde{\mathbf{V}}=\mathbf{G}\mathcal{U}.$$

$${\delta }_{\text{min}}>27\ln \left(\frac{4n}{ϵ}\right).$$

$$\Vert \mathbf{L}-\mathcal{L}\Vert \le 3\sqrt{\frac{3\ln (4n/ϵ)}{{\delta }_{\text{min}}}},$$

$$f:\mathbf{A}\mapsto \sum _{\ell =1}^R{\alpha }_{\ell }{\mathbf{L}}^{\ell },$$

$${\mathbf{x}}_{t+1}=\mathbf{L}{\mathbf{x}}_t,$$

$${\mathbf{x}}_{t+1}=f(\mathbf{A}){\mathbf{x}}_t,$$

$$\hat{\mathbf{\Sigma }}(T):=\frac{1}{s}\sum _{i=1}^s[{\mathbf{x}}_T^{(i)}-{\bar{\mathbf{x}}}_T][{\mathbf{x}}_T^{(i)}-{\bar{\mathbf{x}}}_T{]}^{\top },$$

$$\mathbf{S}(T):=f(\mathcal{A}{)}^{2T}.$$

$$\Vert \hat{\mathbf{\Sigma }}(T)-\mathbf{S}(T)\Vert \le 2TM+B(T),$$

$$M:=3\sqrt{\frac{3\ln (4n/ϵ)}{{\delta }_{\text{min}}}},B(T):=C\left(1+\sum _{i=2}^k\min \{1,|{\mu }_i|+M{\}}^{2T}+(n-k)M^{2T}\right)\sqrt{\frac{\ln s}{s}},$$

$$\Vert \hat{\mathbf{\Sigma }}(T)-\mathbf{S}(T)\Vert \le 2T{M}_f+B_f(T),$$

$$M_f$$

$$B_f(T)$$

$$q\le \frac{64k}{n{\tau }^2{\xi }_k^2}(2T{M}_f+B_f(T){)}^2,$$

$$\mathbf{\Sigma }(T)$$

$$[\mathcal{L}{]}_{ij}=\{\begin{array}{ll}\frac{a}{n}\cdot \frac{k}{a+(k-1)b},& \text{if }{g}_i=g_j,\\ \frac{b}{n}\cdot \frac{k}{a+(k-1)b},& \text{if }{g}_i\ne g_j,\end{array}$$

$$[{\mathcal{L}}^{2T}{]}_{ij}=\{\begin{array}{ll}\frac{k-1}{n}{\left(\frac{a-b}{a+(k-1)b}\right)}^{2T}+\frac{1}{n}& \text{if }{g}_i=g_j,\\ -\frac{1}{n}{\left(\frac{a-b}{a+(k-1)b}\right)}^{2T}+\frac{1}{n}& \text{if }{g}_i\ne g_j.\end{array}$$

$${\lambda }_i(\mathbf{S}(T))=\{\begin{array}{ll}1,& \text{if }i=1,\\ {\left(\frac{a-b}{a+(k-1)b}\right)}^{2T},& \text{if }i=2,\dots ,k.\end{array}$$

$$h(\mathbf{\Omega })=\rho -\frac{1^{\top }\mathcal{A}1}{n^2}=0.$$

$$\mathbf{S}(T)$$

$$\lambda (\mathbf{S}(T))$$

$$h(\mathbf{\Omega })$$

$$\hat{\mathbf{\Sigma }}(T)$$

$$\lambda (\hat{\mathbf{\Sigma }}(T))$$

$$h(\mathbf{\Omega })$$

$$({\lambda }_2({\mathcal{L}}^{2T}){)}^{1/2T}=\frac{a-b}{a+(k-1)b}.$$

$$a=({\lambda }_2({\mathcal{L}}^{2T}){)}^{1/2T}\rho n+\rho n.$$

$$\hat{a}=({\lambda }_2(\hat{\mathbf{\Sigma }}(T)){)}^{1/2T}\rho n+\rho n.$$

Full text

{AMS}

68R10, 62M15, 60B20, 15A29, 15A18

1 Introduction

We consider the problem of inferring a statistical model of a network with unknown adjacency matrix $\bm{A}\in\{0,1\}^{n\times n}$, based on snapshots of multiple trajectories of a dynamical process that evolves over the network. Specifically, we are interested in the following setup (a more detailed description will be provided in Section 3): We observe $s$ independent samples $\bm{x}^{(1)}_{T},\ldots,\bm{x}^{(s)}_{T}$ of a dynamical process at a fixed time $T$. Each sample is generated by the discrete-time system

[TABLE]

initialized at some unknown $\bm{x}_{0}^{(i)}$, independently for each $i$. Each $\bm{x}_{0}^{(i)}$ is modeled as a random vector with unit covariance matrix ($\mathbb{E}[\bm{x}_{0}^{(i)}[\bm{x}_{0}^{(i)}]^{\top}]=\bm{I}$) drawn from some sub-Gaussian distribution (e.g., a standard multivariate Gaussian). Furthermore, we assume that the function $f$ describing the dynamical network process is known and corresponds to a class of diffusive processes. Finally, $\bm{A}$ is assumed to be drawn independently according to a stochastic block model (SBM), so that $\bm{A}$ is the same, albeit unknown, for each sample $\bm{x}_{T}^{(i)}$. Based on such data, our task is to identify the generative model underpinning the network: we wish to infer the parameters and the latent node partition of the SBM, without observing the network edges, by relying solely on the observations $\{\bm{x}_{T}^{(i)}\}_{i=1}^{s}$ at the nodes.

1.1 Motivation

In many applications we are confronted with the following scenario: we observe snapshots of the state of a system at particular times, and based on these observations we wish to infer the nature of the (dynamical) interactions between the entities that we observe. Examples abound and can be found across different fields [54, 26, 27].

Consider, e.g., measuring some node activities such as the expression of opinions at $s$ different, sufficiently separated instances of time in a social network. While the underlying social connections are typically not observable, the latent (unobserved) network will be essentially static, and may be assumed to have emerged according to some homophilic interactions describable by an SBM. The measured opinion profiles will be shaped by both some (random) initial personal opinions of the individuals as well as certain network effects. Similar examples may arise when measuring (structurally stable) biological systems responding to different stimuli. For instance, we might have access to a set of measurements from certain regions of the brain and would like to infer the coupling between these regions. As another example, we might be interested in inferring a network of genetic interactions from multiple snapshots of gene expression, e.g., as obtained from single cell RNA-sequencing [43]. As taking these expression measurements will destroy the system under consideration, we are left with a few snapshots of data, from which we want to reconstruct (a model of) the underlying interactions. Finally, we may consider the monitoring of ecological populations from abundance data [13], where we often can only access snapshot information at specific instances of time, but are interested in knowing the network of interactions between the species.

The above are canonical examples of problems where the system of interest is composed of a large number of different entities and can be conveniently conceptualized as a network. In this abstraction, nodes represent the individual entities and the couplings between these nodes take place along the edges of the network. In many settings, we may have reason to believe that an observed dynamical process emerges from the interaction between nodes over such a network. While we may have some knowledge about the dynamics of each individual node, the edges of the network – which enable the interactions between nodes – are often unknown or only partially known. In such cases it is often desirable to learn the structure of the network based on observational data. This is a general class of network inference problems, which is a fundamental constituent of network analysis, but has received limited attention.

On the other hand, there are many situations where inferring exactly the couplings in the network e.g., identifying all edges, is out of reach, for a number of reasons. First, we may not have access to sufficiently many samples. Second, the data and the network itself may be subject to some fluctuations with associated uncertainties. Finally, we may only be able to partially observe the dynamical system, so that the full model is non-identifiable.

In such cases, it is then reasonable to aim at an approximate description of the system at hand, in terms of a statistical model that captures the coarse-grained features of the underlying network. This is precisely the type of problem addressed in this paper, for the special case where the dynamics are of a linear, diffusive type, and the network model to be inferred is an SBM.

1.2 Contributions

We present a rigorous study of blind model identification for random networks. Our main contributions are as follows:

First, we present general algorithmic ideas on how to solve a blind identification problem using spectral methods, and introduce simple algorithms to solve the concrete problems considered in this paper (see Sections 4 and 5, and Algorithms 1 and 2).

Second, we establish a concentration inequality about the sample covariance matrix that is calculated on the basis of the observed samples (see Theorems 4.1 and 4.2).

Finally, we discuss how this concentration result translates into explicit bounds for the recovery error of the latent partitions (see Theorem 4.3), and provide bounds on the estimation errors for the parameters of the SBM, for the specific example of a planted partition model (see Theorem 5.3).

1.3 Related work

Many different network inference notions have featured in the literature. For instance, the estimation of “functional” couplings based on statistical association measures such as correlation or mutual information, is an analysis often encountered in the context of biological systems [19, 36, 14, 52]. Another type of network inference is causal inference, as discussed in, e.g., [5, 37, 53], in which one aims to infer directed graphs that encode direct causal influences between the entities. In contrast, the type of network inference we consider here is what may be called “topological” inference: given a system of dynamical units, we want to infer their direct “physical” interactions. Some recent surveys include [54, 62, 30, 12]. More specifically, for an evolving dynamical process we might want to infer the underlying adjacency matrix. This may be viewed as a system identification problem, as considered in Control Theory [26, 27]. However, in the system identification literature, one can usually control the inputs to the system and observe the full system response. This is different from our setting, where there is no control, and we only observe a limited number of snapshots.

Problems of network identification that are somewhat similar to ours have received wide interest in the literature recently. In many cases, formulations based on the optimization of a sparsity-promoting cost function are used [22, 8, 16], or formulations based on sparse dictionary learning [65, 6]. A number of methods also employ spectral characteristics of the observed outputs to identify (aspects of) the network topology [17, 31, 32]. In particular, the case of diffusive dynamics, similar to ours, is considered in [50] and [51], which aim at inferring the network topology based on the (cross-) power spectral density of a consensus process driven by noise, and a node knockout strategy. Methods combining spectral identification with optimization techniques are studied in [48, 59, 46, 49].

In all of the above works, however, the focus is on inferring the exact network topology, i.e., the complete specification of the network in terms of its adjacency or Laplacian matrices. In contrast, our goal is not to infer a specific network coupling, but rather to identify a statistical model thereof. Stated differently, we wish to infer a distribution over networks that is compatible with the available observations. This scenario is particularly relevant if our ability to reconstruct the network is limited by the amount of samples that we can obtain. Closest to our approach is the conference paper [45], in which a related observational model is discussed, with a new adjacency matrix drawn for each sample in (1), and which provides results about the asymptotic consistency of partition recovery. In contrast, we consider the case of a constant but unknown adjacency matrix, provide non-asymptotic results in the number of samples and, in addition, we address the problem of parameter estimation.

Similar to this paper, the works [61, 60] study the ‘blind solution’ of a ratio-cut problem (a common heuristic for partitioning graphs into dense, assortative clusters) based on observations of the output of a low-pass graph filter, excited by a low-rank excitation signal. They provide results on how well the ratio cut-problem can be approximated on arbitrary, fixed graphs, under the assumption of a low-rank excitation structure. In contrast, here we aim to infer an explicit generative graph model and can thus provide probabilistic performance guarantees about the recovery of the true planted partition. The work by Hoffmann et al. [18] and Peixoto [39] discuss hierarchical Bayesian models and inference schemes for related problems, but do not provide a theoretical analysis or statistical guarantees. Other related problems include the study of multi-reference alignment and synchronization problems (see, e.g., [40, 4]), or the study of state compression in Markov chains [67]. See also [33] for an example in which network reconstruction is attempted based on dynamics, but without knowledge of temporal information.

Finally, the mathematical machinery that we use has connections to ideas from spectral methods for clustering on graphs and SBM recovery [58, 42, 41, 1], and is related to the statistics literature on so-called “spiked” random matrices, with a low-rank signal component. In particular, our development relies on recent results regarding the efficiency of estimating low-rank covariance matrices [7].

1.4 Outline

The rest of the paper is structured as follows. After briefly reviewing relevant preliminaries and notation in Section 2, in Section 3 we formally introduce the identification problem that we will study. We then present our main results and algorithmic ideas on recovering the latent partitions of the SBM in Section 4. Building on those insights, in Section 5 we discuss procedures for learning the parameters of the SBM. Section 6 contains the proofs of our main technical results, and in Section 7 we present computational experiments to numerically assess the performance of our methods. We conclude with a brief discussion, and outline possible avenues for future work.

2 Preliminaries

2.1 Networks

An undirected network $\mathcal{G}$ consists of a node set $\mathcal{V}$ of known cardinality $n$, and an edge set $\mathcal{E}$ of unordered pairs of elements of $\mathcal{V}$. The edges can be conveniently represented as entries of a symmetric adjacency matrix $\bm{A}\in\{0,1\}^{n\times n}$, such that $A_{ij}=A_{ji}=1$ if and only if $(i,j)\in\mathcal{E}$. We define the diagonal degree matrix $\bm{D}=\operatorname{diag}{(\bm{A}\bm{1})}$, where $\bm{1}$ denotes the vector of all ones and $\text{diag}(\bm{x})$ denotes the diagonal matrix whose entries are given by the components of the vector $\bm{x}$. We define the normalized adjacency matrix as $\bm{L}=\bm{D}^{-1/2}\bm{A}\bm{D}^{-1/2}$. Note that $\bm{L}$ is a shifted version of the normalized Laplacian matrix of the graph, and is sometimes referred to in the literature simply as the “normalized Laplacian” [42, 41].

2.2 The stochastic block model

The SBM is a latent variable model that defines a probability measure over the set of undirected networks of fixed size $n$, represented by an adjacency matrix $\bm{A}\in\{0,1\}^{n\times n}$. In an SBM, the network is assumed to be partitioned into $k$ nonempty groups of nodes, so that each node $i$ is endowed with a latent group label $g_{i}\in\{1,\ldots,k\}$. Given these latent group labels, each link $A_{ij}$ is an independent Bernoulli random variable that takes value $1$ with probability $\Omega_{g_{i},g_{j}}$ and value [math] otherwise; that is, the probability of a link between two nodes depends only on their group labels:111For simplicity, our model allows self-loops.

[TABLE]

To compactly describe the model, we collect link probabilities between the different groups in a symmetric affinity matrix $\bm{\Omega}=[\Omega_{ij}]\in[0,1]^{k\times k}$. Furthermore, we define a partition indicator matrix $\bm{G}\in\{0,1\}^{n\times k}$, with entries $G_{ij}=1$ if node $i$ is in group $j$, and $G_{ij}=0$ otherwise. Based on these definitions, we can write the expected adjacency matrix under the SBM, given the partition and the group labels as222Note the slight abuse of notation here with respect to the conditional expectation. We treat the partition indicator matrix $\bm{G}$ as a degenerate random variable with a deterministic, but unknown value.

[TABLE]

As $\bm{\mathcal{A}}$ encodes all the parameters of an SBM, we will in the sequel refer to such a stochastic block model as $\mathcal{M}(\bm{\mathcal{A}})$. Throughout the rest of the paper, and in order to avoid degenerate cases, we will always assume that $\bm{\mathcal{A}}$ has rank $k$. In particular, this excludes the case where two rows in $\bm{\mathcal{A}}$ are co-linear, i.e., the case where the connection probabilities of two groups are identical or a linear scaling of each other.333Our main results rely on the analysis of a $k$-dimensional dominant subspace in the observed data, which is induced by the structure of $\bm{\mathcal{A}}$. If $\bm{\mathcal{A}}$ had rank $p<k$, then our arguments would need to be adjusted accordingly.

The SBM leads to a network whose expected adjacency matrix has a low-rank block structure (rank $k$). Due to the special structure induced by Equation 2, there exists a simple relationship between the partition indicator matrix $\bm{G}$ and the eigenvectors of the normalized version of the expected adjacency matrix $\bm{\mathcal{L}}$, defined by

[TABLE]

The relationship between $\bm{G}$ and the eigenvectors of $\bm{\mathcal{L}}$ is given in Lemma 2.1 below, and proved in the Appendix. We denote the eigenvalues of $\bm{\mathcal{L}}$ by $\mu_{i}$, and assume that they are sorted in order of decreasing magnitude, so that $|\mu_{1}|\geq|\mu_{2}|\geq\cdots\geq|\mu_{k}|>|\mu_{k+1}|=\cdots=|\mu_{n}|=0$.

Lemma 2.1**.**

Consider the matrix $\bm{V}=[\bm{v}_{1},\ldots,\bm{v}_{k}]$ of eigenvectors of $\bm{\mathcal{L}}$ corresponding to the $k$ nonzero eigenvalues. Then, there exists a $k\times k$ orthogonal matrix $\bm{\mathcal{U}}$ such that

[TABLE]

Moreover, if we define the row-normalized eigenvector matrix $\widetilde{\bm{V}}=\operatorname{diag}(\bm{V}\bm{V}^{\top})^{-1/2}\bm{V}$, then

[TABLE]

Given an SBM with an expected adjacency matrix $\bm{\mathcal{A}}$ as described above, let us denote the expected degree of node $i$ by $\delta_{i}:=\mathbb{E}[\bm{D}_{ii}]=[\bm{\mathcal{A}}\bm{1}]_{i}$, and let $\delta_{\text{min}}$ be the smallest expected degree. We will restrict attention to SBMs with a sufficiently large $\delta_{\text{min}}$.

Definition 2.2** (SBMs with concentration property).**

Given some $\epsilon>0$, we define $\mathcal{M}_{n}(\epsilon)$ as the class of SBMs $\mathcal{M}\left(\bm{\mathcal{A}}\right)$ with $n$ nodes, for which

[TABLE]

The condition on the minimal degree essentially restricts $\mathcal{M}_{n}(\epsilon)$ to correspond to models such that all nodes that belong to the same group belong to the same connected component, with high probability.444Note that for the special case of a single group, we are dealing with an Erdős-Rényi random graph, and it is well known that if $\delta_{\text{min}}>c{\ln(n)}$, with $c>1$, the random graph is connected with high probability, as $n\rightarrow\infty$. This is what will allow us to obtain reliable group membership information from the link patterns. Indeed, without such an assumption, and even if we knew the adjacency matrix $\bm{A}$, the exact recovery (with high probability) of the partition would be impossible [1].

For the above class of models the following concentration result can be established.

Lemma 2.3** (Concentration for normalized adjacency matrix).**

Fix some $\epsilon>0$. Consider a normalized adjacency matrix $\bm{L}:=\bm{D}^{-1/2}\bm{A}\bm{D}^{-1/2}$, computed from an adjacency matrix $\bm{A}$ generated according to an SBM model $\mathcal{M}(\bm{\mathcal{A}})$. If $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$, then the following bound holds with probability at least $1-\epsilon$:

[TABLE]

*where we recall that ${\bm{\mathcal{L}}:=\mathbb{E}[\bm{D}]^{-1/2}\mathbb{E}[\bm{A}]\mathbb{E}[\bm{D}]^{-1/2}}$. *

Lemma 2.3 follows from a result by Chung and Radcliffe [9]. The proof, which parallels [9], is given in the Appendix for completeness. In the following, we will be primarily interested in the case where $\delta_{\text{min}}=\omega\left(\ln(n/\epsilon)\right)$. In this case, when $n$ is large, the condition in Definition 2.2 will be satisfied.555More precisely, the condition $\delta_{\text{min}}=\omega\left(\ln(n/\epsilon)\right)$ means that $\delta_{\text{min}}/{\ln(n/\epsilon)}\rightarrow\infty$, as $n\rightarrow\infty$. Lemma 2.3 then guarantees a tight spectral concentration for large networks, i.e., that $\left\|\bm{L}-\bm{\mathcal{L}}\right\|$ will be small, with high probability.

2.3 Sub-Gaussian random vectors

For convenience, we recall here the definition of sub-Gaussian random variables and vectors.

Definition 2.4** (Sub-Gaussian random variables and vectors).**

*A zero-mean random variable $x\in\mathbb{R}$ is sub-Gaussian if there exists a constant $\sigma>0$ such that $\mathbb{E}[\exp(tx)]\leq\exp(t^{2}\sigma^{2}/2)$, for all $t\in\mathbb{R}$. A zero-mean random vector $\bm{x}\in\mathbb{R}^{n}$ is sub-Gaussian if for any non-random $\bm{u}\in\mathbb{R}^{n}$, the random variable $z=\bm{u}^{\top}\bm{x}$ is sub-Gaussian. *

Note that any Gaussian random variable or vector is sub-Gaussian. We also define the sub-Gaussian norm as follows.

Definition 2.5** (Sub-Gaussian norm [7]).**

*The sub-Gaussian norm of a sub-Gaussian random variable $x$ is $\|x\|_{\Psi_{2}}:=\sup_{k\geq 1}k^{-1/2}(\mathbb{E}|x|^{k})^{1/k}$. For a sub-Gaussian random vector $\bm{x}$, the sub-Gaussian norm $\|\bm{x}\|_{\Psi_{2}}$ is defined as ${\|\bm{x}\|_{\Psi_{2}}:=\sup_{\bm{u}\in\mathbb{R}^{n}\backslash\{0\}}\|\bm{u}^{\top}\bm{x}\|_{\Psi_{2}}/\|\bm{u}\|}$. *

3 The blind identification problem

In this section, we introduce the specific inference problem that we will study. We first provide a brief discussion on the class of dynamical models we consider, and then describe the precise problem setup.

3.1 Dynamical model

We consider the dynamical model described in Eq. 1, where the mapping $f$ is of the form

[TABLE]

for some known constants $\alpha_{\ell}$, and a known degree $R$, where $\bm{L}=\bm{D}^{-1/2}\bm{A}\bm{D}^{-1/2}$ is the (unknown) normalized adjacency matrix of the (unobserved) underlying graph. We assume that the system matrix satisfies $\|f(\bm{A})\|\leq 1$, and corresponds to a (marginally) stable dynamical system.

The above mapping $f$ can be considered to describe some local interactions between the nodes, in that it allows for interactions between nodes that are up to $R$ steps away in the graph. The mapping Eq. 6 $f$ may alternatively be considered from the perspective of graph signal processing [35, 30] as a graph filter of order $R$, that acts (instantaneously) on the state variables of the dynamical system Eq. 1 at each time step.

To make our discussion more concrete we will in the following concentrate on the dynamical system

[TABLE]

as our main example. This discrete-time system, and by extension the dynamics in (8), is closely related to diffusion processes on networks [29, 47]. Indeed, $\bm{L}$ is related via a similarity transformation to the (column-stochastic) transition matrix $\bm{A}\bm{D}^{-1}$ (or diffusion kernel) of an unbiased random walk on the network. Such random walk dynamics have been used to model a variety of processes, such as the spreading of rumors and opinions, and have found applications in other areas of network analysis such as ranking nodes or detecting communities [29]. Another interpretation of the model is in terms of (discrete-time) consensus processes, which are distributed averaging processes that have been considered extensively in the Control Theory literature [57, 21, 34].

It is well known that $\|\bm{A}\bm{D}^{-1}\|=1$ and accordingly $\|\bm{L}\|=1$, since $\bm{L}$ and $\bm{A}\bm{D}^{-1}$ are related via a similarity transform. Hence, the assumption $\|f(\bm{A})\|\leq 1$ is fulfilled by our example dynamics. We will make use of this fact in our calculations later.

Remark 3.1**.**

*By allowing for a large enough $R$ we can view the observation $\bm{x}_{T}$ at time $T$ as a signal $\bm{x}_{0}$ filtered by a graph filter with mapping $f^{\prime}(\bm{A})$, where $f^{\prime}(\bm{A})$ is obtained by composing $T$ copies of the original function $f$. However, we refrain from considering this ‘static’ filter perspective here, as this would disguise a trade-off between the sampling time and the number of samples required for our identification problem, which will become apparent in the following sections. *

3.2 Problem description

Let us provide a detailed breakdown of the observation model we consider.

We have a discrete-time dynamical system with state vector $\bm{x}_{t}\in\mathbb{R}^{n}$, which evolves according to

[TABLE]

where $f$ is of the form Eq. 6. 2. 2.

We assume that the normalized adjacency matrix $\bm{L}$ that appears in $f(\bm{A})$ is generated by first drawing a graph from an SBM with $k$ blocks, (unknown) partition indicator matrix $\bm{G}$, and some (unknown) affinity matrix $\bm{\Omega}$, and then normalizing its adjacency matrix $\bm{A}$ to obtain $\bm{L}$ and $f(\bm{A})$. 3. 3.

We obtain $s$ samples $\bm{x}_{T}^{(i)}$, $i=1,\ldots,s$, at a known time $T$, associated to $s$ different trajectories. 4. 4.

The (unobserved) initial vectors $\bm{x}_{0}^{(i)}$ for the different samples are i.i.d., drawn from a zero-mean sub-Gaussian vector distribution, where the components of the vector are uncorrelated, i.e., $\mathbb{E}[\bm{x}_{0}^{(i)}[\bm{x}_{0}^{(i)}]^{\top}]=\bm{I}$.

Under this setup, we consider the following two problems, which we refer to as the partition recovery and parameter estimation problems.

Problem 3.2**.**

*(SBM partition recovery from dynamics). Given $s$ independent samples, each one drawn at a fixed time $T$ according to the above described model, recover the latent partition of the nodes into blocks. *

Problem 3.3**.**

*(SBM parameter estimation from dynamics). Given $s$ independent samples, each one drawn at a fixed time $T$ according to the above described model, recover the affinity matrix $\bm{\Omega}$ describing the link probabilities in the underlying SBM. *

Note that, as the assignment of labels to groups is arbitrary, we can only hope to recover the block structure and the affinity matrix $\bm{\Omega}$ up to a permutation of the group labels. Similarly, as will become clear in the following sections, the parameters of the SBM will in general not be identifiable, unless some additional assumptions are made on the underlying model or some additional side information is available on the structure of the network.

Moreover, it should be apparent that the difficulty of the above problems will depend on the specific realization of the initial conditions $\bm{x}_{0}^{(i)}$ and of the adjacency matrix $\bm{A}$. Thus, we cannot expect accurate inference for every possible realization. Instead we look for algorithms and results that, under certain assumptions, guarantee accurate inference with high probability.

We make the following technical assumption about the data generation process. {assumption}[Bounded moments assumption] We assume that there exists an absolute constant $c_{0}>0$ such that $\mathbb{E}[(\bm{u}^{\top}\bm{x}_{0}^{(i)})^{2}]\geq c_{0}\|\bm{u}^{\top}\bm{x}_{0}^{(i)}\|^{2}_{\Psi_{2}}$, for all $\bm{u}$.

As discussed in [7], Assumption 3.3 effectively bounds all moments of $\bm{x}_{0}^{(i)}$ as a function of the second moment of $\bm{x}_{0}^{(i)}$. As a special case, this assumption is satisfied by a standard Gaussian vector, with $c_{0}=\pi/2$ [7]. 3.3 will enable us to invoke (via Lemma 6.2) the results of [7] on the estimation of low-rank covariance matrices, thereby facilitating our analysis.

Remark 3.4**.**

*If the number of available samples becomes large ($s\rightarrow\infty$), we could in principle try to infer the exact adjacency matrix $\bm{A}$. However, we are interested in a situation where $s$ is comparatively small ($s\ll n$), and our goal is to infer the SBM model that generated $\bm{A}$ (and thus $\bm{L}$), that is, the latent partition and the affinity matrix $\bm{\Omega}$. *

4 Partition recovery

In this section we first present our main results and an algorithm for solving the partition recovery problem (Problem 1) and then provide a broader discussion of the underlying algorithmic ideas.

4.1 Main results

Our algorithm for partition recovery relies on spectral clustering, and is displayed as Algorithm 1.

It builds upon the idea that the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$ at time $T$, defined by

[TABLE]

where $\bm{\bar{x}}_{T}=\frac{1}{s}\sum_{i=1}^{s}\bm{x}_{T}^{(i)}$ is the sample mean, can serve as a proxy for what we call the ensemble covariance matrix

[TABLE]

Because the initial vector has identity covariance, it is seen that $\bm{S}(T)$ would be the covariance matrix of the vector $\bm{x}_{T}^{(i)}$ if the dynamics were evolving according to the matrix $f(\bm{\mathcal{A}})$, i.e., if we were to replace the adjacency matrix of the sampled network with its expected value. Intuitively, the matrix $f(\bm{\mathcal{A}})$ encodes the partition structure of the SBM, and since $\bm{S}(T)$ is simply a matrix power of $f(\bm{\mathcal{A}})$, it has the same dominant eigenspace. Being able to compute a close proxy for $\bm{S}(T)$ enables us to gain information about the partition structure. Thus, as long as the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$ is close to the ensemble covariance matrix $\bm{S}(T)$, we can expect good recovery performance.

To gain some first insight, Theorem 4.1 below addresses the relation between these two covariance matrices, for our example dynamics $f(\bm{A})=\bm{L}$.

Theorem 4.1**.**

Let $\bm{L}$ be the normalized adjacency matrix constructed from a graph drawn from an SBM $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$ with $n$ nodes and minimal degree $\delta_{\text{min}}$. Let $\bm{x}_{0}^{(i)}$ $i=1,\ldots,s$, be i.i.d. zero-mean sub-Gaussian vectors that satisfy 3.3, with covariance matrix $\bm{I}$.

Then, for the dynamics of the form Eq. 7, and with probability at least $1-5s^{-1}-2\epsilon$, we have:

[TABLE]

where $\widehat{\bm{\Sigma}}(T)$ is defined in (9), $\bm{S}(T)$ is defined in (10),

[TABLE]

$C$* is some absolute constant, and where $\mu_{i}$ are the eigenvalues of $\bm{\mathcal{L}}$ ordered by decreasing magnitude. *

The proof of Theorem 4.1 is given in Section 6.1.666The probability $1-5s^{-1}-2\epsilon$ in the statement of the theorem can be improved to $1-5s^{-1}-\epsilon$ by a more careful use of the union bound in the proof. The same comment applies to Theorems 4.2 and 4.3 below.

For some interpretation, note that the term $2TM$ in (11) is independent of the number of samples. As we will see below, this term arises because the graph on which the dynamics takes place is only a proxy for the expected graph under the SBM, and is drawn only once. We may think of this error term as a ‘baseline’ insurmountable error that we will incur no matter how many samples we collect. In light of the condition on $\delta_{\rm min}$ in Definition 2.2, we have $M<1$, but nevertheless the term $2TM$ in the error bound increases (linearly) with time. Hence, the later we make our observations (larger $T$), the worse the first term in our bound becomes. From the perspective of our model (8), this is intuitive: since $\lambda_{i}(\bm{L})\leq 1$ for all $i$, the relevant signal that we observe will be attenuated over time, making it more difficult to extract information about the embedded partition from our observations when $T$ is large.

In contrast to the first term, the second term $B(T)$ in (11) depends both on the time $T$ and the number of samples $s$. The sample dependence is related to the convergence rate of the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$ to the true covariance. The time-dependence of $B(T)$ is again related to $M$, but in a non-linear fashion. Since $M<1$, the error term $B(T)$ decreases with time.

The above discussion reveals a trade-off. If $T$ is large, then the first term, $2TM$, is large and dominates. If $T$ is small, then the second term, $B(T)$, can be large, especially if the number of samples, $s$, is moderate. In general, an intermediate value of $T$ will be yielding a most favorable upper bound.

It is also insightful to consider the behavior of our bound for both the asymptotic regime in which we have a large number of samples $s\rightarrow\infty$ (even though this is not the setting of interest to us), as well as in the limit of large networks ($n\rightarrow\infty$). If we consider the case $s\rightarrow\infty$, and in the regime where $\delta_{\text{min}}=\omega(\ln(4n/\epsilon))$, then the error term is optimal, exactly in line with the requirements of “exact recovery” of the partition for an SBM where the network (the matrix $\bm{A}$) is observed [1]. Similarly, assuming a sufficiently large $\delta_{\text{min}}$ as $n\rightarrow\infty$, the dependency on the number of samples is aligned with the optimal rate for the estimation of low-rank covariance matrices [7].

For the general case, where $f$ is given by Eq. 6, an analogous result can be established.

Theorem 4.2**.**

Let $\bm{L}$ be the normalized adjacency matrix constructed from a graph drawn from an SBM $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$ with $n$ nodes and minimal degree $\delta_{\text{min}}$. Let $\bm{x}_{0}^{(i)}$ $i=1,\ldots,s$, be i.i.d. zero-mean sub-Gaussian vectors that satisfy 3.3, with covariance matrix $\bm{I}$.

Then, for the dynamics of the form Eq. 8, and with probability at least $1-5s^{-1}-2\epsilon$, we have:

[TABLE]

where $\widehat{\bm{\Sigma}}(T)$ is defined in (9),

[TABLE]

$C$* is some absolute constant, and where $\zeta_{i}$ are the eigenvalues of $f(\bm{\mathcal{A}})$ ordered by decreasing magnitude. *

The proof of this result is given in Section 6.2.

Based on the concentration result in Theorem 4.2, we can establish the following guarantee for the partition recovery problem. This result is phrased in terms of the misclassification rate $q:=|\mathbb{M}|/n$, where $\mathbb{M}$ is the set of “misclassified nodes” formally defined in  Section 6.3.777Importantly, this definition of misclassified nodes assumes that we can solve the $k$-means procedure exactly. In practice, however, only a $(1+\epsilon_{0})$-approximation may be obtained efficiently [24].

Theorem 4.3**.**

Let $\xi_{1}\geq\cdots\geq\xi_{k}$ be the top $k$ eigenvalues of the ensemble covariance matrix $\bm{S}(T)=f(\bm{\mathcal{A}})^{2T}$, which we assume to be positive. Let $\bm{W}$ and $\bm{V}$ denote the matrices formed by the top $k$ eigenvectors of $\widehat{\bm{\Sigma}}(T)$ and $\bm{S}(T)$, respectively, and denote their row-normalized versions by $\widetilde{\bm{W}}$ and $\widetilde{\bm{V}}$. Define $\tau$ to be the minimum 2-norm of any row of $\bm{W}$ and $\bm{V}$. Then, with probability $1-5s^{-1}-2\epsilon$ the following bound for the misclassification rate $q$ holds:

[TABLE]

*where $B_{f}(T)$ is as defined in Theorem 4.2. *

The proof is given in Section 6.3.2. Theorems 4.1, 4.2, and 4.3 provide performance guarantees that show the interplay between the number of groups $k$, the network size $n$, and the spectral properties of the network, as encapsulated in $\mu_{k}$ and $\|\widehat{\bm{\Sigma}}(T)-\bm{S}(T)\|$. Before continuing, let us add a few remarks on some aspects of these results.

The time $T$ at which the samples are obtained is not a required input to the algorithm (cf. Algorithm 1). However, as is apparent from Theorems 4.1, 4.2, and 4.3, the partition recovery performance is influenced by the sampling time. In particular, since the eigenvalues $\xi_{k}$ of the ensemble covariance matrix $f(\bm{\mathcal{A}})^{2T}$ can only decrease with time, the bound we can obtain for the misclassification rate is weaker in comparison to the concentration result for the covariance matrix (Theorems 4.1 and 4.2). 2. 2.

In Algorithm 1, we require the number of blocks $k$ to be an input to our algorithm. In principle, we can try to infer this parameter from the data, e.g., by looking at some appropriate statistics of the eigenvalue gaps; see for example [63, 58, 44], or the discussion in [7] on the detection of eigenvalue jumps in covariance matrix estimation. However, we will not carry out this line of analysis in this paper. 3. 3.

The same approach can be used for other types of (linear) operators derived from a low-rank statistical network model such as the SBM. For instance one may consider inference for other latent graph models, such as degree-corrected SBMs [10, 23].

Remark 4.4**.**

*Note that in Theorem 4.3, we assume that the first $k$ eigenvalues of $\bm{S}(T)$ are positive. Under the premise that we are dealing with a $k$-group SBM $\bm{\mathcal{A}}$ with rank $k$, this will be generically true. However, if $\text{rank}(\bm{\mathcal{A}})=p<k$, then the above Theorem would need to be adjusted, as only the first $p$ eigenvectors would carry information about the latent block structure, leading to possible non-identifiability. The case $\text{rank}(\bm{\mathcal{A}})<k$ could arise, for example, if we consider a block model with equal-sized groups, in which two rows of the affinity matrix $\bm{\Omega}$ are co-linear. For a related discussion on some of these issues see, e.g., [15]. Less trivial examples are also possible, e.g., if the rank of $\bm{\mathcal{A}}$ is $k$ but $\text{rank}(f(\bm{\mathcal{A}}))<k$. *

4.2 Algorithm and proof ideas

We now provide a more technical discussion of the proposed algorithm for partition recovery from a high-level perspective. The proofs in Section 6 will make these ideas mathematically precise. To simplify our discussion below, we will describe the algorithmic ideas in terms of the simple dynamics $f(\bm{A})=\bm{L}$. The general case is completely analogous.

The intuition underpinning our algorithm is that while we cannot observe the edges of the network, the network structure affects the second moments of the observed samples in ways that we can exploit.

Our analysis relies on a comparison of three different covariance matrices: (a) the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$, defined in (9); (b) the ensemble covariance matrix $\bm{S}(T)$, defined in (10); and (c) the covariance matrix $\bm{\Sigma}(T)$, defined by

[TABLE]

where the first equality uses the zero-mean property of $\bm{x}_{T}$, and the last equality uses the symmetry of $\bm{L}$.

As a first step, we establish that the sample covariance $\widehat{\bm{\Sigma}}(T)$ serves as a good estimator of $\bm{\Sigma}(T)$. This is done by leveraging the fact that the effective rank $r_{e}(\bm{\Sigma}(T))=\text{trace}(\bm{\Sigma}(T))/\|\bm{\Sigma}(T)\|$ of $\bm{\Sigma}(T)$ is low, and then showing that $\widehat{\bm{\Sigma}}(T)$ indeed concentrates around $\bm{\Sigma}(T)$, even when the number of samples is relatively small (Lemma 6.1).

Second, we observe that the eigenvectors of $\bm{\Sigma}(T)$ are the same as the eigenvectors of $\bm{L}$, although they now correspond to the modified eigenvalues $\lambda_{i}(\bm{\Sigma}(T))=\lambda_{i}^{2T}(\bm{L})$. This transformation does not change the ordering of the eigenvalues in terms of their absolute values and thus the dominant eigenspaces of $\bm{\Sigma}(T)$ and $\bm{L}$ are the same. In particular this implies that the dominant eigenspace of $\bm{\Sigma}(T)$ allows us to approximate the dominant eigenspace of $\bm{L}$.

Finally, we show that $\bm{L}$ concentrates around $\bm{\mathcal{L}}$ (cf. Lemma 2.3), so that $\bm{\Sigma}(T)$ also concentrates around $\bm{S}(T)$. In particular, the dominant eigenspace of $\bm{\Sigma}(T)$ allows us to approximate the dominant eigenspace of $\bm{S}(T)$ and of $\bm{\mathcal{L}}$, which in turn allows us to estimate the partition structure.

From the above discussion we see that Algorithm 1 is subject to two sources of randomness: (a) the normalized adjacency matrix $\bm{L}$ is an approximation of its expected version $\bm{\mathcal{L}}$; and (b) the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$ is an approximation of the true covariance matrix $\bm{\Sigma}(T)$. Theorems 4.1 and 4.2 allow us to bound these two effects, leading to the performance guarantees in Theorem 4.3.

Remark 4.5** (Computational considerations).**

*In Algorithm 1, we first form the sample covariance and then analyze its spectral properties via an eigenvalue decomposition. Equivalently, we can employ an SVD of the normalized and centered sample matrix $\bm{X}=\frac{1}{\sqrt{s}}[\bm{\widetilde{x}}_{T}^{(1)},\ldots,\bm{\widetilde{x}}_{T}^{(s)}]$, where $\bm{\widetilde{x}}_{T}^{(i)}=\bm{x}_{T}^{(i)}-\bm{\bar{x}}_{T}$. Using the SVD $\bm{X}=\bm{U}\bm{\Theta}\bm{Q}^{\top}$ it is easy to see that $\widehat{\bm{\Sigma}}(T)=\bm{X}\bm{X}^{\top}=\bm{U}\bm{\Theta}^{2}\bm{U}^{\top}$ and hence all the information required by the algorithm is contained in the SVD of $\bm{X}$. While mathematically equivalent, working directly with the SVD of $\bm{X}$ avoids forming $\widehat{\bm{\Sigma}}(T)$ explicitly, and thus reduces the memory requirements when $s\ll n$. *

5 Parameter estimation

In this section we turn to the parameter estimation problem. We first provide an overview of the general algorithmic idea with the help of an example. We then discuss certain obstacles to its solution, including some insurmountable non-identifiability issues. In Section 5.2 we then provide a generic “meta”-algorithm for the solution of the parameter identification problem. Finally, in Section 5.3, we provide some theoretical guarantees for its performance, for a planted partition model, which is a special case of the SBM.

5.1 Algorithm ideas

To estimate the parameters of the underlying SBM we can use a similar strategy as for the recovery of the partitions. The ensemble covariance matrix $\bm{S}(T)=f(\bm{\mathcal{A}})^{2T}$ is a transformation of the normalized adjacency matrix. Hence, it contains relevant information not only about the partitions, but also about the model parameters. Let us illustrate this fact with our example dynamics $f(\bm{A})=\bm{L}$ on a graph drawn from of a planted partition model, a simple form of an SBM.

Example 5.1**.**

Consider a special case of the SBM, known as the planted partition model. Conditional on the group labels, the planted partition model with $n$ nodes is described by probabilities $a/n$ and $b/n$ (with $a\neq b$) of linking to a node in the same group or a different group, respectively.888In Theorem 5.3 below, we will assume that $a,b$ and $n$ are such that $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$. In the context of this example, this means $a$ and $b$ are assumed to be of order at most $\ln(n)$. The affinity matrix for such a model, with $k$ blocks, is of the form $\bm{\Omega}=\frac{(a-b)}{n}\bm{I}_{k}+\frac{b}{n}\bm{1}_{k}\bm{1}_{k}^{\top}$, where $\bm{I}_{k}$ and $\bm{1}_{k}$ denote the $k$-dimensional identity matrix and the vector of all ones, respectively. Let us consider such a planted partition matrix with equally-sized groups. This leads to an expected matrix $\bm{\mathcal{L}}$ with entries

[TABLE]

where, we recall, $g_{i}$ represents the latent group label of node $i$.

Based on these calculations we can find expressions for ${\bm{S}}(T)=\bm{\mathcal{L}}^{2T}$:

[TABLE]

This example suggests the following generic strategy for the parameter recovery problem. First, we can derive the functional dependencies of the entries of $\bm{S}(T)$ on the latent parameters. Let us denote these relations by $F(\bm{\Omega},\bm{G})=\bm{S}(T)$. In the above example, $F(\bm{\Omega},\bm{G})$ is given by Equation 15. Note that it is possible to derive such relationships $F(\bm{\Omega},\bm{G})=\bm{S}(T)$ for any SBM with generic parameters $\bm{\Omega}$ and $\bm{G}$, even though the functional relationships will be more complicated. Second, in principle, we can then use the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$ as an approximation of $\bm{S}(T)$ and estimate the parameters by leveraging the functional relationships $F(\bm{\Omega},\bm{G})\approx\widehat{\bm{\Sigma}}(T)$. However, there are some difficulties to overcome to render this general idea operational.

First, the derived equations might be insufficient to identify the model parameters. Indeed, notice that the equations in (15) do not enable us to identify both model parameters $a$ and $b$ of the planted partition model. This is an artifact of the symmetry properties of the normalized adjacency matrix $\bm{L}$. In particular, this degeneracy arises from an invariance in the dynamical model (8). Indeed, if in Example 5.1 we replace $(a,b)$ by $(\gamma a,\gamma b)$, where $\gamma$ is a positive constant, the expected matrix $\bm{\mathcal{L}}$ remains the same. Thus the model parameters are not recoverable, even if we had access to the exact matrix $\bm{\mathcal{L}}$. In the context of our specific model (8), we thus have to make some additional assumptions to obtain a viable parameter estimation approach.999For a planted partition model with two equally-sized blocks there is additional symmetry in the problem in that we can exchange $(a,b)$ with $(b,a)$ and obtain exactly the same set of equations. In this case an additional assumption, e.g., of the form $a>b$ is thus necessary for identifiability. We remark that for other dynamical models such additional assumptions may not be necessary.

A second issue is that the obtained equations will be typically inconsistent if we just plug in the sample covariance matrix. Thus, instead of simply inverting the functional relations $F(\bm{\Omega},\bm{G})\approx\widehat{\bm{\Sigma}}(T)$, we will have to find approximate solutions. There are many ways of obtaining an approximate solution in such cases, e.g., by formulating an optimization problem $\min_{\bm{\Omega}}\|F(\bm{\Omega},\bm{G})-\widehat{\bm{\Sigma}}(T)\|$. Which approach would have better statistical properties in the face of the noise present in $\widehat{\bm{\Sigma}}(T)$ is an interesting question. We do not go into such details, but provide instead a generic “meta”-algorithm in the next section.

A third possible issue with the above approach is that it requires the knowledge of the group labels. In the context of Example 5.1, based on (15) we can estimate the parameters $a$ and $b$ directly from the entries of $\bm{S}(T)$ – or rather its estimated version $\widehat{\bm{\Sigma}}(T)$ – but only if the group labels $g_{i}$ of the nodes are known. This implies that we would have to first recover the latent partition before we can estimate the parameters based on these relations. Although this is a reasonable procedure in practice, the recovery performance will depend strongly on how good the partition recovery is in the first place.

However, in certain cases we can avoid the latter problem as follows. In particular, we can express the eigenvalues of $\bm{S}(T)$ as a function of the model parameters, and obtain a set of equations $\lambda(\bm{S}(T))=F_{\lambda}(\bm{\Omega})$ for the parameters $\bm{\Omega}$. Notice that these equations are independent of the group labels of the nodes. The reason is that a permutation of the nodes defines a similarity transformation of $\bm{\mathcal{L}}$, and hence the eigenvalues of $\bm{\mathcal{L}}$ – and accordingly those of $\bm{S}(T)$ – are not affected by the ordering of the rows and columns of the matrix.

To make these ideas concrete let us reconsider our earlier example.

Example 5.2** (continued).**

In terms of spectral analysis, the nonzero eigenvalues of $\bm{S}(T)$ are:

[TABLE]

As expected, these equations are independent of the group labels. However, as outlined in our discussion above, as the equations are redundant, on their own they are insufficient to identify all the parameters, but some additional side-information is needed. Nevertheless, the fact that these equations are independent of the group labels makes this “eigenvalue-based” parameter inference strategy more amenable to theoretical analysis, as we will see in Section 5.3.

5.2 A meta-algorithm for parameter estimation

Despite the difficulties identified in the previous section, we can still obtain a viable parameter estimation approach. To avoid issues of non-identifiability, we will have to assume some additional, partial prior knowledge on the model parameters. More concretely, suppose that we know of some additional constraints on $\bm{\Omega}$ (or alternatively $\bm{\mathcal{A}}$), of the form $h(\bm{\Omega})=0$, for some known function $h$. Such constraints could be, e.g., on the expected density $\rho$ of the network, in which case $h(\bm{\Omega})$ would be of the form

[TABLE]

As discussed earlier, based on our dynamical model we can obtain relations of the form ${\bm{S}(T)=F(\bm{\Omega},\bm{G})}$ for the entries of the ensemble covariance matrix and $\lambda(\bm{S}(T))=F_{\lambda}(\bm{\Omega})$ for the eigenvalues. Taking these equations and the constraints $h(\bm{\Omega})$ together, we obtain a system of (generally redundant) equations of the form

[TABLE]

Our approach is to “solve” this system of equations, while replacing $\bm{S}(T)$ by its estimate $\widehat{\bm{\Sigma}}(T)$, as outlined in Algorithm 2. Observe that the sampling time $T$ enters explicitly in the functional relations that we exploit, and therefore, unlike the partition recovery problem (cf. Algorithm 1), $T$ must be known.

As discussed in the previous section, the word “solve” in the description of Algorithm 2 is to be understood as some not fully specified approach that yields an approximate solution of the system (17).

5.3 Theoretical guarantees for a planted partition model

In this section, we provide some theoretical results for a parameter estimation strategy in which we will ignore Equation 17a and concentrate on recovering the parameters solely based on the eigenvalues and the constraints. To this effect, we will leverage our concentration result in Theorem 4.1 to provide a bound between the estimated eigenvalues $\lambda_{i}(\widehat{\bm{\Sigma}}(T))$ and their ensemble counterparts $\lambda_{i}(\bm{S}(T))$. These bounds then translate immediately into ($T$-dependent) error-bounds on the parameter recovery.

For simplicity, we illustrate our approach by considering our example dynamics $f(\bm{\mathcal{A}})=\bm{L}$ on a graph drawn from a two-group planted partition model (see Examples 5.1 and 5.2) subject to a density constraint on the network and the assumption that $a>b$. Similar results could be derived for other dynamics as well, using analogous arguments, but they would depend on the (possibly more complicated) parametrization of the model.

Recall from Example 5.1 that the second eigenvalue of $\bm{\mathcal{L}}^{2T}$ satisfies

[TABLE]

Given also prior knowledge of the density $\rho=\frac{a+(k-1)b}{nk}$, we can recover the parameter $a$, for the case $k=2$ as follows:

[TABLE]

Then, to estimate $a$ from samples, we can replace $\bm{\mathcal{L}}^{2T}$ in (18) by the sample covariance matrix $\widehat{\bm{\Sigma}}(T)$. For this estimation strategy we can derive the following error bound.

Theorem 5.3**.**

Consider an SBM $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$ with two equally-sized blocks, and a parameter matrix of the form $\bm{\Omega}=\frac{(a-b)}{n}\bm{I}_{2}+\frac{b}{n}\bm{1}_{2}\bm{1}_{2}^{\top}$, with $a>b$. Assume that $\rho=(a+b)/(2n)$ is given. Consider the estimator of $a$ given by

[TABLE]

Then, using the notation of Theorem 4.1, the estimation error $\eta:=|\hat{a}-a|/a$ satisfies the upper bound

[TABLE]

*with probability at least $1-5s^{-1}-2\epsilon$. *

The proof of Theorem 5.3 is given in Section 6.4.

Note that the numerator in (20) indicates that the sampling time $T$ plays a fundamental role in the estimation accuracy of $\hat{a}$. Moreover, this dependency is similar to the one encountered in Theorem 4.3 for partition recovery. In particular, there is a term that increases linearly with $T$, i.e., $2TM$, and another term that may decrease with $T$, i.e., $B(T)$ (cf. the discussion after Theorem 4.1). The denominator in (20) can also be explained intuitively. A larger $\lambda_{2}(\bm{\mathcal{L}})$ is associated with a marked difference between $a$ and $b$ [cf. (16)], which translates into better delineated blocks within our SBM. In this scenario, estimating $a$ should be easier, and this is captured in (20).

Remark 5.4**.**

*The above outlined eigenvalue-based parameter inference strategy is independent of the partition recovery problem and thus enables a clean theoretical analysis of the estimation error. Had we considered instead the general approach based on Equation 17a, we would be solving simultaneously the partition recovery and the parameter estimation problems, and the estimation error would be affected by partition recovery errors. While a theoretical analysis of such an interplay may be feasible in principle, we do not pursue it further here, but instead illustrate it with the help of numerical experiments in Section 7. *

6 Theoretical analysis

In the following, we provide the theoretical analysis of the spectral inference algorithms (Algorithms 1 and 2) discussed above. We first prove the concentration result in Theorem 4.1, which underlies the success of our spectral algorithms. The proof of this result follows from a series of lemmas that we establish in Sections 6.1.2 and 6.1.1. Afterwards, we show in Section 6.3 how Theorem 4.1 implies the partition recovery performance guarantees stated in Theorem 4.3. We finally prove Theorem 5.3 in Section 6.4.

6.1 Proof of Theorem 4.1: Concentration of the sample covariance matrix

In the next two subsections we prove Theorem 4.1 through a series of lemmas. The proof uses the triangle inequality

[TABLE]

and separate bounds for the two terms on the right-hand side.

6.1.1 Bounding the first term in (21)

We now establish a bound for the distance ${\|\widehat{\bm{\Sigma}}(T)-\bm{\Sigma}(T)\|}$ between the true and the sample covariances.

Lemma 6.1**.**

Under the conditions in Theorem 4.1, we have

[TABLE]

*with probability at least $1-5s^{-1}-\epsilon$, where $B(T)$ is as defined in Theorem 4.1. *

For the proof of Lemma 6.1 we make use of the following result established in [7].

Lemma 6.2**.**

[7, Theorem 2.1]** Let $\bm{x}\in\mathbb{R}^{n}$ be a zero-mean sub-Gaussian vector satisfying 3.3, with covariance matrix $\bm{\Sigma}$. Let $\bm{x}^{(1)},\ldots,\bm{x}^{(i)}$ be i.i.d samples drawn from $\bm{x}$. Consider the sample covariance matrix ${\widehat{\bm{\Sigma}}:=\frac{1}{s}\sum_{i=1}^{s}[\bm{x}^{(i)}-\bm{\bar{x}}][\bm{x}^{(i)}-\bm{\bar{x}}]^{\top}}$, where ${\bm{\bar{x}}:=\frac{1}{s}\sum_{i=1}^{s}\bm{x}^{(i)}}$ is the sample mean. Then, with probability at least $1-5s^{-1}$ we have

[TABLE]

*where $r_{e}(\bm{\Sigma})={\sum_{i}\lambda_{i}(\bm{\Sigma})}{/}{\|\Sigma\|}$ is the effective rank of $\bm{\Sigma}$, $c$ is a constant that depends on $c_{0}$ in Assumption 3.3, and $\|\cdot\|_{F}$ denotes the Frobenius norm. *

Notice that we always have $\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|\leq\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|_{F}$ and, for the scenario considered here, $\|\bm{\Sigma}(T)\|=1$, as by construction the dominant eigenvalue of $\bm{L}$ equals $1$. Thus, Lemma 6.2 provides a useful tool to prove Lemma 6.1, provided that we can find an upper bound on the effective rank $r_{e}(\bm{\Sigma}(T))$.

Lemma 6.3** (Bound on the effective rank).**

Under the conditions stated in Theorem 4.1, the effective rank satisfies the bound

[TABLE]

*with probability at least $1-\epsilon$. *

Proof 6.4**.**

Recall that we denote the eigenvalues of $\bm{L}$ by $\lambda_{i}$ and the eigenvalues of $\bm{\mathcal{L}}$ by $\mu_{i}$. Recall our convention that the eigenvalues have been sorted in descending order, according to their magnitude, i.e., $|\lambda_{1}|\geq|\lambda_{2}|\geq\cdots\geq|\lambda_{n}|$ and $|\mu_{1}|\geq|\mu_{2}|\geq\cdots\geq|\mu_{n}|$, respectively. From the definition of the effective rank and the fact that $\|\bm{\Sigma}(T)\|=1$ it follows that

[TABLE]

We can decompose this sum into three parts as follows

[TABLE]

In order to relate the (random) eigenvalues $\lambda_{i}$ to quantities that can be expressed in terms of the model parameters, we make use of the triangle inequality ${|\lambda_{i}|\leq|\mu_{i}|+|\lambda_{i}-\mu_{i}|}$. Weyl’s Theorem (see e.g., [20]) states that $|\lambda_{i}-\mu_{i}|\leq\|\bm{L}-\bm{\mathcal{L}}\|$ for all $i$. Furthermore, Lemma 2.3 implies that $\|\bm{L}-\bm{\mathcal{L}}\|\leq M$, with probability at least ${1-\epsilon}$. Finally, we also know that $|\lambda_{i}|\leq 1$ for all $i$, since $\|\bm{L}\|=1$ (more generally, we assumed $\|f(\bm{A})\|\leq 1$). We can thus deduce that $|\lambda_{i}|\leq\min\{1,|\mu_{i}|+M\}$ for all $i$, with probability at least $1-\epsilon$. Using in addition the fact that $\mu_{i}=0$ for $i>k$ and $M<1$ if $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$, we conclude that

[TABLE]

*with probability at least $1-\epsilon$. *

Proof 6.5** (Proof of Lemma 6.1).**

*By combining Lemmas 6.2 and 6.3, Lemma 6.1 follows. *

6.1.2 Bounding the second term in (21)

We derive a bound for the difference $\|\bm{\Sigma}(T)-\bm{S}(T)\|$, between the true and the ensemble covariances. We start with the following lemma.

Lemma 6.6**.**

For any integer $t>0$, we have

[TABLE]

Proof 6.7**.**

We will show that for any integer $t>0$,

[TABLE]

Equation 25* then follows directly from repeated application of (26). To show (26), notice that*

[TABLE]

*where the last equality follows from the fact that $\|\bm{\mathcal{L}}\|=\|\bm{L}\|=1$. *

Proof 6.8** (Concluding the proof of Theorem 4.1).**

Recall that $\bm{S}(T)=\bm{\mathcal{L}}^{2T}$ and $\bm{\Sigma}(T)=\bm{L}^{2T}$, and that $\left\|\bm{L}-\bm{\mathcal{L}}\right\|\leq M$, with probability at least $1-\epsilon$ (Lemma 2.3). Lemma 6.6 then implies that

[TABLE]

*with probability at least $1-\epsilon$. We bound the two terms on the right-hand side of (21) using (22) and (27), respectively, and hence establish the validity of Theorem 4.1. *

6.2 Proof of Theorem 4.2

Proof 6.9**.**

The proof is analogous to the proof of Theorem 4.1. Instead of using a bound for $\|\bm{L}-\bm{\mathcal{L}}\|$, we develop a bound for $\|f(\bm{A})-f(\bm{\mathcal{A}})\|$, and obtain that, with probability at least $1-5s^{-1}-\epsilon$,

[TABLE]

*where the last inequality follows from Lemma 6.6. The rest of the argument remains the same, but with $M_{f}$ replacing $M$. *

6.3 Implications for partition recovery

In the following, we explain how the above concentration results translate into guarantees for the partition recovery problem. Section 6.3.1 provides intermediate results on the spectral geometry of the problem and discusses the precise definition of misclassified nodes. The proof of Theorem 4.3 is then given in Section 6.3.2.

6.3.1 Partition recovery problem

In order to characterize recovery performance, we first clarify what we mean by a node being clustered correctly. To do this we follow the literature on SBMs [42, 41]. Recall that within our partition recovery algorithm, we perform a $k$-means clustering procedure on the rows of the matrix $\widetilde{\bm{W}}$ of (row-normalized) dominant eigenvectors associated with the estimated covariance $\widehat{\bm{\Sigma}}(T)$. This procedure assigns each row $i$ to a centroid vector $\bm{k}_{i}$ corresponding to a cluster.

Let us denote the ‘true’ centroid vector that a node $i$ is assigned to according to the ensemble covariance $\bm{S}(T)$ by $\bm{\kappa}_{i}$, where $i=1,\ldots,n$. From the definition of $\bm{S}(T)$ it follows that each $\bm{\kappa}_{i}^{\top}$ corresponds to a row of the matrix $\widetilde{\bm{V}}$ of the (row-normalized) top $k$ eigenvectors of $\bm{\mathcal{L}}$. Moreover, we know from Lemma 2.1 that there are exactly $k$ different population centroids and that the association of rows to centroids corresponds to the true group structure of the underlying SBM. However, if some of the eigenvalues of $\bm{S}(T)$ are repeated, the eigenvector matrix $\widetilde{\bm{V}}$ will not be uniquely defined, even though the subspace spanned by these vectors will be unique. As we are interested in the distance between the vectors $\widetilde{\bm{V}}$ and their approximation via $\widetilde{\bm{W}}$, our definition of misclassification has to include an orthogonal matrix $\bm{\mathcal{Q}}\in\mathbb{R}^{k\times k}$ that minimizes the orthogonal Procrustes distance ${\|\widetilde{\bm{W}}-\widetilde{\bm{V}}\bm{\mathcal{Q}}\|_{F}}$ (see also [42, 41]). In particular, this minimality property of $\bm{\mathcal{Q}}$ will be essential for Lemma 6.15.

Definition 6.10** (Set of misclassified nodes [41]).**

The set of misclassified nodes $\mathbb{M}$ is defined as

[TABLE]

In words, this means that we call a node $i$ misclassified if the centroid $\bm{k}_{i}$ to which it has been assigned is closer to some other population centroid $\bm{\kappa}_{j}$ than it is to its true population centroid $\bm{\kappa}_{i}$. Note that this definition of misclassified nodes adopted from the literature on SBMs [42, 41] implicitly assumes that i) we can solve the $k$-means clustering problem exactly and ii) the empirical centroids based on $\widetilde{\bm{W}}$ will be centered around the (appropriately rotated) population centroids. We refer to the literature for additional discussion on the consistency of spectral clustering under the SBM [25], and limit theorems for the eigenvectors of random graphs (see, e.g., [2] for a survey on random dot product graphs).

By accounting for the special geometry of the eigenvectors $\widetilde{\bm{V}}$, we can establish the following necessary condition for a node to be misclassified.

Lemma 6.11**.**

[41*]** If a node $i$ is misclassified, then ${\|\bm{k}_{i}^{\top}-\bm{\kappa}_{i}^{\top}\bm{\mathcal{Q}}\|\geq 1/\sqrt{2}}$. *

Proof 6.12**.**

First, we know from Lemma 2.1 that $\widetilde{\bm{V}}=\bm{G}\bm{\mathcal{U}}$ for some orthogonal matrix $\bm{\mathcal{U}}$. Hence, all population centroids $\bm{\kappa}_{i}$ are of unit length and orthogonal to each other. In particular, for all $j$ and $i$ belonging to different blocks/groups, we have

[TABLE]

where $\bm{g}_{i}$ is the true population cluster indicator vector corresponding to node $i$ (see Lemma 2.1).

Second, by the triangle inequality and the orthogonality of $\bm{\mathcal{Q}}$, we know that for all $j$ and $i$ in different blocks, we have

[TABLE]

Plugging in the result from (30) into (31), we obtain that $\|\bm{k}_{i}^{\top}-\bm{\kappa}_{j}^{\top}\bm{\mathcal{Q}}\|\geq\sqrt{2}-\|\bm{k}_{i}^{\top}-\bm{\kappa}_{i}^{\top}\bm{\mathcal{Q}}\|$. From Definition 6.10, for a misclassified node $i$ there exists some $j$ such that

[TABLE]

*where the last inequality is true by our argument above. The result follows after some simple algebraic manipulations. *

To arrive at our result for the misclassification rate in Theorem 4.3, we further follow closely the literature on SBMs [42, 41, 1], and establish the following two lemmas (Lemmas 6.13 and 6.15). The first lemma provides us with a relationship between the difference of the cluster centroids obtained from the $k$-means clustering and the true population cluster on the one side, and the difference between the eigenvectors of $\widehat{\bm{\Sigma}}(T)$ and $\bm{S}(T)$ on the other side.

Lemma 6.13**.**

Denote the matrix of cluster centroids obtained by performing k-means on the rows of $\widetilde{\bm{W}}$ as ${\bm{K}:=[\bm{k}_{1},\ldots,\bm{k}_{n}]^{\top}\in\mathbb{R}^{n\times k}}$. Then,

[TABLE]

Proof 6.14**.**

By the triangle inequality,

[TABLE]

*From the definition of $k$-means, however, we know that $\|\widetilde{\bm{W}}-\bm{K}\|_{F}\leq\|\widetilde{\bm{W}}-\widetilde{\bm{V}}\bm{\mathcal{Q}}\|_{F}$ as $\bm{K}$ is chosen in order to minimize the distances of the centroids to the individual rows. Hence, the claim follows. *

The second lemma builds on the Davis-Kahan theorem [11, 64], and bounds the difference between the row-normalized eigenvectors of $\widehat{\bm{\Sigma}}(T)$ and $\bm{S}(T)$.

Lemma 6.15** (Davis-Kahan).**

Following the definitions in Theorem 4.3, we have

[TABLE]

*where we recall (see Theorem 4.3) $\xi_{k}$ is the $k$-th top eigenvalue of the ensemble covariance matrix, and $\tau$ is be the minimum 2-norm of any row of $\bm{W}$ and $\bm{V}$. *

Proof 6.16**.**

The fact that

[TABLE]

holds for the unnormalized versions of $\bm{W}$ and $\bm{V}$ is a restatement of the Davis-Kahan Theorem as can be found in, e.g., [64, Theorem 2].

By definition, we have $\widetilde{\bm{V}}=\operatorname{diag}(\bm{V}\bm{V}^{\top})^{-1/2}\bm{V}$ and $\widetilde{\bm{W}}=\operatorname{diag}(\bm{W}\bm{W}^{\top})^{-1/2}\bm{W}$. Let us define $\tau_{i}:=\min\{[\bm{V}\bm{V}^{\top}]_{ii}^{1/2},[\bm{W}\bm{W}^{\top}]_{ii}^{1/2}\}$, which is the smaller of the 2-norms of the $i$-th rows of $\bm{V}$ and $\bm{W}$, respectively. Based on this definition of $\tau_{i}$ we see that $\tau=\min_{i}\tau_{i}$. Now, as row-normalization corresponds to diagonal scaling we obtain

[TABLE]

6.3.2 Proof of Theorem 4.3

Based on the results established thus far, we can provide a proof of the bound of the misclassification rate $q:=|\mathbb{M}|/n$ as stated in Theorem 4.3. Our proof parallels results for the misclassification rate if the graph is observed [41].

Proof 6.17**.**

We define the set $\mathbb{M}_{2}:=\{i:\|\bm{k}_{i}^{\top}-\bm{\kappa}_{i}^{\top}\bm{\mathcal{Q}}\|\geq 1/\sqrt{2}\}\supset\mathbb{M}$, which includes all misclassified nodes according to Lemma 6.11. From Lemmas 6.13 and 6.15, and Theorem 4.1, it follows that

[TABLE]

6.4 Implications for parameter estimation

Here we show how our concentration results allow us to prove Theorem 5.3.

Proof 6.18**.**

Note that we can write the relative estimation error as:

[TABLE]

The first inequality follows from Jensen’s inequality since the function $f(x)=x^{1/2T}$ is concave. The second inequality is a direct application of Weyl’s inequality. Hence, using Theorem 4.1 it follows that with probability at least $1-5s^{-1}-2\epsilon$,

[TABLE]

7 Numerical experiments

In this section we present some computational experiments with our algorithms for partition recovery and parameter estimation, to demonstrate their performance in practice.

7.1 Partition recovery

To gain intuition about the qualitative features of the partition recovery problem, we first consider the following example SBM consisting of $5$ groups with varying sizes $\bm{n}_{g}=[200,300,200,300,280]^{\top}$, and affinity matrix

[TABLE]

The partition recovery performance is shown in Figure 1A for a varying number of samples $s$, as a function of the observation time $T$. Here, the recovery performance is measured via the overlap score $Z$ defined as:

[TABLE]

where $z_{\text{chance}}$ is the probability of guessing the group assignments correctly and $z_{\text{actual}}$ is the fraction of nodes that are correctly classified. The guessing probability is computed as $z_{\text{chance}}=max(\bm{n}_{g})/n$, where $\max(\bm{n}_{g})$ is the size of the largest group. Note that $z_{\text{chance}}$ is simply $\frac{1}{\#\text{groups}}$ for groups of equal size. The actual overlap $z_{\text{actual}}$ is computed by solving a linear $k\times k$ optimal assignment problem between the inferred group labels and the planted partition.

Overall, we can infer the partitions remarkably well using our simple spectral inference strategy, even though we never observe the actual network. The results in Figure 1A display a couple of noteworthy features consistent with our technical results, which are worth examining. Specifically, we see that there is a pronounced effect of the observation time: both for small as well as very large $T$, the performance is markedly worse than for intermediate times. Thus, if we have a fixed budget for the number of samples, choosing the ‘right’ observation time can improve recovery performance. Obtaining more samples, in contrast, has a more moderate effect and is most beneficial for small times.

At an intuitive level, these results can be explained by viewing the network’s effects as a linear filter. On the one hand, for short observation times, the original white noise input is only filtered weakly and thus contains still a large noise component, which is detrimental to recovery performance. On the other hand, as the output approaches its stationary state, the information that the signal carries about the network structure is completely washed out. For intermediate times, however, there is a trade-off between attenuating the noise within the observed signal, and not losing too much information about the partition structure. Theorems 4.1 and 4.3 quantify precisely how this trade-off is governed by the spectral features of the network.

Using a planted partition model with an affinity matrix of the form $\bm{\Omega}=\frac{(a-b)}{n}\bm{I}_{k}+\frac{b}{n}\bm{1}_{k}\bm{1}_{k}^{\top}$ (see Example 5.1), we further investigate (see Figure 1B-D). numerically how the partition recovery performance is affected by the group sizes, the number of groups, and the signal-to-noise ratio (SNR) defined by [1],

[TABLE]

We note that that $\text{SNR}=1$ corresponds to the well-known detectability limit in the sparse regime of the SBM, for asymptotically large networks, with $n\rightarrow\infty$ (see [1] for a detailed review).

Figure 1B shows that for a fixed number of groups (here $k=2$), the group sizes have a noticeable effect on recovery performance. On the other hand, if we fix the number of nodes, changing the number of groups does not have a noticeable effect. This hints at the fact that the number of nodes $n$ (and thus the number of nodes per group) is one of the key drivers of the difficulty of partition recovery. Finally, in Figure 1D we plot the influence of the SNR on the detection performance. As expected, a smaller SNR results in a dramatically worse performance. Remarkably, for $\text{SNR}=4$, using only $s=50$ samples leads to an almost perfect partition recovery, for a suitable observation time.

7.2 Model parameter estimation

We now demonstrate numerically that we can infer the model parameters of the SBM using only the observed outputs. For simplicity, we focus on the case of the planted partition model with $k=2$ groups of equal size, and assume that the average expected network density $\rho$, or equivalently the average expected mean degree $d_{\text{av}}=\rho n$, is known a priori (see also the discussion in the context of Example 5.1 and Section 5.3).

Since the estimation error is not directly related to the partition recovery performance, let us first focus on the eigenvalue-based parameter estimation strategy, using the estimator $\hat{a}$ for the parameter $a$ given by Equation 19 (see Section 5.3).

Figure 2 shows the results of applying this parameter estimation scheme in numerical experiments. To quantify the parameter recovery performance we compute the relative recovery error $\eta=|a-\hat{a}|/a$. Similar to partition recovery, we see in Figure 2A that the observation time has a strong effect on the recovery performance and we also see a marked decrease in the error from more than $10\%$ initially, to around $3\%$ for longer observation times.

Figure 2B-C demonstrates the effect of the group size and the SNR on recovery performance. Similar to the partition recovery problem, we observe that smaller group sizes (a smaller network) and a larger SNR lead to better performance.

Let us now consider a different parameter estimation strategy based on first recovering the groups in the SBM and then trying to recover the parameters based on the entries of the sample covariance matrix (see discussion in previous sections). We concentrate again on the planted partition model with two equally-sized groups and a given density $\rho=(a+b)/(2n)$ as above.

Specifically, we estimate the parameter using the following approach. First, we estimate the group labels $g_{i}$ of the nodes using Algorithm 1 and sort the indices of the sample covariance $\widehat{\bm{\Sigma}}(T)$ accordingly. Next, we take the $t$-th root of $\widehat{\bm{\Sigma}}(T)$ (via an SVD), and then compute the average magnitude of all off-diagonal terms that relate different groups, to obtain the following statistic:

[TABLE]

From Equation 15) we know that

[TABLE]

which implies that $x\approx-\frac{1}{n}\left(\frac{a-b}{a+b}\right)^{2}+\frac{1}{n}$. We then use the following estimator for the parameter $a$:

[TABLE]

Figure Figure 3 displays the performance of this estimator in terms of the error $\eta_{\text{part}}=|\hat{a}_{\text{part}}-a|/a$, in comparison to the relative error $\eta_{\lambda}=|\hat{a}-a|/a$ of the eigenvalue-based estimation strategy.

7.3 Real-world networks

In this section, we examine the performance of our methods on real-world networks that are not known to be generated from an SBM. In this final set of numerical experiments, we test two real-world networks, on which we assume a diffusion dynamics of the form (7) takes place.

The first network we study is the well-known Karate-club network which was initially charted by Zachary [66], describing the social interactions among 34 individuals of a Karate-club that eventually split into two factions. We use this split as our ‘ground-truth’ partition to compare against. Although there is a priori no reason for why the split should be correlated with the graph structure [38], this finding is now well-supported in the literature.

The results of applying our procedure for a range of different sampling times and number of samples available are shown in Figure 4A. Interestingly, we find that already with three samples, we can achieve a recovery performance that is as good as applying a corresponding spectral clustering procedure to the full network — which in this case corresponds to making one classification error with respect to the reported fission of the club.

As a second example, we study a social interaction network consisting of the friendship relationships between students at Caltech, a dataset that was originally studied by Traud et al. [55, 56] in a comparative study of a larger set of such networks. While there is again no clear ‘ground truth’ partition of the network, some additional node covariates such as gender or the major of the students are examined as part of their study [55, 56]. Traud et al. [55, 56] observed that in the context of the social network of Caltech students the node covariate ‘house membership’ (dormitories) correlates strongly with a putative community structure in the network, even though the house-membership cannot explain the structure of the network completely.

Motivated by these finding, we study the partition recovery performance with respect to the partition of the nodes induced by the house-membership data in Figure 4B-C. Here we limit ourselves to studying the largest connected component of the network, where we remove all the nodes for which no house-membership label is available. This results in a network of $594$ nodes split between 8 houses (groups). As the structure of the network is not completely explained by the house-membership data, we recover the induced partition only with an overlap score of $0.65$ (see Figure 4B). We also show our recovery performance relative to the partition we would obtain from the spectral clustering procedure that uses $\bm{L}$ (i.e., knowledge of the complete network) instead of $\widehat{\bm{\Sigma}}(T)$ for clustering. As can be seen in Figure 4C, our recovered partition is almost perfectly correlated with the result of the clustering-based procedure, highlighting that our method indeed performs as expected and that the low score with respect to the house-membership partition is largely due to a lack of a correlation between the house membership variable and the network block structure.

8 Discussion

Graph-based tools have become prevalent for the analysis of a range of different systems across the sciences. Network inference, the problem of determining the topology of a network based on a set of nodal observations, has accordingly become a paramount tool to enable this type of analysis. In this work, we have advocated a fresh look onto this problem for applications in which exact inference is impossible, e.g., due to sampling constraints. Instead of inferring the exact network in terms of its adjacency matrix, we infer a generative model of the network structure, thereby learning a probability distribution over networks. Importantly, for a low-dimensional statistical model such as the SBM discussed here, this procedure can be implemented using far fewer samples than would be required to infer the exact network structure.

The proposed method is especially relevant for scenarios in which we aim to analyze the network further from a statistical perspective, e.g., if our ultimate goal is to perform community detection on the network. An added benefit of having fitted a generative model is that we can employ the model for further tasks afterwards, e.g., to generate surrogate network data.

Most of our algorithmic ideas can be extended to more general functions $f$ which inherit the low-rank structure of the (expected) SBM, provided an appropriate concentration result for the system matrix under consideration can be derived. For instance, this includes those functions describing a low-pass graph filter, as considered in [61, 45], but also functions of the type $f(\bm{A})=\sum_{\ell}\alpha_{\ell}\bm{A}^{\ell}$ for some suitably chosen parameters $\alpha_{\ell}$.

To adjust the proposed algorithms to such setups, we would have to establish and prove appropriate relationships between the spectral properties of the system matrix on one side and the partition structure and parameters of the model on the other. Specifically, step 3 of Algorithm 1 uses the particular structure of the eigenvectors of $\bm{L}$ (see Lemma 2.1). Of course, the set of equations relating the model parameters to the statistics of the observed outputs will in general be different and will depend on the specific model chosen. In fact, some of the non-identifiability problems we encountered here in terms of recovering the model parameters may not arise for other classes of system matrices $f(\bm{A})$.

Our work also opens up several additional avenues for future research. This includes the analysis of the inference of more realistic network models based on a setup as discussed here, ranging from degree-corrected SBMs [23] all the way to graphons [28, 3]. Likewise, it will be interesting to extend the above results to more general dynamical processes, including nonlinear dynamics. Another appealing line of future work would be to consider scenarios in which only a part of the network is observed [32], or additional (control) inputs are possible.

Appendix A Omitted proofs

A.1 Proof of Lemma 2.1

In order to prove Lemma 2.1, we provide a slightly more general discussion about the relationship between the eigenvectors of $\bm{\mathcal{L}}$ and the partition indicator matrix $\bm{G}$. Lemma 2.1 then follows directly from Corollary A.6 below.

Definition A.1** (Normalized affinity matrix).**

Consider an SBM with affinity matrix $\bm{\Omega}$ and partition indicator matrix $\bm{G}$. Let us denote the vector of group sizes by $\bm{n}_{g}=\bm{G}^{\top}\bm{1}$ and define $\bm{N}_{g}:=\operatorname{diag}(\bm{n}_{g})=\bm{G}^{\top}\bm{G}$. We then define the normalized affinity matrix as

[TABLE]

Note that $\bm{\mathcal{L}}$ can be written as

[TABLE]

Lemma A.2**.**

*Let $\bm{\mathcal{L}}$ be the expected normalized adjacency matrix of an SBM with a corresponding normalized affinity matrix $\bm{\Gamma}$ with eigenvalues $\lambda_{i}$ and corresponding eigenvectors $\bm{u}_{i}$. Then, any eigenvector $\bm{u}_{i}$ of $\bm{\Gamma}$ gives rise to an eigenvector $\bm{v}_{i}=\bm{G}\bm{N}_{g}^{-1/2}\bm{u}_{i}$ of $\bm{\mathcal{L}}$ with the same eigenvalue. *

Proof A.3**.**

A simple computation shows that

[TABLE]

Corollary A.4**.**

*The partition indicator matrix $\bm{G}$ can be written as a linear combination of the eigenvectors $\bm{v}_{1},\ldots,\bm{v}_{k}$ corresponding to the $k$ nonzero eigenvectors of $\bm{\mathcal{L}}$. *

Proof A.5**.**

From the previous lemma we know that the eigenvectors of $\bm{\mathcal{L}}$ can be written as $\bm{v}_{i}=\bm{G}\bm{N}_{g}^{-1/2}\bm{u}_{i}$, where $\bm{u}_{i}$ are the eigenvectors of $\bm{\Gamma}$. Assembling these relationships into matrices we obtain

[TABLE]

Since $\bm{U}$ is an orthogonal matrix and $\bm{N}_{g}^{1/2}$ is diagonal, both matrices are invertible. It then follows that

[TABLE]

Corollary A.6**.**

*Consider the row-normalized eigenvector matrix $\widetilde{\bm{V}}=\operatorname{diag}(\bm{V}\bm{V}^{\top})^{-1/2}\bm{V}$. There exists an orthogonal transformation $\bm{\mathcal{U}}$ such that $\widetilde{\bm{V}}=\bm{G}\bm{\mathcal{U}}$. *

Proof A.7**.**

Note that

[TABLE]

Combining this equation with (40) and the fact that $\bm{G}$ is an indicator matrix it follows that

[TABLE]

*Since $\bm{U}$ corresponds to the eigenvector basis of the symmetric normalized affinity matrix $\bm{\Gamma}$, we know that $\bm{U}$ is orthogonal. *

A.2 Proof of Lemma 2.3

Proof A.8**.**

We define the matrices $\bm{\mathcal{D}}:=\mathbb{E}[\bm{D}]$ and ${\bm{H}:=\bm{\mathcal{D}}^{-1/2}\bm{A}\bm{\mathcal{D}}^{-1/2}}$ for notational convenience. From the triangle inequality we have

[TABLE]

To obtain our desired concentration inequality, we bound the two terms on the right-hand side of (42) separately.

For the first term we note that since $\|\bm{L}\|\leq 1$,

[TABLE]

We can now analyze the term $\|\bm{I}-\bm{\mathcal{D}}^{-1/2}\bm{D}^{1/2}\|$ as follows:

[TABLE]

where we have used the fact that the spectral norm of a diagonal matrix is its maximal diagonal entry. To make use of the above result, we need to obtain a concentration result for the degrees. Recalling that ${\delta}_{i}=\mathbb{E}(d_{i})=\mathcal{D}_{ii}$, we can use Chernoff’s inequality to obtain a bound of the form:

[TABLE]

If we choose $t_{0}=\sqrt{\ln(4n/\epsilon)/\delta_{\text{min}}}$, the above bound holds for all $i$, and by a simple transformation we obtain $\mathbb{P}(|d_{i}/\delta_{i}-1|\geq t_{0})\leq\epsilon/2n$ for all $i$. Using the union bound, we thus have with probability at least $1-\epsilon/2$:

[TABLE]

Therefore, plugging this result into (A.8), with probability $1-\epsilon/2$ it holds that

[TABLE]

For the second term $\|\bm{H}-\bm{\mathcal{L}}\|$, we can use the following matrix concentration inequality.

Lemma A.9** (Bernstein Matrix inequality [9]).**

Let $\bm{X}_{i}$ be independent $n\times n$ dimensional random Hermitian matrices such that ${\|\bm{X}_{i}-\mathbb{E}[\bm{X}_{i}]\|\leq C\quad\forall i}$, and $\bm{X}=\sum_{i}\bm{X}_{i}$. Then for any $a>0$:

[TABLE]

*where $v^{2}=\|\sum_{i}\mathrm{var}(\bm{X}_{i})\|=\|\sum_{i}\mathbb{E}[\bm{X}_{i}^{2}]-\mathbb{E}[\bm{X}_{i}]^{2}\|$. *

To apply this result, note that we can decompose ${[\bm{H}-\bm{\mathcal{L}}]}=\sum_{i\leq j}\bm{X}_{ij}$, where $\bm{X}_{ij}$ are the Hermitian matrices:

[TABLE]

Here, $\bm{E}^{ij}$ is the canonical basis matrix with entries $[\bm{E}^{ij}]_{ij}=1$ and [math] otherwise. As $\mathbb{E}[\bm{X}_{ij}]=0$ we compute $\text{var}(\bm{X}_{ij})=\mathbb{E}[\bm{X}_{ij}^{2}]$ as follows:

[TABLE]

From this we can bound $v^{2}$ as

[TABLE]

where we have used the bound $1/\delta_{\text{min}}\geq 1/\delta_{j}$.

Based on the above calculations, we select ${a=\sqrt{3\ln(4n/\epsilon)/\delta_{\text{min}}}=\sqrt{3}t_{0}}$, which leads to our desired bound for the second term

[TABLE]

where the first inequality follows from the fact that $C\leq 1/\delta_{\min}$ (cf. Lemma A.9), and for the second inequality we used the fact that $\mathcal{M}(\bm{\mathcal{A}})\in\mathcal{M}_{n}(\epsilon)$ to determine a lower bound on $\delta_{\min}$. Finally, by combining the results in (44) and (A.8) we obtain that with probability $1-\epsilon$,

[TABLE]