Network Inference from Consensus Dynamics with Unknown Parameters

TL;DR

This paper develops algorithms to infer the structure of weighted undirected networks from consensus dynamics data with unknown parameters, using spectral analysis and convex optimization, and demonstrates their superior performance through experiments.

Contribution

It introduces novel algorithms that leverage spectral properties and convex optimization to infer network Laplacians under parameter uncertainty, with theoretical guarantees and empirical validation.

Findings

Algorithms outperform existing methods in network recovery accuracy.

Theoretical guarantees support the effectiveness of the proposed algorithms.

Numerical experiments validate the approach on synthetic and real-world networks.

Abstract

We explore the problem of inferring the graph Laplacian of a weighted, undirected network from snapshots of a single or multiple discrete-time consensus dynamics, subject to parameter uncertainty, taking place on the network. Specifically, we consider three problems in which we assume different levels of knowledge about the diffusion rates, observation times, and the input signal power of the dynamics. To solve these underdetermined problems, we propose a set of algorithms that leverage the spectral properties of the observed data and tools from convex optimization. Furthermore, we provide theoretical performance guarantees associated with these algorithms. We complement our theoretical work with numerical experiments, that demonstrate how our proposed methods outperform current state-of-the-art algorithms and showcase their effectiveness in recovering both synthetic and real-world networks.

Equations

$${\mathcal{L}}_c=\{\mathbf{L}|L_{ij}=L_{ji}\le 0\text{ for }i\ne j,L1=0\}.$$

$$\mathbf{x}[t]=\mathbf{x}[t-1]-{\alpha }_t\mathbf{Lx}[t-1]=(\mathbf{I}-{\alpha }_t\mathbf{L})\mathbf{x}[t-1].$$

$$\mathbf{y}=h(\mathbf{L})\mathbf{\xi },\text{where }h(\mathbf{L})=\prod _{t=1}^T(\mathbf{I}-{\alpha }_t\mathbf{L}).$$

$${\mathbf{y}}_k=h_k(\mathbf{L}){\mathbf{\xi }}_k,\text{where }{h}_k(\mathbf{L})=\prod _{t=1}^{T_k}(\mathbf{I}-{\alpha }_t^{(k)}\mathbf{L}),$$

$$\mathbb{E}[{\mathbf{y}}_k{\mathbf{y}}_k^{\top }]=h(\mathbf{L})\mathbb{E}[\mathbf{\xi }{\mathbf{\xi }}^{\top }]h(\mathbf{L}{)}^{\top }={\sigma }^{2}h(\mathbf{L}{)}^2,$$

$${\mathbf{C}}_y=\mathbb{E}[{\mathbf{y}}_k{\mathbf{y}}_k^{\top }]={\sigma }^{2}h(\mathbf{L}{)}^2={\sigma }^2\mathbf{V}h(\mathbf{\Lambda }{)}^2{\mathbf{V}}^{\top }.$$

$${\lambda }_i({\mathbf{C}}_y)={\sigma }^{2}h({\lambda }_i{)}^2={\sigma }^2\prod _{t=1}^T(1-{\alpha }_t{\lambda }_i{)}^2.$$

$${\mathbf{S}}_M=\frac{1}{M}\sum _{k=1}^M{\mathbf{y}}_k{\mathbf{y}}_k^{\top }$$

$${\hat{\sigma }}^2={\lambda }_{\max }({\mathbf{S}}_M)={\lambda }_1({\mathbf{S}}_M).$$

$$\sqrt{\frac{{\lambda }_i({\mathbf{S}}_M)}{{\hat{\sigma }}^2}}=h({\hat{\lambda }}_i)=\prod _{t=1}^T(1-{\alpha }_t{\hat{\lambda }}_i),$$

$$\Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2\le C_{\sigma ,\delta }\sqrt{N/M},$$

$$\Vert {\hat{\mathbf{v}}}_i-{\mathbf{v}}_i{\Vert }_2\le \frac{2^{3/2}\Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2}{\min ({\lambda }_{i-1}({\mathbf{C}}_y)-{\lambda }_i({\mathbf{C}}_y),{\lambda }_i({\mathbf{C}}_y)-{\lambda }_{i+1}({\mathbf{C}}_y))},$$

$$\Vert {\hat{\mathbf{v}}}_i-{\mathbf{v}}_i{\Vert }_2\le \frac{2^{3/2}{C}_{\sigma ,\delta }\sqrt{N/M}}{\min ({\lambda }_{i-1}({\mathbf{C}}_y)-{\lambda }_i({\mathbf{C}}_y),{\lambda }_i({\mathbf{C}}_y)-{\lambda }_{i+1}({\mathbf{C}}_y))}.$$

$$|{\hat{\lambda }}_i-{\lambda }_i|<C_{\sigma ,{\alpha }_t,{\lambda }_i}\Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2,$$

$$|{\hat{\lambda }}_i-{\lambda }_i|<C_{\sigma ,{\alpha }_t,{\lambda }_i,\delta }\sqrt{N/M}.$$

$$|{\lambda }_i({\mathbf{S}}_M)-{\lambda }_i({\mathbf{C}}_y)|\le \Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2$$

$$|\sqrt{{\lambda }_i({\mathbf{S}}_M)}-\sqrt{{\lambda }_i({\mathbf{C}}_y)}|$$

$$\le \frac{\Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2}{\sqrt{{\lambda }_i({\mathbf{C}}_y)}}.$$

$$|h({\hat{\lambda }}_i)-h({\lambda }_i)|\le c_{\sigma ,{\alpha }_t,{\lambda }_i}\Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2,$$

$${\min }_{\lambda \le {\alpha }_{\max }^{-1}}|h^{\prime }(\lambda )|\ge {\sum }_{k=1}^T{\alpha }_k{\prod }_{t\ne k}(1-{\alpha }_t{\alpha }_{\max }^{-1}):=c_{\alpha }.$$

$$|h({\hat{\lambda }}_i)-h({\lambda }_i)|>\underset{\lambda \le {\alpha }_{\max }^{-1}}{\min }|h^{\prime }(\lambda )|\cdot |{\hat{\lambda }}_i-{\lambda }_i|\ge c_{\alpha }|{\hat{\lambda }}_i-{\lambda }_i|.$$

$$|{\hat{\lambda }}_i-{\lambda }_i|<\frac{c_{\sigma ,{\alpha }_t,{\lambda }_i}}{c_{\alpha }}\Vert {\mathbf{S}}_M-{\mathbf{C}}_y{\Vert }_2.$$

$${\mathbf{L}}^{\ast }=\underset{\mathbf{L}}{\mathrm{argmin}}d(\mathbf{L},\hat{\mathbf{L}})+\beta \Vert \mathrm{vec}(\mathbf{L}){\Vert }_1,\text{s.t. }\mathbf{L}\in {\mathcal{L}}_c,$$

$${\mathbf{a}}^{\ast }=\underset{\mathbf{a}}{\mathrm{argmin}}\Vert \mathbf{Pa}-\hat{\mathbf{b}}{\Vert }_2^2+4\beta \Vert \mathbf{a}{\Vert }_1,\text{ s.t. }\mathbf{a}\ge 0.$$

$$\underset{M\to \infty }{\lim }[{\mathbf{V}}^{\top }{\mathbf{S}}_M\mathbf{V}{]}_{ij}=\underset{M\to \infty }{\lim }{B}_{ij}^{(M)}=0.$$

$${\mathbf{B}}^{(M)}=\frac{1}{M}\sum _{k=1}^M{h}_k(\mathbf{\Lambda }){\mathbf{V}}^{\top }{\mathbf{\xi }}_k{\mathbf{\xi }}_k^{\top }\mathbf{V}{h}_k(\mathbf{\Lambda })=\frac{1}{M}\sum _{k=1}^M{\mathbf{w}}_k{\mathbf{w}}_k^{\top },$$

$$B_{ij}^{(M)}=\frac{1}{M}\sum _{k=1}^M[{\mathbf{w}}_k{]}_i[{\mathbf{w}}_k{]}_j.$$

$$\mathbb{P}\left[|B_{ij}^{(M)}|\ge l\right]\le 2\exp \left(-\frac{M{l}^2}{2{\nu }_{\ast }^2}\right)$$

$$\underset{M\to \infty }{\lim }\mathbb{P}\left[|B_{ij}^{(M)}|\ge l\right]\le \underset{M\to \infty }{\lim }2\exp \left(-\frac{M{l}^2}{4{\sigma }^2}\right)=0,$$

$$B_{11}^{(M)}>B_{22}^{(M)}>\cdots >B_{NN}^{(M)},$$

$$B_{ii}^{(M)}=\frac{1}{M}\sum _{k=1}^M([{\mathbf{w}}_k{]}_i{)}^2,$$

Full text

I Introduction

Networks have become a fundamental tool to model systems across Science and Engineering, with applications ranging from physical to socio-economic and biological domains [2, 3, 4]. In certain cases, we may have relational data that quantifies the couplings between the system entities directly. For example, in transportation networks, we can measure traffic flows between different points in space. However, in many instances, the true couplings between the system entities are unknown and have to be inferred from data collected from the system entities. This task, which we refer to as network inference, is thus a fundamental step prior to any further network analysis.

Network inference has been studied from several perspectives in the literature, and different models have been adopted that associate the network topology with the observed data. Historically, there are two main lines of research, which are based on statistical models and physically-motivated models [5, 6], respectively. A well-known statistical model is the graphical model [7], where the network (graph) encodes conditional independence relations among random variables defined on the system entities (nodes). By employing the graphical model, the problem of inferring the network is thus converted to a particular estimation of the joint probability distribution of these random variables [5, 6]. Associated algorithms include the graphical LASSO [8, 9, 10], which incorporates a graph sparsity prior into the maximum likelihood estimator to recover the precision matrix of these random variables. Physically-motivated models assume that the observed data is generated by some physical process on the network such as diffusion [11, 12, 13, 14], and the network recovered is expected to explain the generative process of the observations.

Graph signal processing (GSP) [15, 16, 17], a fast growing research area that seeks to extend concepts and methods in classical digital signal processing to graphs, offers a new perspective to the problem of network inference. Works leveraging GSP tools include: (i) Models based on signal smoothness [18, 19], where methods to infer the network topology from smooth signals are proposed via minimizing a regularized graph Laplacian quadratic form; (ii) Models based on causal dependency [20], where algorithms are put forth to recover the network structure capturing the dependencies among time series; and (iii) Models based on network diffusion [21, 22, 23, 24], which solve the problem of network inference from a stationary graph process by leveraging the spectral information contained in the observations. Extensions to non-stationary diffusion processes have also been recently proposed [25, 26]. For a thorough review on the topic of network inference, see [5, 6].

In this paper, we consider the problem of network inference from snapshot observations of consensus dynamics. Consensus has been one of the most popular and well-studied dynamics on networks [27, 28, 29] due to both its analytic tractability and its simplicity in approximating several fundamental behaviors. For example, in socio-economic domains, consensus provides a model for opinion formation in social networks. For engineering systems, it has been considered as a basic building block for an efficient distributed computation of global functions in networks of sensors, robots, or other agents.

We consider snapshot data since this is a common scenario in modern observational datasets. For instance, in the study of gene-expression via single-cell RNA sequencing [30], we can only obtain snapshot data, as the process of obtaining a sample destroys the cell under consideration. To obtain multiple samples, we thus have to replicate the experiments with similarly prepared cells. Another example is the monitoring of ecological populations from abundance data [31]. While the underlying process may well be continuous in time, very often we can only access snapshot information at specific instances of time.

Contributions: We study the problem of inferring a network topology from snapshots of discrete-time consensus dynamics with parametric uncertainty. Our specific contributions can be summarized as follows: (i) We formulate three problems (Problems 1-3 in Section II) with increasing degree of uncertainty about the parameters of the dynamics, and propose algorithms to solve them; (ii) We provide provable bounds on the performance of key steps within these algorithms; and (iii) We illustrate the performance of the proposed methods and compare it with state-of-the-art solutions for synthetic and real-world settings.

The proposed observation model strikes a balance between being specific and versatile. On the one hand, the advocated model is more specific than those assuming signal smoothness [18, 19] or generic diffusion processes [21, 22]. This enables us to obtain better network recovery performance when the model (approximately) holds. On the other hand, our formulation with unknown filter parameters (i.e., Problem 3) includes a wider range of settings compared to those assuming a single underlying stationary graph process [21, 22, 23] or specific functional forms for the dynamics [24]. For our Problems 1 and 2 which present scenarios with less parameter uncertainty compared to the more general Problem 3, the work [23] is closest to our work. However, the proposed solutions differ significantly to those in [23] as discussed throughout the paper and illustrated numerically in Section VI.

Paper outline: The remainder of this article is structured as follows. Preliminary concepts related to GSP and consensus are reviewed in Section II. In Section III, we give a formal account of our problem setup and introduce three concrete problem formulations. In Section IV we provide a detailed analysis of the first problem along with the algorithm that we propose for its solution. Section V builds upon insights gained by studying the first problem, and provides solutions to the other two problems considered. Numerical experiments based on both synthetic and real-world data are presented in Section VI, and closing remarks are included in Section VII.

Notation: The entries of a matrix $\mathbf{X}$ and a vector $\mathbf{x}$ are denoted by $X_{ij}$ and $x_{i}$, respectively; to avoid confusion, the alternative notation $[\mathbf{X}]_{ij}$ and $[\mathbf{x}]_{i}$ will be used occasionally, when dealing with indexed families of matrices and vectors. Operations $(\cdot)^{\top}$, $(\cdot)^{{\dagger}}$, $\mathbb{E}(\cdot)$ and $\mathbb{P}(\cdot)$ represent transpose, pseudo-inverse, expected value and probability, respectively. $\mathbf{0}$, $\mathbf{1}$ and $\mathbf{I}$ refer to the all-zero vector, the all-one vector, and the identity matrix, where the sizes are clear from context. $\text{diag}(\mathbf{x})$ denotes a diagonal matrix whose $i$th diagonal entry is $x_{i}$. $\mathrm{vec}(\mathbf{X})$ stacks the columns of $\mathbf{X}$ into a single column vector. For any set $\mathcal{S}$, $\mathbf{x}_{\mathcal{S}}$ denotes the vector formed by the entries of $\mathbf{x}$ indexed by $\mathcal{S}$.

II Preliminaries

We briefly introduce basic GSP concepts (Section II-A) as well as the mathematical formulation of discrete-time consensus (Section II-B).

II-A Fundamentals of graph signal processing

Graphs and graph signals. Consider a weighted and undirected graph $\mathcal{G}$ with $N$ nodes, whose structure is encoded by the weighted adjacency matrix $\mathbf{A}\in\mathbb{R}^{N\times N}$. If nodes $i$ and $j$ are connected, the edge weight $A_{ij}=A_{ji}>0$ reflects the strength of the connection. If there is no edge between nodes $i$ and $j$, we have $A_{ij}=A_{ji}=0$.

A graph signal defined on $\mathcal{G}$ can be represented as a vector $\mathbf{x}=[x_{1},\cdots,x_{N}]^{\top}\in\mathbb{R}^{N}$, where $x_{i}\in\mathbb{R}$ denotes a scalar signal value associated with node $i$.

Graph shift operator and graph filters. A graph shift operator $\mathbf{S}\in\mathbb{R}^{N\times N}$ [15], is a matrix whose off-diagonal sparsity pattern is identical with the adjacency matrix: $S_{ij}$ can only be non-zero if $A_{ij}\neq 0$ or $i=j$. Typical choices for $\mathbf{S}$ are the adjacency matrix [15], the graph Laplacian [16] and their respective generalizations.

A graph filter is a map between graph signals and is defined as a matrix function $h(\cdot)$ of a graph shift operator $\mathbf{S}$. An important class of graph filters are linear and shift-invariant (LSI) graph filters. An LSI filter can be expressed as a polynomial of $\mathbf{S}$, i.e. $h(\mathbf{S})=\sum_{l=0}^{T}h_{l}\mathbf{S}^{l}$, where $T$ and $\{h_{l}\}$ denote the filter degree and filter coefficients, respectively [15, 32]. For a given input signal $\mathbf{x}$, the output of the graph filter is given by $\mathbf{y}=h(\mathbf{S})\mathbf{x}$.

The set of combinatorial graph Laplacians. In this paper, we concentrate on the combinatorial graph Laplacian (CGL) as our graph shift operator, defined as $\mathbf{L}=\text{diag}(\mathbf{A1})-\mathbf{A}$. The set of all CGL matrices can be written as

[TABLE]

Since $\mathbf{L}$ is a real and symmetric matrix, its eigendecomposition can be written as $\mathbf{L}=\mathbf{V\Lambda V}^{\top}$, where $\mathbf{V}$ is a unitary matrix whose columns are the eigenvectors of $\mathbf{L}$, and $\mathbf{\Lambda}=\text{diag}(\boldsymbol{\lambda})$ collects the eigenvalues. Notice that (1) implies that $\mathbf{L}$ is diagonally dominant, which ensures that $\mathbf{L}$ is positive semi-definite.

Throughout the paper we assume the eigenvalues of $\mathbf{L}$, $0=\lambda_{1}<\lambda_{2}<\cdots<\lambda_{N}$, are distinct. This assumption is not fundamental from a technical viewpoint, but simplifies the presentation of our results. In particular, it implies that $\mathcal{G}$ is connected.

II-B Discrete-time consensus dynamics

We consider discrete-time linear consensus dynamics [27, 28, 29], evolving on a graph $\mathcal{G}$ with Laplacian $\mathbf{L}$:

[TABLE]

Here the vector $\mathbf{x}[t]$ is a (time-varying) graph signal, whose entries $x_{i}[t]$ correspond to the opinion of agent $i$ in the network at time $t$. We assume that $0<\alpha_{t}<{\lambda_{N}^{-1}}$ for all $t$, such that the dynamics is stable [28].

The above equations describe a dynamics in which agent $i$ updates its opinion according to a linear combination of (i) its previous opinion and (ii) a weighted discrepancy with its neighbors at the previous time point. The parameter $\alpha_{t}$ is called the diffusion rate and describes the weight given to the discrepancy term in the update at time $t$. Under the dynamics (2), the opinions of all agents coincide asymptotically, i.e., $\lim\limits_{t\to\infty}\mathbf{x}[t]=c\mathbf{1}$ where $c$ is a real constant.

III Problem formulation

We study the problem of inferring the topology of a network from the observation of $M$ consensus dynamics at a single point in time. In contrast to other identification tasks, we assume that only such snapshot information is available, i.e., no trajectories of the states or detailed knowledge of the initial condition is available.

Within this setting, we study three different identification problems. The difference between these problems lies in the amount of knowledge that we assume about the dynamics. In our first problem formulation (Problem 1), we assume that the diffusion rates and observation times are known. Subsequently, we relax these assumptions, and allow for unknown – albeit constant – diffusion rates and observation times (Problem 2). Finally, in Problem 3 we tackle the most general case where we observe multiple consensus dynamics with unknown and possibly different diffusion rates and observation times.

From the perspective of GSP, our task can also be expressed as learning the graph Laplacian $\mathbf{L}$ from the observations of a set of filtered graph signals $\mathbf{y}$. Let $\bm{\xi}:=\mathbf{x}[0]$ denote the graph signal at time zero and $\mathbf{y}:=\mathbf{x}[T]$ the state of the dynamics at a specific observation time $T>0$. From (2) we know that $\mathbf{y}$ and $\bm{\xi}$ are related as

[TABLE]

Thus the dynamical description given in (2) can be equivalently phrased in terms of (3) from a GSP perspective. In this paper we adopt this GSP perspective and assume that we observe a set of filtered graph signals $\{\mathbf{y}_{k}\}_{k=1}^{M}$ of the form (3). Before turning to the specific problem formulations we state our main assumptions.

III-A Assumptions

To ensure that there is some non-trivial information about $\mathbf{L}$ contained in the observations $\{\mathbf{y}_{k}\}_{k=1}^{M}$, we assume that the dynamics has not reached asymptotic consensus. Also, as discussed in Section II-B, we assume diffusion rates small enough to ensure convergence.

Assumption 1** **(Finite-time consensus dynamics)

The observation time $T$ is finite, i.e. $T<\infty$, and the diffusion rates satisfy $0<\alpha_{t}<{\lambda_{N}^{-1}}$ for all $t$.

Second, we assume that we cannot control or observe the initial input $\bm{\xi}$, but have some knowledge about its distribution. Specifically, we make the following assumption about the unknown input signal throughout this paper.

Assumption 2** **(White Gaussian input)

The initial condition $\bm{\xi}$ in (3) is a Gaussian random graph signal $\bm{\xi}\sim\mathrm{N}(\mathbf{0},\sigma^{2}\mathbf{I})$, where $\sigma^{2}$ denotes the input power.

III-B Formal problem statements

Consider a set of $M$ independent, identically distributed initial conditions $\{\bm{\xi}_{k}\}_{k=1}^{M}$ with unknown input power $\sigma^{2}$. In our most general formulation, each of these initial conditions evolve according to a consensus dynamics such that

[TABLE]

for $k=1,\cdots,M$. In our first problem, we focus on the specific case where all the dynamics are identical, and the only unknown parameter is the input power $\sigma^{2}$.

Problem 1** **(unknown input noise level)

Given the set of $M$ outputs $\{\mathbf{y}_{k}\}_{k=1}^{M}$ from (4) for an identical dynamics $h_{k}(\mathbf{L})\equiv h(\mathbf{L})$, estimate the Laplacian ${\mathbf{L}}$ for known observation time $T$ and diffusion rates $\{\alpha_{t}\}$.

The second problem formulation may be interpreted as a relaxation of Problem 1 where we do not know the observation time $T$ nor the diffusion rate $\alpha$, which we assume to be identical for all time steps in the context of Problem 2.

Problem 2** **(unknown graph filter)

Given the set of $M$ outputs $\{\mathbf{y}_{k}\}_{k=1}^{M}$ from (4) for an identical dynamics $h_{k}(\mathbf{L})\equiv h(\mathbf{L})$, estimate the Laplacian ${\mathbf{L}}$ for unknown observation time $T$ and constant diffusion rate $\alpha_{t}\equiv\alpha$.

Finally, we consider the more challenging and general setting where there are $M$ different consensus dynamics with unknown parameters evolving on a single network with Laplacian $\mathbf{L}$.

Problem 3** **(unknown set of filters)

Given the set of $M$ outputs $\{\mathbf{y}_{k}\}_{k=1}^{M}$ from (4) for different dynamics $h_{k}(\mathbf{L})$, estimate the Laplacian ${\mathbf{L}}$ for unknown observation times $T_{k}$ and diffusion rates $\{\alpha_{t}^{(k)}\}$.

The three problems above may be interpreted in terms of a hierarchy of assumptions on the knowledge of the parameters. While Problem 1 is already challenging as we only rely on snapshot data and have no knowledge of the exact initial conditions, in Problems 2 and 3 we make even weaker assumptions on the diffusion parameters $T_{k}$ and $\{\alpha^{(k)}_{t}\}$. However, in many applications we are precisely confronted with such strongly underdetermined problems: while we may have some (approximate) model of the functional form of the dynamics (such as a consensus or a diffusion), the specific parameters of such dynamics as well as their initial conditions are often not known. Hence, ideally we would like to infer the network together with the parameters of the dynamics.

To see the practical relevance of our problem setup, let us consider an illustrative example for Problem 3 in the context of social sciences. Assume that we observe, at a specific point in time, the opinion profile of all agents in a social network represented by $\mathbf{L}$ regarding $M$ independent topics, each of which evolved according to a consensus dynamics as described in (4). The discussion about each of the $M$ topics, which we index by $k\in\{1,\cdots,M\}$, may have started at a different point in time – corresponding to unknown durations $T_{k}$. Moreover, the interactions between the agents may have been heterogeneous across topics and time – associated with unknown diffusion rates $\alpha^{(k)}_{t}$. Our goal is to identify the underlying social network topology $\mathbf{L}$ from the observation of $M$ opinion profiles $\{\mathbf{y}_{k}\}$ at a given time point. For further illustration of these problems, including real-world data implementations, see Section VI-D.

Note that in Problems 1 and 2, $\{\mathbf{y}_{k}\}_{k=1}^{M}$ are independent realizations of the output of a single, possibly unknown (Problem 2) graph filter $h(\mathbf{L})$ under a white noise input. This is precisely the definition of a stationary graph process [33, 34, 35]. Hence, Problems 1 and 2 can also be expressed as inferring $\mathbf{L}$ from independent samples of an identical stationary graph process $\mathbf{y}$. This has important implications for the solution of these problems that we can exploit. In contrast, in Problem 3 each observation corresponds to a realization of a different stationary process, and thus requires an adjusted solution strategy.

Remark 1** **(Related work [24])

In [24], the authors consider a set of heat diffusion processes which can be modeled by graph filters in the form of $h(\mathbf{L})=\exp(-\tau\mathbf{L})$ with different nonnegative constants $\tau$. The authors have studied the problem of inferring the graph structure from a set of observations $\mathbf{y}_{k}$, each of which corresponds to a different linear combination (i.e., a mixture) of these heat diffusion processes. In contrast, in our problem settings, each observation $\mathbf{y}_{k}$ corresponds to the output of a single, albeit parametrically uncertain, consensus dynamics. For different observations, however, the corresponding consensus dynamics can be identical (Problems 1 and 2) or different (Problem 3).

Comparing the graph filter for consensus dynamics in (3) with the graph filter representing heat diffusion, it is easy to show that they have the same asymptotic behavior, i.e., $\lim_{T\to\infty}\prod_{t=1}^{T}(\mathbf{I}-\alpha_{t}\mathbf{L})=\lim_{\tau\to\infty}\exp(-\tau\mathbf{L})=\mathbf{11}^{\top}/N$. Indeed, we have $\exp(-\tau\mathbf{L})=\lim_{j\rightarrow\infty}(\mathbf{I}-\tau\mathbf{L}/j)^{j}$, and thus when $T$ and $\tau$ take finite values, the graph filters are markedly different.

In the next two sections we propose different algorithms to solve Problems 1, 2 and 3. Our algorithms are optimization-based and can be considered as two-step procedures. First, we extract spectral information about the CGL $\mathbf{L}$ from realizations of the output signals, i.e., we collect information about the eigenvectors and eigenvalues of $\mathbf{L}$. Second, we construct a convex optimization problem to recover the CGL, leveraging the extracted spectral information and incorporating the CGL constraints as well as a sparsity prior.

IV Topology inference from a single consensus process

The focus of this section is on Problem 1, where we are confronted with a single stationary graph process $h_{k}(\mathbf{L})\equiv h(\mathbf{L})$. Accordingly, it follows from (4) and Assumption 2 that

[TABLE]

where the last equality uses the fact that $\mathbf{L}$ is symmetric.

From (5), we can see that information about $\mathbf{L}$ is encoded in the covariance of the output $\mathbf{y}$. Specifically, as $h(\mathbf{L})=\mathbf{V}h(\mathbf{\Lambda})\mathbf{V}^{\top}$, where $[h(\mathbf{\Lambda})]_{ii}=h(\lambda_{i})$ for all $i$, the covariance of the output can be written as

[TABLE]

The following results are an immediate consequence.

Proposition 1

The output covariance matrix $\mathbf{C}_{y}$ and the graph Laplacian $\mathbf{L}$ share the same set of eigenvectors. Moreover, the eigenvalues of $\mathbf{C}_{y}$ are given by the following transformation of the eigenvalues of $\mathbf{L}$:

[TABLE]

Proof : This follows directly from (6). $\square$

Proposition 2

The input power $\sigma^{2}$ of the initial condition $\bm{\xi}$ and the eigenvalues of the graph Laplacian $\mathbf{L}$ can be directly recovered from the spectrum of the output covariance $\mathbf{C}_{y}$.

Proof : Since by assumption $0=\lambda_{1}<\lambda_{2}<\cdots<\lambda_{N}$, it follows from 1 that $\sigma^{2}=\lambda_{1}(\mathbf{C}_{y})>\lambda_{2}(\mathbf{C}_{y})>\cdots>\lambda_{N}(\mathbf{C}_{y})$. Hence, the input power $\sigma^{2}$ equals the largest eigenvalue of $\mathbf{C}_{y}$.

Furthermore, as $0<\alpha_{t}<{\lambda_{N}^{-1}}$ for all $t$, we know that ${\lambda_{N}}<\alpha_{\max}^{-1}$, where $\alpha_{\max}$ is the maximum value of all $\alpha_{t}$. This implies that $\lambda_{i}<\alpha_{\max}^{-1}$ for all $i$. Defining the scalar function $f(\lambda)=\sigma^{2}{h(\lambda)^{2}}$, it can be seen that $f(\lambda)$ is a monotonically decreasing function of $\lambda$ for $\lambda\leq\alpha_{\max}^{-1}$. Hence, the eigenvalues of $\mathbf{L}$ can be uniquely recovered from those of $\mathbf{C}_{y}$ using (7). $\square$

Having established that the output covariance $\mathbf{C}_{y}$ captures all essential information about $\mathbf{L}$, we now discuss how well we can approximate $\mathbf{C}_{y}$ from a finite number of samples. Given a set of $M$ independent samples $\{\mathbf{y}_{k}\}_{k=1}^{M}$ from an identical distribution with zero-mean and covariance $\mathbf{C}_{y}$, we can compute the sample covariance matrix

[TABLE]

as an estimate of the covariance $\mathbf{C}_{y}$.

The sample covariance $\mathbf{S}_{M}$ is a real symmetric and positive semi-definite matrix, and we can write its eigendecomposition as $\mathbf{S}_{M}=\mathbf{U}\mathbf{\Sigma}\mathbf{U}^{\top}$, where the diagonal matrix $\mathbf{\Sigma}$ collects the eigenvalues of $\mathbf{S}_{M}$ denoted by $\lambda_{1}(\mathbf{S}_{M})\geq\lambda_{2}(\mathbf{S}_{M})\geq\cdots\geq\lambda_{N}(\mathbf{S}_{M})\geq 0$. The input power $\sigma^{2}$ can be estimated as

[TABLE]

We estimate the eigenvectors of $\mathbf{L}$ as $\hat{\mathbf{V}}=\mathbf{U}$ and the eigenvalues of $\mathbf{L}$ by replacing $\mathbf{C}_{y}$ with $\mathbf{S}_{M}$ in (7) to obtain

[TABLE]

where $\hat{\lambda}_{i}$ denotes the estimate of $\lambda_{i}$. Notice that in (10) we leverage the fact that $h(\lambda)$ is a non-negative function of $\lambda$ when $\lambda\leq\alpha_{\max}^{-1}$. It follows from (9) that $\sqrt{\lambda_{i}(\mathbf{S}_{M})/\hat{\sigma}^{2}}\in[0,1]$, and $\hat{\lambda}_{i}$ is the unique real root in the range $[0,\alpha_{\max}^{-1}]$ of the polynomial in (10). After obtaining $\hat{\mathbf{V}}$ and $\hat{\boldsymbol{\lambda}}=[\hat{\lambda}_{1},\cdots,\hat{\lambda}_{N}]^{\top}$, we can estimate $\mathbf{L}$ as $\hat{\mathbf{L}}=\hat{\mathbf{V}}\text{diag}(\hat{\boldsymbol{\lambda}})\hat{\mathbf{V}}^{\top}$. We refer to this method as InverseFilter, and it is summarized in Algorithm 1. The same idea has been used in the algorithm proposed in [23].

The theoretical performance guarantees of the InverseFilter approach are stated in Theorems 1 and 2. Both theorems rely on the following instrumental lemma.

Lemma 1

Consider $\mathbf{y}_{k}$ generated as in (4) for $1\leq k\leq M$ and assume that $M\geq N$. Then for every $\delta>0$, with probability at least $1-\delta$ one has that

[TABLE]

where $C_{\sigma,\delta}$ is a constant which depends on $\sigma$ and $\delta$.

Proof : Define a vector $\mathbf{z}\in\mathbb{R}^{N}$ satisfying $\|\mathbf{z}\|_{2}=1$. Since $\mathbf{y}_{k}$ is a zero-mean Gaussian random vector, we have that $\mathbb{E}[\mathbf{z}^{\top}\mathbf{y}_{k}]=0$ and $\mathbb{E}[e^{\mathbf{z}^{\top}\mathbf{y}_{k}}]=e^{\frac{1}{2}\mathbf{z}^{\top}\mathbf{C}_{y}\mathbf{z}}\leq e^{\frac{\sigma^{2}}{2}}$. Hence, $\mathbf{z}^{\top}\mathbf{y}_{k}$ is sub-Gaussian satisfying $\mathbb{P}\left[|\mathbf{z}^{\top}\mathbf{y}_{k}|>t\right]<2e^{-\frac{t^{2}}{2\sigma^{2}}}$ 111The concepts of sub-Gaussian random variable and its tail bounds are reviewed in Appendix A; see Definition 1 and Lemma 2.. The result follows from Proposition 2.1 in [36]. $\square$

Theorem 1** **(Eigenvectors of InverseFilter)

For all $i$, the eigenvectors $\hat{\mathbf{v}}_{i}$ estimated by InverseFilter satisfy

[TABLE]

where we set $\lambda_{0}(\mathbf{C}_{y})=\infty$ and $\lambda_{N+1}(\mathbf{C}_{y})=-\infty$. Moreover, if $M\geq N$ we have that for every $\delta>0$, with probability at least $1-\delta$,

[TABLE]

Proof : The first result is a restatement of the Davis-Kahan theorem that follows from Corollary 1 in [37]. The second result is obtained by combining the first result with the bound in Lemma 1. $\square$

Notice that the estimated eigenvectors $\hat{\mathbf{v}}_{i}$ and the true ones ${\mathbf{v}}_{i}$ are defined up to a sign inversion. Thus, for the bounds in Theorem 1 to hold both signs have to be picked in a consistent manner. More precisely, we set the sign of our estimated eigenvectors such that $\hat{\mathbf{v}}_{i}^{\top}\mathbf{v}_{i}\geq 0$ for all $i$.

The estimation error on the eigenvalues obtained via InverseFilter can also be bounded, as we show next.

Theorem 2** **(Eigenvalues of InverseFilter)

Assuming that $\sigma$ is known, one has that

[TABLE]

where $1\leq i\leq N$, $\hat{\lambda}_{i}$ is obtained using (10) ($\hat{\sigma}$ is set as the ground truth $\sigma$ here), and $C_{\sigma,\alpha_{t},\lambda_{i}}$ is a constant which depends on $\sigma$, $\{\alpha_{t}\}_{t=1}^{T}$ and $\lambda_{i}$. Moreover, if $M\geq N$ we have that for every $\delta>0$, with probability at least $1-\delta$

[TABLE]

Proof : According to Theorem 6.4.3 in [38], we have that

[TABLE]

for $1\leq i\leq N$, which implies that

[TABLE]

It follows from (7), (10) and (IV) that

[TABLE]

where $c_{\sigma,\alpha_{t},\lambda_{i}}=\frac{1}{\sigma\sqrt{\lambda_{i}(\mathbf{C}_{y})}}$ is a constant depending on $\sigma$, $\{\alpha_{t}\}_{t=1}^{T}$ and $\lambda_{i}$. The derivative of the filter function $h(\lambda)$ is $h^{\prime}(\lambda)=-\sum_{k=1}^{T}\alpha_{k}\prod_{t\neq k}(1-\alpha_{t}\lambda)$, and this implies that

[TABLE]

Since both of $\lambda_{i}$ and $\hat{\lambda}_{i}$ are no more than $\alpha_{\max}^{-1}$, we have

[TABLE]

Combining this with (14), we have that

[TABLE]

Set $C_{\sigma,\alpha_{t},\lambda_{i}}={c_{\sigma,\alpha_{t},\lambda_{i}}}/{c_{\alpha}}$, and the proof of (11) is completed. Finally, (12) follows by combining (11) with the result in Lemma 1. $\square$

Theorems 1 and 2 reveal that, for the InverseFilter method, the estimation error bounds of the eigenvalues and eigenvectors of $\mathbf{L}$ are both proportional to $\|\mathbf{S}_{M}-\mathbf{C}_{y}\|_{2}$, i.e., the estimation error of the covariance $\mathbf{C}_{y}$. Moreover, from Lemma 1 it follows that $\|\mathbf{S}_{M}-\mathbf{C}_{y}\|_{2}$ is bounded by $\sqrt{N/M}$ up to a scalar multiple when $M\!\geq\!N$. Hence, under the same assumption, the estimation errors of the eigenvalues and eigenvectors of $\mathbf{L}$ obtained by InverseFilter decrease as $1/\sqrt{M}$ with the number of observations. This implies that InverseFilter is a consistent estimator of the true Laplacian ${\mathbf{L}}$. In addition, the estimation errors of the eigenvectors depend on the gaps between the eigenvalues of the covariance matrix. More precisely, if a given eigenvalue has a large gap with the rest, its corresponding eigenvector can be estimated with high accuracy.

IV-A Solution to Problem 1: Leverage Laplacian structure

One drawback of the InverseFilter method is that the estimate $\hat{\mathbf{L}}$ obtained is generally not a valid CGL. This is especially a problem when the sample size is small, significantly deteriorating the estimation accuracy. To alleviate this problem we enforce the constraint that our estimate should (i) be a valid Laplacian belonging to the set $\mathcal{L}_{c}$ in (1), and (ii) correspond to a sparse graph. To this end, we solve the following convex optimization problem,

[TABLE]

where $d(\cdot,\cdot)$ is a convex function which can reflect how close two matrices are, such as a Bregman divergence [39].

In this paper, we consider common choices for $d(\mathbf{L},\hat{\mathbf{L}})$ such as $\|\mathbf{L}-\hat{\mathbf{L}}\|_{\mathrm{F}}$ and $\|\mathbf{L}-\hat{\mathbf{L}}\|_{2}$ as well as their squares. In practical problems, the underlying graph is usually sparse, hence we add a regularization term $\beta\|\mathrm{vec}(\mathbf{L})\|_{1}$ to promote the graph sparsity, where $\beta$ is a nonnegative regularization parameter. The $\ell_{1}$ norm is used as a convex surrogate of the non-convex $\ell_{0}$ pseudo-norm (i.e., the number of non-zero entries in a vector).

The constraint $\mathbf{L}\in\mathcal{L}_{c}$ in (15) leverages the structural information and guarantees that the estimate ${\mathbf{L}}^{*}$ is a valid CGL. Note that for $\beta=0$, problem (15) can be interpreted as finding the nearest CGL to $\hat{\mathbf{L}}$. Hence, we refer to the method recovering $\mathbf{L}^{*}$ (for general $\beta$) as NearestCGL, and summarize it in Algorithm 2. Finally, notice that although we discuss our problem based on the consensus model in (4), the proposed NearestCGL as well as InverseFilter can also be applied to other graph filters whose corresponding functions are nonnegative and one-to-one.

Such functions guarantee that the eigenvalues of $\mathbf{L}$ can be uniquely identified from the eigenvalues of the covariance $\mathbf{C}_{y}$ [cf. (6)].

The proposed algorithm can scale to larger networks by translating the optimization problem (15) into a well-studied canonical formulation for which efficient solvers exist. To this end, let us select $d(\mathbf{L},\hat{\mathbf{L}})=\|\mathbf{L}-\hat{\mathbf{L}}\|_{\mathrm{F}}^{2}$. For convenience, we define the vector $\boldsymbol{\ell}:=\mathrm{vec}(\mathbf{L})$ and similarly $\hat{\boldsymbol{\ell}}:=\mathrm{vec}(\hat{\mathbf{L}})$. We further define two index sets $\mathcal{I}=\{(i-1)N+i\,\,|\,\,1\leq i\leq N\}$ and $\mathcal{J}=\{(i-1)N+j\,\,|\,\,1\leq i<j\leq N\}$. Then $\boldsymbol{\ell}_{\mathcal{I}}$ and $\boldsymbol{\ell}_{\mathcal{J}}$ denote the vectors collecting the diagonal entries and the entries in the lower triangular portion of $\mathbf{L}$, respectively. We also define a nonnegative vector $\mathbf{a}:=-\boldsymbol{\ell}_{\mathcal{J}}$ which contains the weights of all possible edges as well as $\mathbf{b}:=[\boldsymbol{\ell}_{\mathcal{I}}^{\top},\sqrt{2}\boldsymbol{\ell}_{\mathcal{J}}^{\top}]^{\top}$ and $\hat{\mathbf{b}}:=[\hat{\boldsymbol{\ell}}_{\mathcal{I}}^{\top},\sqrt{2}\hat{\boldsymbol{\ell}}_{\mathcal{J}}^{\top}]^{\top}$. Due to the constraint $\mathbf{L1}=\mathbf{0}$, it is easy to construct a (sparse) matrix $\mathbf{P}$ such that $\mathbf{Pa}=\mathbf{b}$. It follows that $\|\mathbf{L}-\hat{\mathbf{L}}\|_{\mathrm{F}}^{2}=\|\mathbf{b}-\hat{\mathbf{b}}\|_{2}^{2}=\|\mathbf{Pa}-\hat{\mathbf{b}}\|_{2}^{2}$ and $\|\mathrm{vec}(\mathbf{L})\|_{1}=4\|\mathbf{a}\|_{1}$. Hence, problem (15) is equivalent to

[TABLE]

Problem (16) is in the form of the $\ell_{1}$-regularized least squares problem with nonnegative constraints. Efficient methods to solve such problems have been well studied. In this paper (Section VI-E), we adopt a solver222The code can be found in https://web.stanford.edu/~boyd/l1_ls/. which uses the truncated Newton interior-point method described in [40].

Remark 2** **(Solution to Problem 2)

A potential (naive) solution to Problem 2 can be obtained by running multiple instances of the NearestCGL method. More precisely, without loss of generality we may absorb the unknown (but constant) $\alpha$ into ${\mathbf{L}}$, thus having only $T$ as an explicit unknown parameter. We may then run NearestCGL for $T\in\{1,2,\cdots,T_{\max}\}$ for some maximum possible observation time $T_{\max}$, to obtain a series of potential Laplacians ${\mathbf{L}}^{*}_{(T)}$. We can then choose the best Laplacian among those by selecting the sparsest one or the one that minimizes $d({\mathbf{L}}^{*}_{(T)},{\hat{\mathbf{L}}_{(T)}})$. For the latter case, the intuition is that $\hat{\mathbf{L}}_{(T)}$ is expected to be far from a valid graph Laplacian if a wrong value of $T$ is selected.

Remark 3** **(Related work [23])

A related method has been proposed in [23] that studies a similar problem with graph filters of the form $h_{\tau}(\mathbf{L})=(\mathbf{L}^{{\dagger}})^{\tau}$ where $\tau\in\mathbb{N}$ and $\tau\geq 1$, and the goal is to jointly estimate $\mathbf{L}$ and $\tau$. The authors propose to perform a line search for $\hat{\tau}=1,2,3,\cdots$, and execute their algorithm (as will be explicitly stated in Section VI-A) for each value of $\hat{\tau}$ to estimate the graph Laplacian. Denote by $\mathbf{L}_{\hat{\tau}}$ the estimate associated with $\hat{\tau}$. For all pairs of $\{\hat{\tau},\mathbf{L}_{\hat{\tau}}\}$, the minimizer of $\|h_{\hat{\tau}}(\mathbf{L}_{\hat{\tau}})^{2}-\mathbf{S}_{M}\|_{\mathrm{F}}$ is then chosen as the solution (assuming unit input power).

The minimization of $\|h_{\hat{\tau}}(\mathbf{L}_{\hat{\tau}})^{2}-\mathbf{S}_{M}\|_{\mathrm{F}}$ behaves similar to minimizing $\|\mathbf{L}_{\hat{\tau}}-h_{\hat{\tau}}^{-1}(\sqrt{\mathbf{S}_{M}})\|_{\mathrm{F}}$. As described for the naive solution above, $\mathbf{L}^{\tau/\hat{\tau}}$ is not a valid graph Laplacian, in general, when $\hat{\tau}\neq\tau$. Hence, the distance between $\mathbf{L}_{\hat{\tau}}$ (which is a valid graph Laplacian guaranteed by the recovery algorithm) and $h_{\hat{\tau}}^{-1}(\sqrt{\mathbf{S}_{M}})$ (an estimate of $\mathbf{L}^{\tau/\hat{\tau}}$) is expected to be large for $\hat{\tau}\neq\tau$. However, when the sample size is so small that $h_{\hat{\tau}}^{-1}(\sqrt{\mathbf{S}_{M}})$ cannot approximate $\mathbf{L}^{\tau/\hat{\tau}}$ well, $h_{\hat{\tau}}^{-1}(\sqrt{\mathbf{S}_{M}})$ might be close to some graph Laplacian, which is not necessarily the true $\mathbf{L}$. This might lead to a small $\|h_{\hat{\tau}}(\mathbf{L}_{\hat{\tau}})^{2}-\mathbf{S}_{M}\|_{\mathrm{F}}$ and, thus, yield a worse solution.

Due to this issue, we propose a different approach to Problem 2. Instead of estimating the graph Laplacian and the filter parameter at the same time, we adopt a two-step procedure in which we first estimate the unknown parameter and then estimate the graph Laplacian with the parameter estimate. The proposed method is empirically better than the naive one previously described. Since the proposed method relies on our solution to Problem 3, in the next section we tackle Problem 3 and defer the solution to Problem 2 to Section V-B.

V Topology inference from multiple consensus processes

The main focus of this section is on Problem 3, where our goal is to infer the network topology $\mathbf{L}\in\mathcal{L}_{c}$ from a set of observations $\{\mathbf{y}_{k}\}_{k=1}^{M}$ each of which corresponds to a different consensus dynamics. Different $\mathbf{y}_{k}$ have different covariances $\sigma^{2}{h_{k}(\mathbf{L})^{2}}$, thus in general the sample covariance $\mathbf{S}_{M}$ defined in (8) will not converge to the covariance of any specific $\mathbf{y}_{k}$ for increasing sample size $M$. As a result, the method proposed in Section IV is not applicable for this setting. However, $\mathbf{S}_{M}$ still contains spectral information about the unknown $\mathbf{L}$ and can be leveraged in the estimation. We first discuss how information about the eigenvectors of $\mathbf{L}$ can be inferred from $\mathbf{S}_{M}$. The result is given in Theorem 3, which states that $\mathbf{S}_{M}$ in (8) and $\mathbf{L}$ are simultaneously diagonalizable, provided that the sample size $M$ is sufficiently large. For notational purposes, we define the matrix $\mathbf{B}^{(M)}=\mathbf{V}^{\top}\mathbf{S}_{M}\mathbf{V}$. The proofs in this section use some properties of random variables which are summarized in Appendix A; see Lemmas 3-6.

Theorem 3** **(Eigenvectors of the sample covariance)

For $M\!\to\!\infty$, the eigenbasis $\mathbf{V}$ diagonalizes $\mathbf{S}_{M}$, i.e., for all $i\neq j$

[TABLE]

Proof : The matrix $\mathbf{B}^{(M)}$ can be rewritten as

[TABLE]

where we set $\mathbf{w}_{k}=h_{k}(\mathbf{\Lambda})\mathbf{V}^{\top}\bm{\xi}_{k}$. For each off-diagonal entry in $\mathbf{B}^{(M)}$ we thus have

[TABLE]

Since $\bm{\xi}_{k}\sim\mathrm{N}(\mathbf{0},\sigma^{2}\mathbf{I})$, it follows from the definition of $\mathbf{w}_{k}$ that $\mathbf{w}_{k}\sim\mathrm{N}(\mathbf{0},\sigma^{2}{h_{k}(\mathbf{\Lambda})^{2}})$ and thus $[\mathbf{w}_{k}]_{i}\sim\mathrm{N}(0,\sigma^{2}{h_{k}(\lambda_{i})^{2}})$ for all $i$ are independent.

According to Lemma 4, $[\mathbf{w}_{k}]_{i}[\mathbf{w}_{k}]_{j}$ is therefore a sub-exponential random variable with parameters ${(\nu_{k},b_{k})=(\sqrt{2}\sigma^{2}h_{k}(\lambda_{i})h_{k}(\lambda_{j}),\sqrt{2}\sigma^{2}h_{k}(\lambda_{i})h_{k}(\lambda_{j}))}$. Hence, $MB^{(M)}_{ij}$ is the sum of $M$ independent zero-mean sub-exponential random variables.

It follows from Lemma 6 that $MB^{(M)}_{ij}$ is sub-exponential with parameters $(\sqrt{M}\nu_{\ast},b_{\ast})$ where $\nu_{\ast}=\sqrt{\sum_{k=1}^{M}\nu_{k}^{2}/M}$ and $b_{\ast}=\max_{k}b_{k}$. Then from Lemma 3, we have

[TABLE]

for $0\leq l\leq\nu_{\ast}^{2}/b_{\ast}$. Recall that each $h_{k}(\lambda_{i})$ is upper bounded by $1$, hence $\nu_{\ast}^{2}\leq 2\sigma^{4}$. Consequently, for small enough $l>0$:

[TABLE]

which completes the proof. $\square$

Theorem 3 guarantees that the eigenbasis $\mathbf{V}$ can be recovered by performing the eigendecomposition of $\mathbf{S}_{M}$ for large enough $M$. While in most practical instances of network inference $M$ will be bounded, Theorem 3 can nevertheless be used as a basis for an inference algorithm even if only a finite number of observations are available; see Section V-A.

We remark that the validity of (17) does not imply that $\lim_{M\to\infty}\mathbf{S}_{M}$ exists. Indeed, our weak assumptions on the diffusion parameters $T_{k}$ and $\alpha^{(k)}_{t}$ – which translate into mild conditions on $h_{k}(\mathbf{L})$ – could result in an $\mathbf{S}_{M}$ that does not converge for increasing $M$. However, even if $\mathbf{S}_{M}$ does not converge to a specific matrix, we may interpret (17) as stating that $\mathbf{S}_{M}$ converges to the set of matrices diagonalized by $\mathbf{V}$. As shown in Theorem 4, $\mathbf{S}_{M}$ also provides insights about the eigenvalues of the unknown $\mathbf{L}$.

Theorem 4** **(Eigenvalues of the sample covariance)

Assume that the eigenvalues of $\mathbf{L}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\top}$ are sorted in an increasing order, i.e., $\lambda_{1}<\lambda_{2}<\cdots<\lambda_{N}$. For every $\delta>0$ there exists a large enough number of observations $M_{\delta}$ such that the diagonal entries of the matrix $\mathbf{B}^{(M)}=\mathbf{V}^{\top}\mathbf{S}_{M}\mathbf{V}$ follow an inverse order of the eigenvalues of $\mathbf{L}$, i.e.,

[TABLE]

with probability at least $1-\delta$ for every $M\geq M_{\delta}$.

Proof : It follows from (18) that

[TABLE]

where $[\mathbf{w}_{k}]_{i}\sim\mathrm{N}(0,\sigma^{2}{h_{k}(\lambda_{i})^{2}})$. According to Lemma 5, $([\mathbf{w}_{k}]_{i})^{2}$ is sub-exponential with mean $\sigma^{2}{h_{k}(\lambda_{i})^{2}}$ and parameters $(\nu_{k}^{i},b_{k}^{i})=(2\sigma^{2}{h_{k}(\lambda_{i})^{2}},4\sigma^{2}{h_{k}(\lambda_{i})^{2}})$. Thus, we have

[TABLE]

and $MB^{(M)}_{ii}$ is the sum of $M$ independent sub-exponential random variables. It follows from Lemma 6 that $MB^{(M)}_{ii}$ is sub-exponential with parameters $(\sqrt{M}\nu_{\ast i},b_{\ast i})$ where $\nu_{\ast i}=\sqrt{\sum_{k=1}^{M}(\nu_{k}^{i})^{2}/M}$ and $b_{\ast i}=\max_{k}b_{k}^{i}$. Then from Lemma 3, we have

[TABLE]

for $0\leq l\leq\nu_{\ast i}^{2}/b_{\ast i}$. The last inequality follows from the fact that each $h_{k}(\lambda_{i})$ is upper bounded by $1$, which results that $\nu_{\ast i}^{2}\leq 4\sigma^{4}$. A direct application of the union bound on (21) yields

[TABLE]

for $0\leq l\leq\nu_{\ast}^{2}/b_{\ast}:=\min_{i}(\nu_{\ast i}^{2}/b_{\ast i})$, from which it immediately follows that

[TABLE]

where we fixed a desired probability level at $1-\delta$. Our goal now is to choose $l$ small enough to ensure that (19) is satisfied and then solve for the corresponding number of observations $M_{\delta}$ in (22) using such an $l$. To do this, first recall that $1=h_{k}(\lambda_{1})>\cdots>h_{k}(\lambda_{N})$. We further assume that $h_{k}(\lambda_{i})>h_{k}(\lambda_{j})+\tau$ when $i<j$ for some $\tau>0$, where $\tau$ does not depend on $M$. It then follows from (20) that $e_{i}>e_{j}+\sigma^{2}\tau^{2}$ for $i<j$. Consider a deviation from the mean $l^{\ast}:=\sigma^{2}\tau^{2}/\gamma$ where $\gamma\geq 2$ is large enough to ensure that $l^{\ast}\leq\nu_{\ast}^{2}/b_{\ast}$. By specializing (22) to $l=l^{\ast}$ and solving for $M$ as a function of $\delta$, we have that for all $M$ such that

[TABLE]

every random variable $B^{(M)}_{ii}$ is at most a distance $l^{\ast}$ from its mean with probability at least $1-\delta$. Since by definition $l^{\ast}<(e_{i}-e_{j})/2$ for $i<j$, the variables $B^{(M)}_{ii}$ are sorted in the same order as their means with high probability, and the proof is completed. $\square$

Theorem 4 reveals that, no matter whether $\mathbf{S}_{M}$ converges to a specific matrix or not, the diagonal entries of $\mathbf{V}^{\top}\mathbf{S}_{M}\mathbf{V}$ follow a specific order, i.e., the inverse order of the eigenvalues of the true CGL $\mathbf{L}$ with high probability. This observation, in combination with Theorem 3, is leveraged in Section V-A to develop a network topology inference algorithm for finite $M$.

V-A Solution to Problem 3: Leverage spectral ordering

The recovery of the CGL $\mathbf{L}$ under the setting of Problem 3 is generally an underdetermined problem. Even when we fix the true eigenbasis $\mathbf{V}$ and the ordering of the eigenvalues, there still exists freedom in choosing the exact eigenvalues as long as the order is preserved. To sort out this ambiguity, which amounts to selecting the eigenvalues following a specific order, we assume that the network topology of interest is optimal in some sense. Following the criterion considered in NearestCGL to solve Problem 1 [cf. (15)], in this paper we focus on the optimality based on sparsity. Notice that alternative criteria can be introduced in the form of a generic convex function by replacing the objective in (23a).

Recall that we denote by $\mathbf{S}_{M}=\mathbf{U}\mathbf{\Sigma}\mathbf{U}^{\top}$ the eigendecomposition of the sample covariance defined in (8) where the eigenvalues in $\mathbf{\Sigma}$ are sorted in decreasing order. Our inferred CGL ${\mathbf{L}}^{*}$ can be found by solving the following convex optimization problem.

[TABLE]

Analogously, to the discussion after (15), the $\ell_{1}$ norm is a convex relaxation of the $\ell_{0}$ pseudo-norm, thus the objective (23a) promotes sparsity in the estimate ${\mathbf{L}}^{\ast}$. Alternatively, $\|\mathrm{vec}(\mathbf{L})\|_{1}$ can be replaced by its iterative reweighted counterpart [41], which has shown to perform better in practice.

The constraint in (23b) forces ${\mathbf{L}}^{*}$ to be a valid CGL. The constraints in (23c) impose that ${\mathbf{L}}^{*}$ must be close to being diagonalized by $\mathbf{U}$, i.e. the eigenbasis of the sample covariance $\mathbf{S}_{M}$. It has been shown in Theorem 3 that $\mathbf{U}$ coincides with $\mathbf{V}$, i.e. the eigenbasis of $\mathbf{L}$, for arbitrarily large sample size $M$. However, for all practical implementations, $M$ is finite and thus, we do not require ${\mathbf{L}}^{*}$ to be diagonalized by $\mathbf{U}$ directly. Rather, we require ${\mathbf{L}}^{*}$ to be close to another matrix which is diagonalized by $\mathbf{U}$. The matrix distance is measured by the convex function $d(\cdot,\cdot)$. The problem can be reduced to a linear program if we choose $d(\cdot,\cdot)$ as the maximum norm. Lastly, the constraints in (23d) leverage Theorem 4 and incorporate the fact that the eigenvalues of $\mathbf{S}_{M}$ and the true CGL $\mathbf{L}$ are inversely ordered by forcing the entries in $\boldsymbol{\gamma}$ to follow an approximately increasing order. The parameter $\eta$ is a positive integer which can be adjusted according to the sample size $M$. When $\eta=1$ we impose a strict order on the eigenvalues for the cases in which $M$ is sufficiently large whereas for $\eta>1$ we impose a laxer order for the cases where $M$ is not large enough. We set $\gamma_{N}=1$ to avoid the trivial solution ${\mathbf{L}}^{*}=\mathbf{0}$. We can also set $\gamma_{N}$ to other positive values since it will only vary the scale of the estimate ${\mathbf{L}}^{*}$. Notice that this scale ambiguity is insurmountable given that a common factor across all unknown diffusion rates $\alpha^{(k)}_{t}$ can be absorbed into the unknown $\mathbf{L}$. We denote the proposed method as OrderedSpecTemp given that it uses sorted eigenvectors (spectral templates) in the recovery of ${\mathbf{L}}^{*}$, and we summarize it in Algorithm 3.

Notice that in (23) we may move the constraint (23c) to the objective function and utilize the $\ell_{1}$ norm as the regularization term. The objective function then becomes $d(\mathbf{L},\mathbf{U}\mathrm{diag}(\boldsymbol{\gamma})\mathbf{U}^{\top})+\beta\|\mathrm{vec}(\mathbf{L})\|_{1}$, where the graph sparsity can be tuned with the regularization parameter $\beta$ in addition to the spectral fitting term. As is the case for existing work [18, 19, 23, 24], the accuracy of prior knowledge, i.e., the graph sparsity level, will influence the quality of the selected regularization parameter and thus affect the algorithm performance.

It should also be noticed that in Algorithm 3 we are able to recover $\mathbf{L}$ without the need of recovering the set of parameters $\alpha_{t}^{(k)}$. Indeed, under the setting of Problem 3, $\alpha_{t}^{(k)}$ can be different over $k$ and different over $t$. This weak assumption makes it challenging (if not impossible) to identify $\{\alpha_{t}^{(k)}\}$. Even for the simplified case in which $\alpha_{t}^{(k)}$ is identical for every $k$ but potentially varying with $t$ (i.e., we consider multiple snapshots of an identical linear time-varying dynamics), jointly estimating $\mathbf{L}$ and $\{\alpha_{t}\}$ is challenging. We leave this research direction for future work.

We want to highlight that, while Algorithm 3 is proposed under the setting of Problem 3, it can also be applied to unknown filters of the form $h_{k}(\mathbf{L})\equiv h(\mathbf{L})$. The only requirement is that the filter $h(\lambda)$ is a nonnegative and decreasing function of $\lambda$, which corresponds to a certain type of low-pass graph filter. Under this condition, the eigenvalues of the covariance $\mathbf{C}_{y}$ follow the inverse order of the eigenvalues of $\mathbf{L}$ (if we fix the order of the eigenvectors). In this case, the constraint incorporating the eigenvalue order information is still valid. This constraint is not as restrictive as it might look at first sight. For example, heat diffusion processes [42, 24] satisfy this condition.

Finally, let us outline a way of solving problem (23) efficiently. We select $d(\cdot,\cdot)$ as the squared Frobenius norm. Note that we have $\|\mathbf{L}-\mathbf{U}\text{diag}(\boldsymbol{\gamma})\mathbf{U}^{\top}\|_{\mathrm{F}}^{2}=\|\mathrm{vec}(\mathbf{L})-(\mathbf{U}\odot\mathbf{U})\boldsymbol{\gamma}\|_{2}^{2}$ where $\odot$ denotes the Khatri-Rao product. Hence problem (23) can also be formulated as

[TABLE]

where we have moved the constraint (23c) to the objective function and used the $\ell_{1}$ norm as the regularizer. We can solve the convex problem (24) via block coordinate descent by alternately optimizing $\mathbf{L}$ and $\boldsymbol{\gamma}$. If we fix $\boldsymbol{\gamma}$ and optimize $\mathbf{L}$, the problem will be the same as (16). If we fix $\mathbf{L}$ and optimize $\boldsymbol{\gamma}$, the problem is in fact in the form of a least-squares minimization with linear inequality constraints and, again, standard solver packages can be used.

Remark 4** **(Related work [21])

Leveraging the ordering information of the eigenvectors – in terms of the associated eigenvalues – is essential to the performance of OrderedSpecTemp. Indeed, in [21], a method is proposed based on the spectral templates without the order information to recover adjacency and normalized Laplacian matrices from diffused signals. Such a method cannot be directly extended to our case in order to recover CGL matrices. Notice that we consider a different graph shift operator compared to [21]. If we do not leverage the order information, i.e. ignoring (23d), and add one more constraint such as $\text{tr}(\mathbf{L})\geq 1$ to problem (23) to avoid the solution ${\mathbf{L}}^{*}=\mathbf{0}$, and consider the infinite sample size case in which we have exact spectral templates, then the solution to (23) will be uninformative. More precisely, we would recover ${L}^{*}_{ii}=\frac{1}{N}$ and ${L}^{*}_{ij}=-\frac{1}{N(N-1)}$ for $i\neq j$. In addition, [21] considers a stationary graph process in which multiple observations are obtained from an identical dynamics, but in Problem 3 we consider a more practical and challenging case where observations are sampled from different dynamics.

V-B Solution to Problem 2: Combine NearestCGL and OrderedSpecTemp

Thus far we have proposed solutions to Problem 1 (namely, NearestCGL) and Problem 3 (namely, OrderedSpecTemp). Based on these, we now propose a solution to Problem 2 that empirically outperforms the one delineated in Remark 2.

First notice that for Problem 2, $\mathbf{L}$ can only be recovered up to a scalar multiple since $\alpha$ is unknown and $\alpha\mathbf{L}=\frac{\alpha}{\alpha_{0}}(\alpha_{0}\mathbf{L})$ for any non-zero scalar $\alpha_{0}$. Since $T$ is unknown, NearestCGL as stated in Algorithm 2 cannot be applied. More precisely, (10) cannot be solved to find the estimated eigenvalues in InverseFilter. However, we may still apply OrderedSpecTemp as explained in Algorithm 3, and we denote its solution as ${\mathbf{L}}^{*}_{\mathrm{ord}}$. To estimate the unknown observation time $T$, we leverage the fact that $T$ is an integer and perform a line search to optimize

[TABLE]

where $\hat{\mathbf{L}}_{t}$ is obtained via InverseFilter (cf. Algorithm 1) by assuming that $T=t$. The ratio of traces in (25) ensures that the scales of ${\mathbf{L}}^{*}_{\mathrm{ord}}$ and $\hat{\mathbf{L}}_{t}$ are comparable. Notice that for this problem setting where the diffusion rates $\alpha_{t}\equiv\alpha$ are constant, InverseFilter can be simplified as

[TABLE]

where $\hat{\sigma}^{2}$ is obtained via (9). Intuitively, if ${\mathbf{L}}^{*}_{\mathrm{ord}}$ were to coincide with the true Laplacian ${\mathbf{L}}$, then we would expect $\hat{T}$ in (25) to return the true observation time $T$. In the absence of such ground truth, OrderedSpecTemp provides an estimate of the Laplacian that relies on the available observation. However, Algorithm 3 does not leverage the fact that all the observations come from a single consensus dynamics with a constant diffusion rate (cf. Problem 2). Hence, once we have estimated $T$, we take advantage of this fact by implementing NearestCGL (cf. Algorithm 2) for $T=\hat{T}$, to obtain our estimated Laplacian ${\mathbf{L}}^{*}$. This hybrid algorithm is summarized in Algorithm 4.

Remark 5

It should be noted that the assumption $\alpha_{t}\equiv\alpha$ is made only in the context of Problem 2. For such dynamics whose diffusion rates do not change over time, we can assume that $\alpha_{t}$ is constant and adopt Algorithm 4. When we do not have such prior information, we can leverage Algorithm 3 which imposes less assumptions on the diffusion parameters. Intuitively, in Problems 1-3, we assume different levels of knowledge about $\{\alpha^{(k)}_{t}\}$ to find a trade-off between model flexibility and algorithm performance. In fact, we can use Algorithm 3 to solve all of the three problems. However, if we do have additional prior information about the model (like in Problems 1 and 2), we can incorporate it into our algorithm to further improve the estimation accuracy. This point is validated in our experiments; see Fig. 2.

VI Numerical experiments

We illustrate the performance of the methods described and compare them with state-of-the-art solutions through a series of numerical experiments. The experiments based on synthetic data are divided by the problem type that is being studied (Sections VI-A-VI-C) whereas in Section VI-D we implement our methods to real-world data. In addition, we discuss the scalability of the proposed algorithms and infer a larger real-world network in Section VI-E.

VI-A Experiments on Problem 1

We compare the performance of InverseFilter (Algorithm 1), NearestCGL (Algorithm 2), and the approach proposed in [23] for estimating the network topology from a single consensus dynamics. The method in [23] assumes that the graph filter is a one-to-one and non-negative function of the CGL $\mathbf{L}$ with one unknown parameter, and proposes an approach to jointly estimate $\mathbf{L}$ and the unknown parameter from filtered signals. We refer to this approach as StructGLasso in this paper333We thank the authors of [23] for kindly sharing their code with us., and it can be directly adopted for our problem setting where the parameters in the graph filter are given. The StructGLasso method solves the following problem to estimate $\mathbf{L}$,

[TABLE]

where $\hat{\mathbf{L}}$ is the estimate obtained via InverseFilter, and $|\mathbf{L}|$ denotes the pseudo-determinant of $\mathbf{L}$. Notice that the idea behind (26) is to incorporate structural information into graphical LASSO [8, 43].

We generate synthetic datasets based on the signal model in Problem 1. We consider random Erdős-Rényi (ER) graphs [44] of size $N=36$ and edge-formation probability $p=0.1$. The edge weights are randomly selected from a uniform distribution in the interval $(0.1,3)$. We consider white Gaussian input, i.e. $\bm{\xi}\sim\mathrm{N}(\mathbf{0},\sigma^{2}\mathbf{I})$, and we set $\sigma=1$. Therefore, the dataset entries $\{\mathbf{y}_{k}\}_{k=1}^{M}$ are randomly drawn from the distribution $\mathrm{N}(\mathbf{0},{h(\mathbf{L})^{2}})$ where $h(\mathbf{L})=\prod_{t=1}^{T}(\mathbf{I}-\alpha_{t}\mathbf{L})$. We set $T=3$ and $\{\alpha_{1},\alpha_{2},\alpha_{3}\}=\{0.7,0.8,0.9\}/\lambda_{\max}(\mathbf{L})$. To measure the performance, we compute the recovery error as $\|{\mathbf{L}}^{*}-\mathbf{L}\|_{\mathrm{F}}/\|\mathbf{L}\|_{\mathrm{F}}$ where $\mathbf{L}$ is the ground truth and ${\mathbf{L}}^{*}$ is the estimate. For NearestCGL and StructGLasso, we consider two cases where (i) there is no regularization ($\beta=0$), and (ii) the regularization parameter $\beta$ is carefully selected. For StructGLasso, $\beta$ is chosen from the following set

[TABLE]

as suggested in [43, 23] where $s_{\max}=\max_{i\neq j}|[\hat{\mathbf{L}}^{{\dagger}}]_{ij}|$ is the maximum off-diagonal entry of $\hat{\mathbf{L}}^{{\dagger}}$ in absolute value. For NearestCGL, we consider two choices for $d(\mathbf{L},\hat{\mathbf{L}})$ in (15), namely $\|\mathbf{L}-\hat{\mathbf{L}}\|_{\mathrm{F}}$ and $\|\mathbf{L}-\hat{\mathbf{L}}\|_{2}$. For NearestCGL with $\|\cdot\|_{\mathrm{F}}$, $\beta$ is selected from $\{0,0.055\!:\!0.0025\!:\!0.085\}$ for $M/N\leq 3$ and from $\{0\!:\!0.01\!:\!0.08\}$ otherwise. For NearestCGL with $\|\cdot\|_{2}$, $\beta$ is selected from $\{0,0.008\!:\!0.001\!:\!0.02\}$ for $M/N\leq 3$ and from $\{0\!:\!0.002\!:\!0.016\}$ otherwise. We perform Monte-Carlo simulations and compute the average recovery error over $20$ realizations. The methods are all implemented using CVX [45] and the results are shown in Fig. 1.

It can be observed that, for small sample size, the regularization improves the estimation performance significantly. For NearestCGL, the matrix distance based on the Frobenius norm outperforms the one based on the spectral norm. Moreover, the proposed NearestCGL with Frobenius norm and a fine tuned regularization parameter outperforms all of the other methods for every sample size considered.

Remark 6** **(Related work [23])

Both our proposed NearestCGL and the algorithm in [23] can be applied to the class of graph filters whose corresponding functions are nonnegative and one-to-one if the type and parameters of the graph filter are known. The first step of these two algorithms are similar (i.e., InverseFilter) while their second steps differ significantly. It can be observed from (15) and (26) that our second step directly relies on $\hat{\mathbf{L}}$ while (26) is based on the pseudo-inverse $\hat{\mathbf{L}}^{{\dagger}}$. Compared to [23], we use different functions $d(\cdot,\cdot)$ to measure the difference between two matrices. We adopt matrix norms, whereas [23] uses a log-determinant Bregman divergence [39]. The specific graph filter type will influence the estimation accuracies of $\mathbf{L}$ and $\mathbf{L}^{{\dagger}}$ obtained in the first step, which causes the performance differences between our algorithm and the one in [23]. Roughly speaking, [23] shows the advantage of its proposed algorithm for a few types of graph filters which explicitly involve the pseudo-inverse of $\mathbf{L}$. This is because, for such graph filters, the estimation accuracy of $\mathbf{L}^{{\dagger}}$ obtained by InverseFilter is generally higher than that of $\mathbf{L}$. One the other hand, our model can be extended beyond consensus to a wider range of graph filters that are explicit functions of $\mathbf{L}$ including heat diffusion $h(\mathbf{L})=\exp(-\tau\mathbf{L})$. In this respect, both approaches seem to complement each other since they work better for different classes of filters.

For conciseness we only present the results for white Gaussian input and ER graphs, while the observed results are preserved when considering white uniform input – i.e., each entry in the input $\bm{\xi}$ is randomly generated from a uniform distribution in the interval $[-\sqrt{3}\sigma,\sqrt{3}\sigma]$ – and other graph models including grid graphs, stochastic block model (SBM), small-world, and Barabási-Albert model [44].

To study the influence in the performance of a deviation from the independent, white Gaussian input assumption (Assumption 2), we consider two types of inputs: (i) we generate the input as a time series of the form $\bm{\xi}_{t}=a\bm{\xi}_{t-1}+(1-a)\mathbf{w}_{t}$ where $\bm{\xi}_{1}$ and $\{\mathbf{w}_{t}\}$ are randomly drawn from $\mathrm{N}(\mathbf{0},\mathbf{I})$; and (ii) we generate the input following the distribution $\mathrm{N}(\mathbf{0},\mathbf{C}_{\boldsymbol{\xi}})$ where $\mathbf{C}_{\boldsymbol{\xi}}$ is drawn from the Wishart distribution $\mathrm{W}_{N}(\mathbf{I},d)$ and normalized by $d$. The other simulation parameters are the same as those used in generating Fig. 1. The results are shown in Figs. 1 and 1, respectively.

It can be observed that, for case (i), the performance of NearestCGL decreases as the parameter $a$ increases. This is as expected since larger values of $a$ introduce more correlation between successive inputs, thus deviating from the independent input assumption ($a=0$). However, regardless of the value of $a$, the recovery error decreases as the sample size $M$ increases. Intuitively, this can be attributed to the fact that $\bm{\xi}_{t_{1}}$ and $\bm{\xi}_{t_{2}}$ are close to independent if they are separated by a long time period. Hence, NearestCGL can be seen to output a consistent estimator under this setting.

For case (ii), the recovery error of NearestCGL decreases as the parameter $d$ increases. This is also natural since the resulting $\mathbf{C}_{\boldsymbol{\xi}}$ is closer to the identity matrix for larger $d$, and the white case $\mathbf{C}_{\boldsymbol{\xi}}=\mathbf{I}$ can be considered as the extreme setting when $d=\infty$. For finite $d$, the recovery error does not converge to zero for large $M$. This follows from (5), where a colored input implies that ${\mathbf{C}}_{y}$ does not share the same set of eigenvectors as ${\mathbf{L}}$. Hence, NearestCGL does not output a consistent estimator in this setting.

VI-B Experiments on Problem 2

In evaluating the performance of the method proposed in Section V-B, we first analyze the accuracy of the estimator $\hat{T}$ of the observation time [cf. (25)]. We consider three parameter settings: (i) $T=5$ and ER graphs of size $N=36$ and edge-formation probability $p=0.2$; (ii) $T=4$ and SBM graphs with $N=36$ nodes and two blocks where the vertex attachment probabilities across blocks and within blocks are $p_{1}=0.1$ and $p_{2}=0.3$, respectively; and (iii) $T=3$ and small-world graphs generated via Watts-Strogatz model with $N=36$ nodes, mean node degree $4$ and rewiring probability $0.2$. For all settings, we adopt unweighted graphs and the diffusion rate is chosen as $\alpha=0.8/\lambda_{\max}(\mathbf{L})$. We set $T_{\mathrm{max}}=10$ in (25). In solving (25), we must first obtain ${\mathbf{L}}^{*}_{\mathrm{ord}}$ via OrderedSpecTemp. Consequently, in (23) we set $\eta=1$, $d(\mathbf{L},\mathbf{K})=\|\mathbf{L}-\mathbf{K}\|_{2}$, and $\epsilon$ equal to the smallest possible value in the set $\{0.002\!:\!0.002\!:\!0.03\}\cup\{0.03+0.005\,r\,|\,r=1,2,\cdots\}$ that guarantees feasibility of (23). One iteration of a reweighted $\ell_{1}$ minimization scheme [41] is adopted. For each setting, we implement $100$ realizations and compute the ratio of successful recovery of $T$. The results are shown in Fig. 2. We can see a sharp increase on the probability of successful identification with increasing the sample size, eventually converging to $1$.

Furthermore, for setting (i), we compare the recovery performance of three methods, namely, OrderedSpecTemp (Algorithm 3), the proposed hybrid algorithm (Algorithm 4), and NearestCGL (Algorithm 2) with the true $T$. Notice that this latter method is not applicable in practice (since $T$ is unknown) but rather serves as a benchmark for ideal performance. In NearestCGL as well as the hybrid method, we simply set $\beta=0$. The results (averaged across $30$ realizations) are shown in Fig. 2. We can see that the performance of the proposed hybrid approach is better than that of OrderedSpecTemp, effectively corroborating that ${\mathbf{L}}^{*}$ is a better estimate than ${\mathbf{L}}^{*}_{\mathrm{ord}}$ in Algorithm 4. Moreover, the hybrid approach is very close to the benchmark of NearestCGL with perfect knowledge, especially in the large sample size regime.

VI-C Experiments on Problem 3

We first test the proposed OrderedSpecTemp method (Algorithm 3) when the eigenbasis $\mathbf{V}$ is perfectly known, i.e., when ${\mathbf{U}}={\mathbf{V}}$ in (23c). This situation corresponds to the infinite sample size regime (cf. Theorem 3). A high recovery rate under this setting is important since it is a necessary condition for acceptable recovery in the finite sample size regime. We compare OrderedSpecTemp with two other methods: (i) SpecTemp+LEigVec: in this method, we assume that the full order information is unavailable and only the index of the leading eigenvector is known, so we construct the problem similar to (23) while in (23d) we only keep the constraint $\gamma_{N}=1$; (ii) The method proposed in [1], where the constraint (23d) is replaced by $\gamma_{i}\leq\gamma_{i+\eta}-\epsilon_{2}\text{ for }i=1,\cdots,N-\eta$, and $\epsilon_{2}>0$ can be chosen freely since it only affects the scale of the recovered Laplacian. One advantage of (23d) over the method in [1] is that (23d) can be directly adopted in the case when the CGL $\mathbf{L}$ has repeated eigenvalues while [1] requires the assumption that $\mathbf{L}$ has distinct eigenvalues.

We consider random unweighted ER graphs of varying sizes $N\in\{10,20,\cdots,50\}$ and edge-formation probabilities $p\in\{0.1,0.2,\cdots,0.5\}$. Since the method in [1] is proposed under the assumption that the CGL has distinct eigenvalues, we generate graphs satisfying $\min_{i\neq j}|\lambda_{i}-\lambda_{j}|\geq 10^{-4}$. We set parameters in (23) as $\mathbf{U}=\mathbf{V}$, $\epsilon=0$ and $\eta=1$. For the method in [1], we set $\epsilon_{2}=1$. For all three methods, three iterations of a reweighted $\ell_{1}$ minimization scheme [41] are adopted. For each graph generated, we consider the recovery to be successful if $\|{\mathbf{L}}^{*}-\mathbf{L}\|_{\mathrm{F}}/\|\mathbf{L}\|_{\mathrm{F}}<0.02$. Note that the estimate $\mathbf{L}^{\ast}$ here is scaled to have the same trace as the true $\mathbf{L}$ in order to compute the recovery error. The results (averaged across $50$ realizations) are shown in Fig. 3. We can see that the proposed OrderedSpecTemp method outperforms the other two methods, and its recovery rate is equal to or near $1$ for all graph settings. Notice that the SpecTemp+LEigVec method works better for sparser graphs. This indicates that the eigenvector order information (ignored by SpecTemp+LEigVec) becomes more important for denser graphs.

We now test OrderedSpecTemp in the finite sample size regime, and compare it with SpecTemp+LEigVec to further study the value of incorporating the eigenvector ordering information. We consider unweighted ER graphs of size $N=36$ and varying edge-formation probabilities $p=\{0.1,0.3,0.5\}$. For each graph, our goal is to recover the CGL $\mathbf{L}$ from the observation of $M$ synthetic consensus dynamics, where we vary $M$ from $N$ to $10^{3}N$. For each dynamics, we draw the input from a standard multivariate Gaussian distribution, selecting $T_{k}$ uniformly at random in $\{3,4,5\}$ and each diffusion rate $\alpha^{(k)}_{t}$ uniformly at random in $(0,1/\lambda_{\max}(\mathbf{L}))$. We set $\eta=1$ for OrderedSpecTemp, use the spectral norm as the matrix distance, and set $\epsilon$ as the smallest possible value in the set $\{0.005\,r\,|\,r=1,2,\cdots\}$ that guarantees feasibility of the optimization problem for both two methods [cf. (23)]. Apart from the recovery error, we also consider the F-score metric, which is commonly used to evaluate binary classification performance and computed as

[TABLE]

where $\mathrm{tp}$, $\mathrm{fp}$ and $\mathrm{fn}$ respectively represent true-positive, false-positive and false-negative detection of edges in the estimate ${\mathbf{L}}^{*}$ with respect to edges in the true $\mathbf{L}$. F-score takes values between [math] and $1$, and $1$ represents perfect classification. The results (averaged across $20$ realizations) are shown in Fig. 4.

As displayed in Fig. 4, for both methods and all graph density settings, we observe a monotonous decrease of the recovery error and a monotonous increase of the F-score with increasing number of samples $M$. This is not surprising since we know that for larger sample size, the eigenbasis of the sample covariance becomes closer to the eigenbasis of the CGL and thereby facilitates recovery. We can also see that, OrderedSpecTemp outperforms SpecTemp+LEigVec and the performance gap increases when the graph becomes dense. This is consistent with the results in Fig. 3, i.e. SpecTemp+LEigVec works well for very sparse graphs, while when the graph density increases, the full order information becomes more valuable in aiding recovery.

VI-D Topology inference from real-world data

We present two different real-world case studies where we recover connections between the states of the U.S. in both cases but based on very different information sources, namely, temperature measurements and congressional roll-call votes.

Temperature network. We apply the proposed OrderedSpecTemp method on a real-world dataset consisting of $M=5844$ average daily temperature measurements collected from $N=45$ states in the U.S. over the years $2000$-$2015$ [46]. This same dataset was analyzed in [23]. The temperature signals are spatially smooth across different states, and the Rocky Mountains region (which is mainly located in the states of Montana, Wyoming, Colorado, and Idaho) has lower average temperature values than its geographical neighborhood (cf. Fig. 7 in [23]).

In solving (23), we set $\eta=1$, $d(\mathbf{L},\mathbf{K})=\|\mathbf{L}-\mathbf{K}\|_{2}$, and $\epsilon$ equal to the smallest possible value (found via five iterations of binary search between [math] and $1$) that guarantees feasibility of (23). The result shown in Fig. 5 is obtained after three iterations of a reweighted $\ell_{1}$ minimization scheme  [41].

As shown in Fig. 5, the edges connecting neighboring states generally have larger weights since temperature values are similar between regions located near to each other. In this sense, temperature values are accurately described by a consensus dynamics where discrepancies between neighbors tend to be reduced. However, there are other factors – besides the geographical distance – that can also influence the similarity of temperature values, such as landform and altitude. It can be observed that the edges between the Rocky Mountain states and their neighbors to the east have relatively small weights. In this way, the recovered network captures the natural barrier for temperature similarity imposed by the mountain formation.

As mentioned, in [23] the authors apply their proposed method – referred to as StructGLasso in this paper – to this same dataset by assuming different graph filter functions and parameters. The network recovered here is comparable to the one inferred in [23] for an exponential decay filter, which was deemed as better revealing the structure of the signal (cf. Fig. 8 (d)-(f) in [23]). This implies that, in the absence of prior knowledge about the specific filter $h(\mathbf{L})$ that is driving the underlying network process, adopting a more versatile model where only a decaying frequency response is assumed [cf. (23d)] can be beneficial. On the contrary, if the filter type is known in advance, then it is reasonable to incorporate it in the recovery method, as advocated in [23].

Senate network. We now apply OrderedSpecTemp on a real-world dataset from congressional roll-call votes in the U.S. [47]. We use the roll-call data of $100$ senators ($2$ per state) in the $114$th congress ($2015$-$2017$), consisting of $M=502$ roll-calls. We quantify yea, nay and other cases (e.g. abstention) as $1$, $-1$ and [math], respectively, to represent each senator’s opinion. Since each state has two senators, we use the sum of their opinions as the graph signal value of their state, resulting in a graph with $N=50$ nodes. We choose the parameters in (23) following the same process explained for the temperature network.

We divide the states into three categories: (i) labeled as D if both senators in this state are from the Democratic Party, (ii) labeled as R if both senators in this state are from the Republican Party, (iii) labeled as M if the senators are from different (mixed) parties. The recovered network contains $391$ edges, and we count the number of edges inside each category and between different categories in the top-$K$ edges sorted by weight for different $K$, as shown in Fig. 5. It can be observed that there are more edges inside the categories D and R, while less edges between them. In Fig. 5, we plot the recovered network with the top-$150$ edges in spring layout using NetworkX [48]. It can be observed that two tight clusters of states emerge (D and R) with a looser cluster of mixed states M connecting these two.

VI-E Scalability of the proposed algorithms

To illustrate the scalability of the proposed algorithms, we infer a large-scale real-world network under the setting of Problem 1 [cf. (16)]. In particular, we consider the U.S. political blogs network444The data can be found in http://www-personal.umich.edu/~mejn/netdata/. which is a directed network of hyperlinks between weblogs on US politics [49]. We preprocessed the graph by: (i) converting the directed graph to an undirected one, (ii) removing self-loops and multi-edges, and (iii) keeping the largest connected component. After these changes, we obtain an undirected and unweighted graph with $1,\!222$ nodes and $16,\!714$ edges. We adopt a graph filter of the form $h(\mathbf{L})=(\mathbf{I}-\alpha\mathbf{L})^{T}$, where we set $\alpha=0.9/\lambda_{\max}(\mathbf{L})$ and $T=15$. We consider different values for the ratio of sample size and graph size, i.e., $M/N=1,2,3,4,5$. The corresponding values for the regularization parameter are selected as $\beta=40,30,25,20,5$, respectively. The results, averaged across 20 realizations, are shown in Fig. 6. We observe that the proposed NearestCGL outperforms the InverseFilter algorithm both in terms of recovery error and F-score. As an indication of scalability, we remark that solving the problem (16) once for this specific graph requires around $17$ seconds of computation in a standard laptop computer.

VII Conclusions

We proposed a set of algorithms for the identification of a network based on observing snapshots of consensus processes considering different levels of parameter uncertainty. To achieve this, we constructed convex optimization problems that output a sparse, valid graph Laplacian which is provably consistent with the spectral information obtained from the observed data. Finally, we showcased the effectiveness of the proposed methods in synthetic and real-world scenarios.

Potential future research avenues include: (i) investigation of different types of restrictions on the available data, such as observation models with snapshot data sampled from partial nodes of the graph; (ii) study of the trade-off between specific network topologies and the required sample size to achieve a desired level of estimation accuracy; (iii) consideration of a richer class of dynamical models, including non-deterministic processes such as switched systems; (iv) joint estimation of several related networks from the concurrent observation of consensus dynamics [50, 51]; and (v) topology inference of directed graphs.

Appendix A Random variables and properties

In this appendix, we review concepts and properties of random variables that we used throughout the paper; for a thorough review, see [52].

Definition 1** **(Sub-Gaussian random variable)

A random variable $x$ with mean $\mu=\mathbb{E}[x]$ is sub-Gaussian if there is a positive parameter $\sigma$ such that

[TABLE]

Lemma 2** **(Sub-Gaussian tail bound [52])

Suppose that $x$ is sub-Gaussian with mean $\mu$ and parameter $\sigma$, then

[TABLE]

Definition 2** **(Sub-exponential random variable)

A random variable $x$ with mean $\mu=\mathbb{E}[x]$ is sub-exponential if there are positive parameters $(\nu,b)$ such that

[TABLE]

Lemma 3** **(Sub-exponential tail bound [52])

Suppose that $x$ is sub-exponential with mean $\mu$ and parameters $(\nu,b)$, then

[TABLE]

Lemma 4** **(​​[1])

Let $x\sim\mathrm{N}(0,\sigma_{1}^{2})$ and $y\sim\mathrm{N}(0,\sigma_{2}^{2})$ be independent random variables. Then their product $z=xy$ is sub-exponential with parameters $(\sqrt{2}\sigma_{1}\sigma_{2},\sqrt{2}\sigma_{1}\sigma_{2})$.

Lemma 5** **(​​[1])

Let $x\sim\mathrm{N}(0,\sigma^{2})$, then its square $z=x^{2}$ is sub-exponential with mean $\mathbb{E}(z)=\sigma^{2}$ and parameters $(2\sigma^{2},4\sigma^{2})$.

Lemma 6** **(​​[52])

Let $\{x_{k}\}_{k=1}^{M}$ be independent sub-exponential random variables with mean $\mu_{k}$ and parameters $(\nu_{k},b_{k})$ respectively. Then the variable $\sum_{k=1}^{M}(x_{k}-\mu_{k})$ is sub-exponential with parameters $(\sqrt{\sum_{k=1}^{M}\nu_{k}^{2}},\,\max\limits_{k}b_{k})$.