State-Space Network Topology Identification from Partial Observations

TL;DR

This paper presents a novel method for recovering network topology from partial observations using extended subspace techniques, bridging control theory and graph signal processing with theoretical guarantees and an convergent algorithm.

Contribution

It introduces a unified approach combining system identification and graph signal processing to recover network structure with provable guarantees from partial data.

Findings

Algorithm successfully recovers network topology from partial observations.

Theoretical guarantees ensure accurate topology identification under certain conditions.

Numerical results validate the effectiveness and applicability of the proposed method.

Abstract

In this work, we explore the state-space formulation of a network process to recover, from partial observations, the underlying network topology that drives its dynamics. To do so, we employ subspace techniques borrowed from system identification literature and extend them to the network topology identification problem. This approach provides a unified view of the traditional network control theory and signal processing on graphs. In addition, it provides theoretical guarantees for the recovery of the topological structure of a deterministic continuous-time linear dynamical system from input-output observations even though the input and state interaction networks might be different. The derived mathematical analysis is accompanied by an algorithm for identifying, from data, a network topology consistent with the dynamics of the system and conforms to the prior information about the underlying structure. The proposed algorithm relies on alternating projections and is provably convergent. Numerical results corroborate the theoretical findings and the applicability of the proposed algorithm.

Equations

$${\partial }_t\mathbf{x}(t)$$

$$\mathbf{y}(t)$$

$$\mathcal{P}:=\{h(\cdot ),c(\cdot )\},$$

$${\partial }_t\mathbf{x}(t)$$

$$\mathbf{y}(t)$$

$$f_{\ast }({\mathbf{S}}_{\ast }):=\frac{1}{2\pi i}{\int }_{{\Gamma }_{f_{\ast }}}{f}_{\mathrm{s},\ast }(z)R(z,{\mathbf{S}}_{\ast })dz,$$

$$R(z,{\mathbf{S}}_{\ast }):=({\mathbf{S}}_{\ast }-z\mathbf{I}{)}^{-1}.$$

$$\mathbf{x}(k+1)$$

$$\mathbf{y}(k)$$

$$\mathbf{A}({\mathbf{S}}_i)$$

$$\mathbf{B}({\mathbf{S}}_j)$$

$$\mathbf{x}(k+1)$$

$$\mathbf{y}(k)$$

$${\mathbf{y}}_{k,\alpha }\triangleq [\mathbf{y}(k{)}^T,\dots ,\mathbf{y}(k+\alpha -1{)}^T{]}^T,$$

$${\mathbf{y}}_{k,\alpha }={\mathcal{O}}_{\alpha }\mathbf{x}(k)+{\mathcal{T}}_{\alpha }{\mathbf{u}}_{k,\alpha }+{\mathbf{n}}_{k,\alpha },$$

$${\mathcal{O}}_{\alpha }\triangleq \left[\begin{array}{c}\mathbf{C}\\ \mathbf{C}\mathbf{A}\\ ⋮\\ \mathbf{C}{\mathbf{A}}^{\alpha -1}\end{array}\right],$$

$${\mathcal{T}}_{\alpha }\triangleq \left[\begin{array}{cccccc}\mathbf{D}& 0& 0& \cdots & \cdots & 0\\ \mathbf{C}\mathbf{B}& \mathbf{D}& \ddots & & \cdots & 0\\ ⋮& \ddots & \ddots & & & ⋮\\ ⋮& & & & & 0\\ \mathbf{C}{\mathbf{A}}^{\alpha -2}\mathbf{B}& \mathbf{C}{\mathbf{A}}^{\alpha -3}\mathbf{B}& \cdots & \cdots & \mathbf{C}\mathbf{B}& \mathbf{D}\end{array}\right],$$

$${\mathbf{u}}_{k,\alpha }\triangleq [\mathbf{u}(k{)}^T,\dots ,\mathbf{u}(k+\alpha -1{)}^T{]}^T,$$

$$\mathbf{Y}=[{\mathbf{y}}_{1,\alpha },\dots ,{\mathbf{y}}_{T,\alpha }],$$

$$\mathbf{Y}={\mathcal{O}}_{\alpha }\mathbf{X}+{\mathcal{T}}_{\alpha }\mathbf{U}+\mathbf{N},$$

$$\mathbf{X}=[\mathbf{x}(1),\dots ,\mathbf{x}(T)],$$

$$\mathbf{U}=[{\mathbf{u}}_{1,\alpha },\dots ,{\mathbf{u}}_{T,\alpha }]\text{and}\mathbf{N}=[{\mathbf{n}}_{1,\alpha },\dots ,{\mathbf{n}}_{T,\alpha }].$$

$$\mathbf{Y}={\mathcal{O}}_{\alpha }\mathbf{X}+{\mathcal{T}}_{\alpha }\mathbf{U}.$$

$${\mathbf{\Pi }}_{{\mathbf{U}}^T}^{\perp }\triangleq \mathbf{I}-{\mathbf{U}}^T(\mathbf{U}{\mathbf{U}}^T{)}^{-1}\mathbf{U}.$$

$$\mathbf{Y}{\mathbf{\Pi }}_{{\mathbf{U}}^T}^{\perp }={\mathcal{O}}_{\alpha }\mathbf{X}{\mathbf{\Pi }}_{{\mathbf{U}}^T}^{\perp }.$$

$$\mathrm{span}(\mathbf{Y}{\mathbf{\Pi }}_{{\mathbf{U}}^T}^{\perp })=\mathrm{span}({\mathcal{O}}_{\alpha }).$$

$$\mathbf{Y}{\mathbf{\Pi }}_{{\mathbf{U}}^T}^{\perp }=(\mathbf{Y}{\mathbf{\Pi }}_{{\mathbf{U}}^T}^{\perp }{)}_N={\mathbf{W}}_{\alpha ,N}{\mathbf{\Sigma }}_N{\mathbf{V}}_{\alpha ,N}^T,$$

$${\mathcal{O}}_{\alpha }={\mathbf{W}}_{\alpha ,N}\mathbf{T},$$

$$\hat{\mathbf{T}}=(\mathbf{J}{\mathbf{W}}_{\alpha ,N}{)}^{†}\mathbf{C},$$

$${\mathbf{J}}_{\mathrm{u}}{\mathcal{O}}_{\alpha }\mathbf{A}={\mathbf{J}}_{\mathrm{l}}{\mathcal{O}}_{\alpha }\in {\mathbb{R}e}^{(\alpha -1)L\times N},$$

$$\hat{\mathbf{A}}=({\mathbf{J}}_{\mathrm{u}}{\mathbf{W}}_{\alpha ,N}\hat{\mathbf{T}}{)}^{†}{\mathbf{J}}_{\mathrm{l}}{\mathbf{W}}_{\alpha ,N}\hat{\mathbf{T}}.$$

$$\mathbf{y}(k)-(\mathbf{u}(k{)}^T\otimes {\mathbf{I}}_L)$$

$$\mathbf{C}{\mathbf{A}}^k\mathbf{x}(0)+$$

Full text

I Introduction

The topology of networks is fundamental to model the interactions between entities and to improve our understanding about the processes defined over them. We can find examples of such processes in transportation networks [1], brain activity [2], and epidemic dynamics or gene regulatory networks [3], to name a few. The coupling between the process and the network topology has led to the extension of signal processing (SP) techniques to tools that take into account the network structure to define signal estimators [4, 5, 6], filters [7, 8, 9], and optimal detectors [10, 11, 12].

While in several scenarios, the network structure is available, in many others, it is unknown and needs to be estimated. This not only for enhancing data processing tasks but also for data interpretability, i.e., the network topology provides an abstraction for the underlying data dependencies. Therefore, retrieving the network structure or the dependencies of the involved members (variables) has become a research topic of large interest [13, 14, 15, 16, 17, 18, 19, 20, 21].

Despite that many works focus on the problem of topology identification [22, 23] or Gaussian graphical modeling (GGM) [24, 25], most of these approaches only leverage a model based on so-called graph filters [26] or enforce a particular structure by a penalized likelihood approach with sparsity constraints [27, 28]. Among the works that consider an alternative interaction model e.g., [21, 17, 29, 30], only few of them considers the network data as states of an underlying process. However, none of the above works study the case where the input, i.e., excitation or probing signal of the process and the process itself evolve according to different topologies.

In many instances, physical systems can be defined through a state-space formulation with known dynamics. An example is the diffusion of molecules in tissue. This process is used to analyze brain functions by mapping the interaction of the molecules with obstacles [31]. Besides the mapping of the brain, the area of neural dynamics presents the problem of network design in transport theory [32]. In such applications, the topology that provides a stable desired response needs to be found, following a differential equation. Finally, we recall the problem of finding the connections between reactants in chemical reaction networks [33]. Here, the evolution of the number of the molecules in a solution is governed by the interaction of the reactants present in it. Hence, to understand the underlying chemical process, the relation between reactants is required. Considering these examples, it is clear that a more general approach, parting from first-principles, to find the underlying connections is required. So, in this work, we focus on the problem of retrieving the network structure of a process modeled through a deterministic continuous-time dynamical system whose system matrices depend on the underlying topology.

We further remark that all above works estimate the network topology under the assumption that the observations are available on all entities. We here depart from this assumption and consider the problem of network topology inference from partial observations. In this way, the focus is not only on retrieving a topology that describes the dynamics of the process under full network access, but also for cases where observations from a subset of the entities are only available. To address these tasks, we devise a framework that allows estimating the network topology from partial observations up to ambiguities defined through an equivalence class of restricted cospectral graphs [34].

I-A Overview and main contributions

Network topology inference from partial measurements is generally an ill-posed problem and has not been addressed by existing methods. In this work, we introduce a first-principles approach based on state-space models for model-driven topology estimation. Our contributions that broaden the state-of-the-art are threefold.

• –

We propose a first-order differential graph formulation for systems whose dynamics are described by deterministic continuous-time linear models. Under this model, it is possible to obtain a first-principles based network topology identification framework that leverages the geometric structure in its state-space description. Further, we provide conditions under which we can retrieve the network topology from sampled observations.

• –

We extend subspace techniques from the system identification literature to network topology inference. These methods leverage the embedded geometrical structure of the state-space model that describes the dynamics. We provide ways to enforce the underlying dynamical structure on the topology identification process thereby guaranteeing system consistency, i.e., we enforce the structure present in the input-output data relations such that the estimated topology achieves the same dynamics as the original process.

• –

We introduce the problem of network topology identification from partial observations. We mathematically analyze this problem and show it is ill-posed. We further describe the ambiguities present when recovering the network topology for measurements that do not uniquely identify the underlying structure. We develop an algorithm to estimate the network structure in the current setting. This algorithm relies on the alternating projections (AP) approach [35] and is provably globally convergent. We prove that under mild conditions the AP method converges locally with a linear rate to a feasible solution. Finally, we extend the topology inference problem from partial observations to scenarios where incomplete or inaccurate process dynamics are present. For these cases, we provide a mathematical analysis and conditions that guarantee the convergence of the AP method in estimating a feasible network topology.

I-B Outline and notation

This paper is organized as follows. Section II formulates the problem of network topology inference for continuous-time dynamical systems from sampled observations. Section III introduces a first-order differential graph model and its state-space description. The system identification framework for the proposed graph-based model is introduced in Section IV. Section V analyzes the ambiguity of network topology identification from partial observations and provides an AP method to find a feasible network structure. Section VI discusses system consistency constraints that can be enforced into the AP method to match the network dynamics. Section VII corroborates the derived theory with numerical results. Finally, section VIII concludes the paper.

We adopt the following notation. Scalars, vectors, matrices and sets are denoted by lowercase letters $(x)$, lowercase boldface letters $(\bm{x})$, uppercase boldface letters $(\bm{X})$, and calligraphic letters $(\mathcal{X})$, respectively. $[\bm{X}]_{i,j}$ denotes the $(i,j)$th entry of the matrix $\bm{X}$ whereas $[\bm{x}]_{i}$ represents the $i$th entry of the vector $\bm{x}$. $\bm{X}^{T}$ and $\bm{X}^{-1}$ are the transpose and the inverse of $\bm{X}$, respectively. The Moore-Penrose pseudoinverse of $\bm{X}$ is denoted by $\bm{X}^{\dagger}$. ${\rm vec}(\cdot)$ is the vectorization operation. ${\rm bdiag}(\bm{X},\bm{Y})$ denotes a block diagonal matrix whose blocks are given by the matrices $\bm{X}$ and $\bm{Y}$. $\bm{I}$ is the identity matrix of appropriate size. $\|\bm{X}\|_{\rm F}$ and $\|\bm{X}\|_{2}$ denote the Frobenius- and $\ell_{2}$-norm of $\bm{X}$, respectively. ${\rm span}(\cdot)$ and ${\rm rank}(\cdot)$ are the span and rank of a matrix, respectively. Finally, we use $[K]$ to denote the set $\{1,2,\ldots,K\}$ and $\mathcal{D}_{K}$ to denote the set of $K\times K$ diagonal matrices.

II Problem Statement

Consider a set of $N$ nodes ${\mathcal{V}}=\{v_{1},\ldots,v_{N}\}$ representing cooperative agents such as sensors, individuals, and biological structures. On top of these agents, a process ${\mathcal{P}}$ is defined that describes the evolution through time of the agent signals ${\bm{x}}(t)\in\mathbb{R}^{N}$. The signal ${\bm{x}}(t)$ is such that the $i$th component $x_{i}(t)$ represents the signal evolution of agent $v_{i}$. The agent interactions w.r.t. the evolution of the signal ${\bm{x}}(t)$ are captured by a graph $G_{x}=({\mathcal{V}},{\mathcal{E}}_{x})$, where ${\mathcal{E}}_{x}$ is the edge set of this graph. We consider that process ${\mathcal{P}}$ is represented by a first-order differential model

[TABLE]

where $\partial_{t}{\bm{x}}(t):=d{\bm{x}}(t)/dt$ and $h(\cdot)$ and $c(\cdot)$ are maps that describe respectively the dynamics of the signal ${\bm{x}}(t)$ and the observables ${\bm{y}}(t)$. In model (1), ${\bm{u}}(t)$ is the (known) system input and ${\mathcal{E}}_{u}$ is the edge set of another graph $G_{u}=({\mathcal{V}},{\mathcal{E}}_{u})$ that captures the interactions between the elements of ${\mathcal{V}}$ for ${\bm{u}}(t)$. Put simply, process (1) describes the evolution of the signal ${\bm{x}}(t)$ under the influence of the maps $h(\cdot)$ and $c(\cdot)$ and the network topologies $G_{x}=({\mathcal{V}},{\mathcal{E}}_{x})$ and $G_{u}=({\mathcal{V}},{\mathcal{E}}_{u})$. For future reference, we will represent both graphs as $G_{*}=\{{\mathcal{V}},{\mathcal{E}}_{*}\}$, where $``*"$ is a space holder for $x$ and $u$.

We can differently define the process ${\mathcal{P}}$ through the set

[TABLE]

which contains the interactions in the system. While process ${\mathcal{P}}$ describes a continuous-time process, we usually have access to sampled realizations of it, i.e., the observables ${\bm{y}}(t)$ are collected on a finite set of time instances $\mathcal{T}~{}:=~{}\{t_{1},t_{2},\ldots,t_{T}\}$. As a result, we might also want to consider discrete-time approximations of (2), where we either have a discrete-time realization of ${\mathcal{P}}$ and/or a finite number of observables $\mathcal{Y}:=\{{\bm{y}}(t)\}_{t\in\mathcal{T}}$.

With this in place, we ask the following question: how can we retrieve the network topologies ${\mathcal{G}}_{x}$ and ${\mathcal{G}}_{u}$ for the agent signal ${\bm{x}}(t)$ and input signal ${\bm{u}}(t)$ given the known process ${\mathcal{P}}$ [cf. (2)], input signal $\bm{u}(t)$, and observables in $\mathcal{Y}$?

In this work, we answer the above question by employing results from Hankel matrices [36] and linear algebra whose foundations lie in the system identification theory [37]. The approach we will consider is based on subspace techniques that are by definition cost function free. This differs from the commonly used techniques in network topology identification [22] where the learned topology heavily depends on the considered cost function (e.g., smoothness or sparsity). If this prior knowledge is incorrect, it might lead to structures not related to physical interactions. The adopted techniques provide theoretical insights in when and how the underlying network structures can be identified. Furthermore, they have two other benefits. First, they impose none parameterization on the dynamical model and, therefore, avoid solving nonlinear optimization problems as in prediction-error methods [38]. Second, they allow identifying $G_{*}=\{{\mathcal{V}},{\mathcal{E}}_{*}\}$ from sampled data, i.e., by having access only to a subset of the observables ${\bm{y}}(t)$.

We finally remark that although more general models can be considered, e.g., higher-order differential models, we restrict ourselves to first-order models to ease exposition and provide a thorough and stand-alone work. Nevertheless, first-order differential models are of broad interest as they include diffusion processes –see [39] and references therein.

III First Order Differential Graph Model

The structure of $G_{*}=\{{\mathcal{V}},{\mathcal{E}}_{*}\}$ is mathematically represented by a matrix $\bm{S}_{*}$ (sometimes referred to as a graph shift operator [40, 41]) that has as candidates the graph adjacency matrix, the graph Laplacian, or any other matrix that captures the connections (relations) of the elements in the network. We then consider that process ${\mathcal{P}}$ is described through the linear continuous-time dynamical system

[TABLE]

where ${\bm{C}}\in\mathbb{R}^{L\times N}$ and ${\bm{D}}~{}\in~{}\mathbb{R}^{L\times N}$ are system matrices related to the observables $\bm{y}(t)$. The matrix function $f_{*}:\mathbb{R}^{N\times N}\rightarrow\mathbb{R}^{N\times N}$ is defined via the Cauchy integral [42]

[TABLE]

where $f_{\rm s,*}(\cdot)$ is the scalar version of $f_{*}(\cdot)$ and is analytic on and over the contour $\Gamma_{f_{*}}$. Here, $R(z,\bm{S}_{*})$ is the resolvent of $\bm{S}_{*}$ given by

[TABLE]

From the definition of the dynamical system in (1), it can be seen that (3) is a proper representation of the family of first-order (linear) differential models, where the system matrices are matrix functions of the matrix representing the graph. Observe that although restricted to the class of linear models, (3) captures different settings of practical interest such as diffusion on networks [39], graph filtering operations [26, 9], random walks [43], and first-order autoregressive graph processes [44].

The corresponding discrete-time state-space system related to (3) is

[TABLE]

where $\tilde{f}_{*}(\cdot)$ is a matrix function (to be specified in the sequel) and ${\bm{x}}(k)\in{\mathbb{R}}^{N}$, ${\bm{u}}(k)\in{\mathbb{R}}^{N}$, and ${\bm{y}}(k)\in{\mathbb{R}}^{L}$ are the discrete counterparts of ${\bm{x}}(t)$, ${\bm{u}}(t)$, and ${\bm{y}}(t)$, respectively. The variables, $\bm{w}(k)$ and $\bm{v}(k)$ represent perturbations in the states and additive noise in the observables, respectively. By defining then the matrices $A(\bm{S}_{x}):=\tilde{f}_{x}(\bm{S}_{x})$ and $\bm{B}({\bm{S}_{u}}):=\tilde{f}_{u}(\bm{S}_{u})$, the connection between the continuous-time (3) and the discrete-time representation (6) is given by

[TABLE]

where $\tau$ is the sampling period and $e^{\bm{X}}=\sum_{k=0}^{\infty}\frac{1}{k!}\bm{X}^{k}$ is the matrix exponential function [42]. Using then (7) and (8), we can compactly write model (6) as

[TABLE]

where we dropped the dependency of the system matrices ${\bm{A}}$ and ${\bm{B}}$ from $\bm{S}_{x}$ and $\bm{S}_{u}$ to simplify the notation. Throughout the paper, we assume that the matrices ${\bm{C}}$ and ${\bm{D}}$ are known, i.e., we know how our observables are related to the states and inputs..

IV State-Space Identification

The family of subspace state-space system identification methods [45] relies on geometrical properties of the dynamical model (9). By collecting a batch of $\alpha$ different observables ${\bm{y}}(t)\in\mathbb{R}^{L}$ into the $\alpha L-$dimensional vector

[TABLE]

we get the relationship

[TABLE]

where

[TABLE]

is the extended observability matrix of the system [46] and

[TABLE]

is the matrix that relates the batch input vector

[TABLE]

with the batch observables $\bm{y}_{k,\alpha}$. The vector $\bm{n}_{k,\alpha}$ comprises the batch noise that depends on the system perturbation $\{\bm{w}(k)\}$ and on the observable noise $\{\bm{v}(k)\}$. The detailed structure of $\bm{n}_{k,\alpha}$ is unnecessary for our framework. The size of the batch $\alpha$ is a user-specified integer and must be larger than the number of states. Assuming the number of nodes is the number of states, this implies $\alpha>N$.

Given then expression (10) and the structures for $\mathcal{O}_{\alpha}$ in (11) and $\mathcal{T}_{\alpha}$ in (12), we proceed by estimating first ${\bm{A}}$ using the algebraic properties of (11) and subsequently ${\bm{B}}$ from the structure of (12) and a least squares problem.

Retrieving the state matrix ${\bm{A}}$. A basic requirement for estimating ${\bm{A}}$ is system observability [46]. Observability allows to infer the system state from the outputs for any initial state and sequence of input vectors. Put differently, we can estimate the entire system dynamics from input-output relations. System (9) is observable if the system matrices $\{{\bm{A}},{\bm{C}}\}$ satisfy ${\rm rank}(\mathcal{O}_{N})=N$.

Consider now a set of $Q\triangleq T+\alpha-1$ input-output pairs $\{\bm{y}(k),\bm{u}(k)\}_{k=1}^{Q}$. By stacking the discrete batch vectors $\bm{y}_{k,\alpha}$ for all observations into the matrix

[TABLE]

and using expression (10), we can generate the Hankel-structured data equation [47]

[TABLE]

where $\bm{X}$ is the matrix that contains the evolution of the states accross the columns, i.e.,

[TABLE]

and where the input $\bm{U}$ and noise $\bm{N}$ are block Hankel matrices defined as

[TABLE]

The structure in (13) is at the core of system identification methods [37] and this arrangement leads naturally to a subspace-based approach to find ${\bm{A}}$. To detail this, consider the noise-free case

[TABLE]

Since the control inputs are known and provided that ${\bm{U}}$ is full row-rank, we can project out $\bm{U}$ from (14) by right-multiplying $\bm{Y}$ with the projection matrix

[TABLE]

Since $\bm{U}\bm{\Pi}_{\bm{U}^{T}}^{\perp}=\bm{0}$, this operation leads to

[TABLE]

Under the assumption that the inputs $\bm{U}$ are sufficiently exciting [47] (the inputs excite all the modes of the system), the matrix $\bm{X}\bm{\Pi}_{\bm{U}^{T}}^{\perp}$ has full row-rank. This implies

[TABLE]

That is, the signal subspace of the projected observables $\bm{Y}\bm{\Pi}_{\bm{U}^{T}}^{\perp}$ coincides with that of the extended observability matrix. Therefore, ${\rm span}(\mathcal{O}_{\alpha})$ can be estimated from $\bm{Y}\bm{\Pi}_{\bm{U}^{T}}^{\perp}$ and, subsequently, the system matrix ${\bm{A}}$ by using the block structure of $\mathcal{O}_{\alpha}$ in (11). To detail this procedure, consider the economy-size singular value decomposition (SVD) of $\bm{Y}\bm{\Pi}_{\bm{U}^{T}}^{\perp}$

[TABLE]

where $(\bm{Y}\bm{\Pi}_{\bm{U}^{T}}^{\perp})_{N}$ is the $N$-rank approximation of $\bm{Y}\bm{\Pi}_{\bm{U}^{T}}^{\perp}$ and equality holds because of the full row-rank assumption of $\bm{X}\bm{\Pi}_{\bm{U}^{T}}^{\perp}$. Then, from condition (17), we have

[TABLE]

for some invertible similarity transform matrix $\bm{T}\in\mathbb{R}^{N\times N}$. Given then $\bm{C}$ known and full column-rank, we can estimate $\bm{T}$ from the structure of $\mathcal{O}_{\alpha}$ [cf. (11)] as

[TABLE]

where $(\cdot)^{\dagger}$ denotes the Moore-Penrose pseudoinverse and $\bm{J}\bm{W}_{\alpha,N}$ denotes the first $L$ rows of $\bm{W}_{\alpha,N}$. If $\bm{C}$ is not full column-rank, the transform matrix $\bm{T}$ is not unique. We will deal with the non-uniqueness of $\bm{T}$ in Section V.

Finally, to get ${\bm{A}}$, we exploit the shift invariant structure of $\mathcal{O}_{\alpha}$ w.r.t. ${\bm{A}}$, i.e.,

[TABLE]

where $\bm{J}_{\rm u}\mathcal{O}_{\alpha}$ and $\bm{J}_{\rm l}\mathcal{O}_{\alpha}$ denote the upper and lower $(\alpha-1)L$ blocks of $\mathcal{O}_{\alpha}$ [cf. (11)], respectively. By substituting then the expression (19) for $\mathcal{O}_{\alpha}$ into (21) and by using the estimate $\hat{\bm{T}}$ (20) for the transform matrix $\bm{T}$, the least squares estimate for ${\bm{A}}$ is given by

[TABLE]

Retrieving the input matrix ${\bm{B}}$. While the input matrix ${\bm{B}}$ can be obtained following a similar approach as for $\bm{A}$ [48] (yet with a more involved shift-invariant structure), we here compute it together with the initial state ${\bm{x}}(0)$ by solving a least squares problem.

To do so, we expand system (9) to all its terms as

[TABLE]

where $\bm{I}_{L}$ is the $L\times L$ identity matrix and vec$(\cdot)$ is the vectorization operator. We then collect the unknowns $\bm{x}(0)$ and ${\rm vec}(\bm{B})$ into the vector $\bm{\theta}=[\bm{x}(0)^{T}\,{\rm vec}(\bm{B})^{T}]^{T}$ and define the matrix

[TABLE]

where we substituted the state transition matrix $\bm{A}$ with its estimate (22) while the other quantities ${\bm{C}},{\bm{D}}$ and ${\bm{u}}(k)$ are known. Finally, we get the input matrix ${\bm{B}}$ by solving

[TABLE]

Given the system matrices $\{\hat{{\bm{A}}},\hat{\bm{B}},{\bm{C}},{\bm{D}}\}$, the state interaction graph ${\mathcal{G}}_{x}$ ($\bm{S}_{x}$) and input interaction graph ${\mathcal{G}}_{u}$ ($\bm{S}_{u}$) can be obtained by enforcing the constraints derived from the information of the physical process, i.e., the model dynamics. Hence, the network structure depends heavily on the estimate of the subspace span of $\mathcal{O}_{\alpha}$.

IV-A Noisy setting

We discuss here a method for estimating ${\bm{A}}$ and ${\bm{B}}$ with perturbations in the state evolution $\bm{x}(t)$ and noise in the observables $\bm{y}(t)$. To tackle these challenges, we lever instrumental variables.

Consider the following partition of the observables and input matrix

[TABLE]

where $\bm{Y}_{1}$ (resp. $\bm{U}_{1}$) and $\bm{Y}_{2}$ (resp. $\bm{U}_{2}$) have respectively $\beta$ and $\gamma$ blocks of size $L$ with $\gamma=\alpha-\beta$. The model for the observables $\bm{Y}_{2}$ is [cf.(13)]

[TABLE]

Following then the projection-based strategy to remove the dependency from $\bm{U}_{2}$, we can write the projected noisy observables as [cf. (16)]

[TABLE]

with $\bm{X}_{2}$ and $\bm{N}_{2}$ being the respective partitions of the state evolution and perturbation matrix.

Observe from (26) that the signal subspace is corrupted with the noise projection term $\bm{N}_{2}\bm{\Pi}_{\bm{U}_{2}^{T}}^{\perp}$. To remove the latter, we follow the instrumental variable approach. Since the noise is uncorrelated with the input ${\bm{U}}_{1}$ and since the noise present in the second block $\bm{N}_{2}$ is uncorrelated with the noise in the first block $\bm{N}_{1}$, it holds that

[TABLE]

Subsequently, we introduce the variable $\bm{Z}_{1}\triangleq[\bm{U}_{1}^{T},\,\bm{Y}_{1}^{T}]$ and consider the matrix

[TABLE]

to do the estimation of the signal subspace. The matrix $\bm{G}_{1}$ is asymptotically “noise free” due to (27). From the economy-size SVD ($N$-rank approximation) of $\bm{G}_{1}$, we get

[TABLE]

Finally, by using $\bm{W}_{\gamma,N}$, the system matrices ${\bm{A}}$ and ${\bm{B}}$ can be estimated from expressions (20), (22), and (24).

Note that the estimator for the signal subspace [cf.(29)] has intrinsic statistical properties. To keep under control the variance and bias, different works [49] have proposed to left and right weigh the matrices in (28) before the SVD. As establishing optimal weighting matrices requires further analysis, we do not detail them here and refer the interested reader to [50].

IV-B Continuous-time model identification

Given the discrete-time system matrices $\{\hat{{\bm{A}}},\hat{{\bm{B}}},{\bm{C}},{\bm{D}}\}$, we can estimate the continuous-time transition matrices of (3). Observe from (7) that estimating $\hat{f}_{x}(\bm{S}_{x})$ from $\hat{{\bm{A}}}$ requires computing only the matrix logarithm of $\hat{{\bm{A}}}$. Nevertheless, as stated by the following proposition, there are conditions process ${\mathcal{P}}$ should meet for this matrix logarithm to $(i)$ exist and $(ii)$ be unique. For clarity, the involved matrices are real since we have matrix functions that map real matrices onto real matrices.

Proposition 1**.**

Let the analytic function $f_{\rm s,x}(\cdot)$ of $f_{x}(\bm{S}_{x})$ in (7) satisfy

• (a)

$e^{f_{\rm s,x}(z)}\not\in\mathbb{R}^{-},\;\forall\;z\in{\rm eig}(\bm{S}_{x})$**

• (b)

$f_{\rm s,x}(z)>-\infty,\;\forall\;z\in{\rm eig}(\bm{S}_{x})$**

where $\mathbb{R}^{-}$ is the closed negative real axis. Then, process $\mathcal{P}$ guarantees that $(i)$ and $(ii)$ are met.

If the conditions of Proposition 1 are met for $f_{{\rm s},x}(\cdot)$, then ${\bm{A}}$ has no eigenvalues on $\mathbb{R}^{-}$. This implies that the principal logarithm111For two matrices $\bm{X}$ and $\bm{Y}$, $\bm{X}$ is said to be the matrix logarithm of $\bm{Y}$ if $e^{\bm{X}}=\bm{Y}$. If a matrix in invertible and has no-negative real eigenvalues, there is a unique logarithm which is called principal logarithm [42]. of ${\bm{A}}$, $\ln({\bm{A}})~{}=~{}\tau f_{x}(\bm{S}_{x})$, exists and is unique. Note that if ${\bm{A}}$ is real, its principal logarithm is also real. Hence, if $f_{{\rm s},x}(\cdot)$ satisfies the conditions of Proposition 1, the continuous-time transition matrices for nonsingular $\hat{\bm{A}}-\bm{I}$ are given by

[TABLE]

The expression for $\hat{f}_{u}(\bm{S}_{u})$ is derived from

[TABLE]

Given $\hat{f}_{x}(\bm{S}_{x})$ and $\hat{f}_{u}(\bm{S}_{u})$, we are left with the estimation of the underlying topologies. The network topology can be estimated with methods that promote sparse representations [15, 51], covariance methods [24, 28] or by solving a root finding problem, i.e., inverse maps. We analyzed this last case in the shorter version of this work [52].

However, for the case of partial observations, the similarity transform can not be uniquely identified making the above methods not applicable for identifying the network topology due to the ambiguities in the system matrices. Therefore, in the next section, we discuss such ambiguities and introduce a AP method for estimating the network topology.

Remark 1**.**

Although we focused principally on continuous-time models, all stated results hold also for purely discrete models by doing the appropriate minor changes for the functional dependencies of the system matrices.

V Partially Observed Network

So far, we considered the observables ${\bm{y}}(t)$ are available for the whole network, i.e., $\bm{C}=\bm{I}$ or more generally ${\rm rank}(\bm{C})=N$. The latter allows a unique estimate for the similarity matrix $\bm{T}$ in (20). We here move to the more involved case where it is not possible to observe the process on all nodes or the observations are not sufficiently rich, i.e., ${\rm rank}({\bm{C}})<N$. This disables to find a unique transform matrix $\bm{T}$ and, instead of retrieving the original system matrices, we retrieve an estimate of a set of equivalent matrices

[TABLE]

which also realize the system in (10). Further, it follows from (32) that (although ${\bm{A}}_{\rm T}\neq{\bm{A}}$ in general) if ${\bm{A}}$ is diagonalizable, i.e., ${\bm{A}}=\bm{Q}_{A}\bm{\Lambda}_{A}\bm{Q}_{A}$, then

[TABLE]

holds, where ${\rm eig}({\bm{A}})$ are the eigenvalues of ${\bm{A}}$. The equality (33) yields since ${\bm{A}}$ and ${\bm{A}}_{\rm T}$ are similar matrices222Two matrices $\bm{X}$ and $\bm{Y}$ are said to be similar if there exists an invertible matrix $\bm{P}$ such that $\bm{X}=\bm{P}\bm{Y}\bm{P}^{-1}$..

In these situations, we cannot remove the ambiguity in the system matrices without additional information. In the sequel, we motivate why this disambiguation problem is particularly hard. We further derive a method to estimate an approximately feasible realization of the network topology related to the signal subspace and to the knowledge of the (bijective) scalar mappings $\{f_{{\rm s},x}(\cdot),f_{{\rm s},u}(\cdot)\}$.

Finally, we stress that despite that in the following part we will concentrate in recovering the state network topology $\mathcal{G}_{x}$, this is not a problem for cases when the input network topology $\mathcal{G}_{u}$ is of interest. For such cases, an additional step can always be performed to retrieve $\bm{S}_{u}$ from the transformed matrix $\bm{B}_{T}$. For instance, from the estimate of $\bm{A}$ based on the state network topology and $\bm{A}_{T}$, the transform $\hat{\bm{T}}$ can be estimated. Subsequently, the inverse operation is applied to ${\bm{B}}_{T}$ to get $\hat{\bm{B}}$ and, therefore, $\mathcal{G}_{u}$. Hence, there is no detriment on focusing in recovering $\mathcal{G}_{x}$.

V-A The graph inverse eigenvalue problem

To start, consider the shift operator $\bm{S}_{x}$ belongs to a set $\mathcal{S}$ which contains all permissible matrices representing the dynamics of the state that lead to ${\bm{A}}_{\rm T}$. The set $\mathcal{S}$ describes the properties of the graph representation matrix, e.g., zero diagonal (adjacency) $[\bm{S}_{x}]_{n,n}=0\,\forall\,n\in[N]$, unitary diagonal (normalized Laplacian) $[\bm{S}_{x}]_{n,n}=1\,\forall\,n\in[N]$, zero eigenvalue related to the constant eigenvector (combinatorial Laplacian) $\bm{S}_{x}\bm{1}=\bm{0}$, symmetry (undirected graphs) $\bm{S}_{x}=\bm{S}_{x}^{T}$.

The ambiguity (33) introduced by the similarity transform, $\bm{T}$, transforms the problem of finding $\bm{S}_{x}$ into

[TABLE]

where $\bm{\lambda}_{x}={\rm eig}(\bm{S}_{x})$ is the vector containing the eigenvalues of $\bm{S}_{x}$ obtained by applying the inverse map to the eigenvalues of $\bm{A}_{T}$. Problem (34) recasts the network topology identification problem to the problem of finding a graph shift operator matrix $\bm{S}_{x}$ that has a fixed spectrum. This problem belongs to the family of inverse eigenvalue problems [53]. In a way, problem (34) is the complement of the network spectral template approach [15]. Here, instead of having an eigenbasis and searching for a set of eigenvalues, we have a set of eigenvalues and search for an eigenbasis.

Problem (34) is ill-posed as, in most cases, its solution is non-unique. Thus, it leads to ambiguities in its solution. In what follows, we characterize this ambiguity in terms of equivalence classes between graphs and provide a method able to find a network topology satisfying the conditions of (34).

V-B Ambiguous graphs: cospectral graphs

We find ambiguities in graph topologies in graphs that belong to an equivalence class [54]. An example of equivalent graphs, in terms of their spectrum, are the isomorphic graphs [34]. For unlabeled graphs, two graphs $\mathcal{G}$ and $\mathcal{G}^{\prime}$ with respective graph shift operator matrices $\bm{S}$ and $\bm{S}^{\prime}$ are equivalent if there exists a permutation matrix $\bm{P}$ such that

[TABLE]

That is, the graph representation matrices are row- and column-permuted versions of each other. The permutation matrix $\bm{P}$ implements the isomorphism. Therefore, if the only information available is the graph eigenvalues $\bm{\lambda}$, the graphs are always indistinguishable up to node reordering. This situation is not at all undesirable as the ordering of the nodes is often not important. However, isomorphic graphs are not only ones sharing the spectrum.

Graphs that share the spectrum are called cospectral (or isospectral) graphs [55]. Note, however, that cospectral graphs are not necessarily isomorphic. Fig. 1 illustrates an example of two cospectral graphs. Graph cospectrallity renders the feasible set of (34) not a singleton and, therefore, we need to settle with any feasible graph satisfying the constraints. That is, the identified topology from (34) will be a graph within the equivalence class of cospectral graphs regarding $\bm{\lambda}$ and $\mathcal{S}$. The following definition formalizes the latter.

Definition 1**.**

(Equivalent cospectral graphs)* Two graphs $\mathcal{G}$ and $\mathcal{G}^{\prime}$ of $N$ nodes are cospectral equivalents with respect to the spectrum $\bm{\lambda}$ and the graph representation set $\mathcal{S}$, if they belong to the set*

[TABLE]

As a result, the feasibility problem (34) reduces to a graph construction problem (graph inverse eigenvalue problem). Put differently, the problem at hand can be rephrased as given a spectrum $\bm{\lambda}$ and a set $\mathcal{S}$, construct a graph shift operator matrix $\bm{S}\in\mathcal{S}$ with spectrum $\bm{\lambda}$. We shall discuss next a method that addresses this construction problem.

V-C Graph construction by alternating projections

Before detailing the graph construction method, we introduce the following assumptions.

1. (A.1)

The set $\mathcal{S}$ is closed. 2. (A.2)

For any $\bm{S}\in\mathbb{R}^{N\times N}$, the projection $P_{\mathcal{S}}(\bm{S})$ of $\bm{S}$ onto the set $\mathcal{S}$ is unique.

The first assumption is technical and guarantees that the set $\mathcal{S}$ includes all its limit points. The second assumption is slightly more restrictive and ensures that the problem

[TABLE]

has a unique solution. Although this assumption might seem restrictive, in most cases we only have access to a convex description of the feasible set $\mathcal{S}$ which satisfies A.2 (as the set is assumed closed) or only the projection onto the convex approximation of the feasible set of network matrices can be performed efficiently. Thus, it is fair to consider that A.2 holds in practice.

For the sake of exposition, we restrict our upcoming discussion to the case of symmetric matrices, i.e., undirected graphs but remark that a similar approach can be followed for directed graphs333This could be done by exchanging the spectral decomposition for the Schur decomposition [56] which decomposes a matrix into unitary matrices and an upper triangular matrix.. Further, denote by $\mathcal{S}^{N}$ the set of symmetric $N\times N$ matrices and by $\mathcal{S}_{+}^{N}$ the set of positive semidefinite matrices. We then recall the following result from [57, Thm. 5.1].

Theorem 1** (adapted).**

Given $\bm{S}\in\mathcal{S}^{N}$ and let $\bm{S}=\bm{Q}\bm{\Lambda}\bm{Q}^{T}$ be the spectral decomposition of $\bm{S}$ with non-increasing eigenvalues $[\bm{\Lambda}]_{ii}\geq[\bm{\Lambda}]_{jj}$ for $i<j$. For a fixed $\bm{\Lambda}_{o}\in\mathcal{D}_{N}$ with non-increasing elements $[\bm{\Lambda}_{o}]_{ii}\geq[\bm{\Lambda}_{o}]_{jj}$ for $i<j$, a best approximant, in the Frobenius norm sense, of $\bm{S}$ in the set

[TABLE]

*is given by *

[TABLE]

where ${\rm O}(N)$ denotes the set of the $N\times N$ orthogonal matrices.

Theorem 1 implies that the projection of $\bm{S}$ onto $\mathcal{M}$ is not necessarily unique. As an example, consider the graph shift operator $\bm{S}$ with repeated eigenvalues. Here, it is not possible to uniquely define a basis for the directions related to the eigenvalues with multiplicity larger than one. Hence, infinitely many eigendecompositions exist that lead to many projections of $\bm{S}$ onto $\mathcal{M}$. Further, note that as every element of $\mathcal{M}$ is uniquely determined by an element of ${\rm O}(N)$, the structure of $\mathcal{M}$ is completely defined by the structure of ${\rm O}(N)$. Therefore, as ${\rm O}(N)$ is a smooth manifold, $\mathcal{M}$ is one as well.

Alternating projections method. As it follows from Assumptions A.1 and A.2 and Theorem 1, we can project any graph shift operator matrix $\bm{S}\in\mathcal{S}^{N}$ onto $\mathcal{S}$ and $\mathcal{M}$. Further, by noticing that the construction problem (34) is equivalent to finding a matrix in

[TABLE]

we can consider the simple, and intuitive, method of alternating projections (AP) [35] to find a point in (37). The AP method finds a point in the intersection of two closed convex sets by iteratively projecting a point onto the two sets. It performs the updates

[TABLE]

starting from a point $\bm{S}_{0}\in\mathcal{M}$. The AP method has guaranteed convergence for convex sets and it does that linearly. However, for alternating projections on a combination of different types of sets (we have a set $\mathcal{S}$ satisfying A.1 and A.2, and a smooth manifold $\mathcal{M}$), additional conditions on both sets are necessary to guarantee convergence.

First, let us formalize the notion of a fixed point for the iterative procedure (38).

Definition 2** (Fixed point).**

A matrix $\bm{S}\in\mathcal{S}^{N}$ is a fixed point of the alternating projections procedure in (38) if there exists an eigendecomposition of $P_{\mathcal{S}}(\bm{S})$,

[TABLE]

with non-increasing elements $[\bm{\Lambda}]_{ii}\geq[\bm{\Lambda}_{jj}]$ if $i<j$ such that

[TABLE]

This definition makes explicit two things. First, it defines $\bm{S}$ as a fixed point if and only if

[TABLE]

Second, for the cases where $\bm{S}$ is a fixed point and $P_{\mathcal{S}}(\bm{S})$ has eigenvalues with multiplicity larger than one, progress can be still made towards a feasible solution when $P_{\mathcal{M}}(P_{\mathcal{S}}(\bm{S}))\not\in\mathcal{S}\cap\mathcal{M}$. To see the latter, consider the case where an alternative eigendecomposition

[TABLE]

is available for the fixed point $\bm{S}$. Assuming that

[TABLE]

with $\bm{Q}$ the eigenbasis that makes $\bm{S}$ a fixed point, we can see that the new point $\tilde{\bm{S}}$ escapes from the fixed point. Further, since successive projections between two closed sets is a nonincreasing function over the iterations [58, Thm. 2.3], we can prove that $\tilde{\bm{S}}$ presents a progress towards a feasible solution in $\mathcal{S}\cap\mathcal{M}$.444 We considered the notion of fixed point to obtain a feasible set of the system matrices (matrices that realize the system) since, beyond their structure, the most important characteristic is their spectrum (set of eigenvalues).

With this in place, the following theorem shows that the AP method converges when used for the graph construction problem.

Theorem 2**.**

Let $\mathcal{S}$ meet Assumptions A.1-A.2 and consider the set $\mathcal{M}$ defined in Theorem 1. Let also $\bm{S}_{0},\bm{S}_{1},\bm{S}_{2},\ldots,$ be a sequence generated by the alternating projections method in (38). Then, there exists a limit point $\bm{S}$ of this sequence that is a fixed point of (38) [cf. 2] satisfying

[TABLE]

Further, if the limit is zero, then $\bm{S}\in\mathcal{S}\cap\mathcal{M}$.

Proof.

See Appendix IX-A555Appendix in the supplemental material.. ∎

The above theorem proves that the AP method retrieves a matrix $\bm{S}\in\mathcal{M}$ that realizes the original system, i.e., it preserves the underlying system modes and is an approximately feasible network representation. Nevertheless, Theorem 2 does not quantify the rate of convergence of such a method. By particularizing results for super-regular sets [59, Thm. 5.17], the following theorem shows that if the problem is feasible, locally, the proposed method converges linearly to a point in (37).

Theorem 3**.**

Let the set of all permissible matrices $\mathcal{S}$ [cf. (34)] be convex and meet Assumptions A.1-A.2. Let also the set $\mathcal{M}$ be defined as in Theorem 1. Denote by $N_{\mathcal{S}}({\bm{S}})$ the normal cone of the closed set $\mathcal{S}$ at a point $\bm{S}$ and, similarly, by $N_{\mathcal{M}}({\bm{S}})$ the normal cone of the set $\mathcal{M}$ at $\bm{S}$. Further, suppose that $\mathcal{M}$ and $\mathcal{S}$ have a strongly regular intersection at $\bar{\bm{S}}$, i.e., the constant

[TABLE]

is strictly less than one with $B$ being a closed unit Euclidean ball. Then, for any initial point $\bm{S}_{0}\in\mathcal{M}$ close to $\bar{\bm{S}}$, any sequence generated by the alternating projections method in (38) converges to a point in $\mathcal{M}\cap\mathcal{S}$ with $R$-linear rate

[TABLE]

Proof.

See Appendix IX-B. ∎

These results guarantee that the AP method converges globally (at least) to a fixed point in $\mathcal{M}$ and locally, i.e., within the neighborhood of the solution (if it exists), to a fixed point in $\mathcal{M}\cap\mathcal{S}$.

V-D Inaccurate and partial information

We now consider the case where the estimated eigenvalues are inexact because of noise or are incomplete because the full eigendecomposition of the system matrices is not feasible. To deal with such cases, we modify the structure of the set $\mathcal{M}$ (in Theorem 1) to reflect the uncertainty and the partial eigendecomposition. The modified set has to be compatible with the structure used in Definition 2 and Theorems 2 and 3 to guarantee the convergence of the AP method. Thus, in the sequel, we focus on proving compactness for the modified versions of $\mathcal{M}$, which suffices to guarantee convergence of the AP method by the result of Theorem 2.

Uncertainty in the system matrices. The following proposition shows that the set $\mathcal{M}_{\epsilon}$, which allows the estimated eigenvalues to lie within an $\epsilon-$uncertainty ball, is compact.

Proposition 2**.**

*Let $\bm{\Lambda}_{o}\in\mathcal{D}_{N}$ with $[\bm{\Lambda}_{o}]_{ii}\geq[\bm{\Lambda}_{o}]_{jj}$ for $i<j$ be fixed. If $0\leq\epsilon<\infty$ is a fixed scalar accounting for uncertainties on the elements of $\bm{\Lambda}_{o}$, then the set *

[TABLE]

is compact.

Proof.

See Appendix IX-C. ∎

The following result provides a best approximant of a matrix $\bm{S}$ in the set $\mathcal{M}_{\epsilon}$ in the Frobenius norm.

Theorem 4**.**

Given $\bm{S}\in\mathcal{S}^{N}$ with eigendecomposition $\bm{S}=\bm{Q}\bm{\Lambda}\bm{Q}^{T}$ and non-increasing eigenvalues $[\bm{\Lambda}]_{ii}\geq[\bm{\Lambda}]_{jj}$ for $i<j$. For a fixed $\bm{\Lambda}_{o}\in\mathcal{D}_{N}$ with $[\bm{\Lambda}_{o}]_{ii}\geq[\bm{\Lambda}_{o}]_{jj}$ for $i<j$, a best approximant of $\bm{S}$ in $\mathcal{M}_{\epsilon}$, in the Frobenius norm sense, is given by

[TABLE]

where

[TABLE]

Proof.

See Appendix IX-D. ∎

Corollary 1**.**

The nonzero entries of $\bm{\Lambda}_{\epsilon}^{*}$ are

[TABLE]

with $\gamma_{i}:=[\bm{\Lambda}]_{ii}-[\bm{\Lambda}_{o}]_{ii}$.

Partial eigendecomposition. In physical systems, a discrete model of $N$ degrees of freedom provides accurate information of about $N/3$ of the system natural frequencies [53, Ch. 5]. In other cases, the full eigendecomposition of the system matrix is not always possible. We, therefore, provide a projection onto a set that only considers a noisy part of the system matrix spectrum is available.

The following theorem provides the main result.

Theorem 5**.**

Let $\bm{\Lambda}_{m}\in\mathcal{D}_{m}$ with $[\bm{\Lambda}_{m}]_{ii}\geq[\bm{\Lambda}_{m}]_{jj}$ for $i<j$ be fixed. If $0\leq\epsilon<\infty$ is a fixed scalar accounting for uncertainties on the elements of $\bm{\Lambda}_{m}$ and $\rho:=\max_{\bm{S}\in\mathcal{S}}\|\bm{S}\|_{2}$, then a best approximant, in the Frobenius norm sense, of $\bm{S}\in\mathcal{S}^{N}$, with $\|\bm{S}\|_{2}\leq\rho$, in the set

[TABLE]

is given by

[TABLE]

Here, $\bm{\sigma}$ denotes the permutation of the subset of $[N]$ that solves the combinatorics problem

[TABLE]

where $\bm{\Lambda}$ is the diagonal matrix of eigenvalues of $\bm{S}$, and $\bar{\bm{\sigma}}$ is the complementary set of $\bm{\sigma}$. The matrix $\bm{Q}$ is given by the (sorted) eigendecomposition of $\bm{S}$, i.e.,

[TABLE]

where $\bm{\Lambda}_{\bm{\sigma}}$ is the permuted version of $\bm{\Lambda}$ and

[TABLE]

Furthermore, $\mathcal{M}_{\epsilon}^{m}$ is compact.

Proof.

See Appendix IX-F. ∎

Corollary 2**.**

The optimal permutation $\bm{\sigma}$ of the indices $[N]$ that solves (44) can be found by solving a minimum-weight bipartite perfect matching problem.

The above results in Theorems 4 and 5 showed that the sets $\mathcal{M}_{\epsilon}$ and $\mathcal{M}_{\epsilon}^{m}$ are compact and provide a best Frobenius-norm approximant for $\bm{S}$ in each case. Therefore, we can apply the AP method in (38) to these scenarios using the appropriate modifications. Finally, since the sets $\mathcal{M}_{\epsilon}$ and $\mathcal{M}_{\epsilon}^{m}$ meet the conditions of Theorem 2, the convergence results for the AP method extend also to the discussed scenarios.

VI System Consistency Constraints

The set $\mathcal{S}$ [cf. (34)] plays an important role in the system topology that the AP method identifies. As such, it should be constrained such that the AP method yields a consistent system, i.e., the AP estimated state network should define an equivalent system to the original one. We briefly discuss here two constraints that can be added to $\mathcal{S}$ to enforce system consistency.

Given that the system matrix [cf. (7)] is a bijective matrix function and by using the same construction as for the shift invariance property in (21), we can build the linear system

[TABLE]

Here, we leverage the invariance of the matrix function to nonsingular transforms, i.e.,

[TABLE]

where $f(\cdot)$ is a matrix function, hence, can be applied to $\bm{A}_{T}$. In this way, we get a linear system that depends on $\bm{S}_{x}$ and enforces the shift invariance condition. Nevertheless, the shift invariance condition does not change the optimality nor the uniqueness of the projection onto $\mathcal{S}$. This is because $\bm{T}$ is a free optimization variable and does not affect the projection distance [cf. (36)].

If other constraints for the transform matrix $\bm{T}$ are known, they can be included when solving for the projection (36). These additional constraints will not impact the projection optimality because they do not change the cost function for $P_{\mathcal{S}}$. For instance, consider the constraint that requires symmetry in the input matrix $\bm{B}$. Then, since $\bm{B}_{T}=\bm{T}\bm{B}$, $\bm{B}$ is symmetric, and $\bm{B}_{T}\bm{T}=\bm{T}\bm{B}\bm{T}^{T}$, we have

[TABLE]

Again, such an additional constraint does not change the projection distance as it only modifies the description of the convex set in which the matrix must be projected.

We can introduce other constraints to the set $\mathcal{S}$ to further restrict the family of feasible network representations. However, these constraints should be analyzed as case-specifc and go beyond the main goal of this work. Next, we corroborate the above theoretical findings with numerical results.

VII Numerical Results

In this section, we present a series of numerical results to illustrate the performance of the proposed methods for different scenarios. We first illustrate how the model and the noise coloring influence the estimation performance of commonly used topology identification methods. Then, we corroborate our theoretical results and, finally, we present results for the topology identification from partially-observed networks.666The code to generate these numerical results can be found in https://gitlab.com/fruzti/systemid_codes

VII-A Discrete model validation

In this section, we corroborate the discrete model (6) in finding a graph topology from continuous-time data generated following the model (3). The underlying graphs are two fixed random regular graphs of $N=50$ nodes with node degree of $d=3$. The data are generated by a continuous-time solver with system evolution matrix $f_{x}(\bm{S}_{x})=-\bm{S}_{x}$ and input matrix $f_{u}(\bm{S}_{u})=-(\bm{S}_{u}+\bm{I})$. The observable matrix $\bm{C}$ is set to identity and $\bm{D}$ is the zero matrix. The input signal is drawn from a standard normal distribution and we set the number of samples to $N^{3}$ with a sampling time of $\tau=10^{-3}$.

We compare the proposed method with the spectral templates techniques in [15]. For the latter, the system matrices $\bm{A}$ and $\bm{B}$ are first obtained from the continuously sampled data (cf. Section III). Then, the eigenvectors of these matrices are used as spectral templates. These results are shown in Fig. 2.

Fig. (2(a)) shows the estimated network topologies for a particular input signal realization, while Fig. (2(b)) compares the respective spectra. We observe that the approach relying on spectral templates overestimates the number of edges and underestimate the graph eigenvalues. However, the graph obtained with spectral templates is a matrix function of the original graph, i.e., there is a function (polynomial) that maps the estimated graph to the original one. This is because $\bm{S}_{x}$ has all eigenvalues with multiplicity one. The proposed technique relying on the discrete model (6) retrieves the eigenvalues and the graph structure perfectly. This result is not surprising since subspace-based system identification is a consistent estimator for the transition matrix and the proposed method use the knowledge of $f_{s,x}$ while the spectral templates does not. For this scenario, we also considered building the graph from the data covariance matrix, but this technique did not lead to satisfactory results. We attribute this misbehavior to the fact that the covariance matrix is not diagonalizable by the graph modes (i.e., eigenvectors of $\bm{S}_{x}$) due to the prescence of the input.

VII-B Instrumental variables approach on social graph

We now evaluate the instrumental variable approach of Section IV-A on the Karate club graph [60]. The graph represents the connections of $N=34$ members through $78$ undirected edges. We consider the discrete model (6) with $\bm{A}$ described by a continuous-time diffusion process $\bm{A}=e^{-\tau\bm{L}_{x}}$. The diffusion rate (or sampling time) is fixed to $\tau=10^{-3}$. The input signal is randomly generated from a standard normal distribution and the power of both the state and the observation noise is $\sigma^{2}=10^{-3}$

We consider three different approaches to estimate the underlying network topology: $i)$ a covariance-based approach, where the covariance matrix is estimated from the observables; $ii)$ the instrumental variable approach combined with the spectral template method from [15]; and $iii)$ the proposed instrumental variable approach by enforcing the dynamics of the continuous system. These results are reported in Fig. 3.

Fig. 3(a) shows the fitting accuracy of the subspaces, while Fig. 3(b) illustrates the fitting of the eigenvalues. In Fig. 3(c), we provide the obtained graphs where the edges with absolute weight less than $10^{-3}$ are omitted. We observe that the system identification flow allows a better graph reconstruction and the proposed method offer the best alignment of the eigenbasis. Further, note that by levering the underlying physical model of the diffusion, we can reconstruct the graph spectrum with high fidelity. We further remark that despite both the basis and the spectrum are aligned, the retrieved graph looks different from the true one. This is because of the ambiguities discussed in Section V-B. However, the obtained graph has the same eigenvalues as the original one and its basis diagonalizes the original network matrix. Notice also that from the three methods, only the one leveraging the model information retrieves a connected graph after thresholding.

Finally, we remark that the task of estimating this topology based purely on a spectral decomposition is hard. This is because the combinatorial Laplacian of the Karate club graph has eigenvalues with a multiplicity larger than one. Thus, there is no unique basis for its eigendecomposition leading to difficulties in reconstructing the underlying topology. This issue is also present even when the topology is the sparsest matrix that generates such dynamics.

VII-C Convergence of the alternating projections method

We analyze here the convergence behavior of the alternating projections method (38). We present results using the sets $\mathcal{M}_{\epsilon}^{m}$ and $\mathcal{S}=\mathcal{L}_{{\rm CVX}}$. The latter is the convex relaxation of the combinatorial Laplacian set. These sets are chosen to illustrate the convergence results as $\mathcal{L}_{{\rm CVX}}$ encompasses the problem of finding Laplacian matrices with given eigenvalues and $\mathcal{M}_{\epsilon}^{m}$ is the most general set proposed in this work. Additional results for the other sets are provided in the supplementary material. For this scenario, we select a regular graph777See supplemental material with $N=30$ and node degree $d=3$ and consider only half of its eigenvalues known, i.e, $m=N/2$. The AP method is analyzed for five different initial points.

These results are shown in Fig. 4. Here, each solid line represents a different starting point. The (blue) dashed line shows the convergence behavior when the starting point is the (diagonal) eigenvalue matrix. These results show two main things. First, the predicted monotone behavior of the error $\|\bm{S}_{k}-P_{\mathcal{S}}(\bm{S}_{k})\|_{{\rm F}}$ holds and stagnates when a limit point is reached by the iterative sequence. Second, the error $\|\bm{S}_{k}-\bm{S}_{k+1}\|_{{\rm F}}$ converges to the desired accuracy (order $10^{-6}$), although not monotonically, The error convergence rate is generally linear and the starting point influences the slope. Finally, we emphasize that even when the set of known eigenvalues lies within an $\epsilon$-ball of uncertainty, the alternating projections method convergences. The convergence is guaranteed by the compactness of the set $\mathcal{M}_{\epsilon}^{m}$.

VII-D Partial observations

In this section, we consider the task of retrieving a graph that realizes a given system from partial observations. We consider two regular random graphs of $N=14$ nodes and three edges per node. The data are generated from a continuous diffusion on the network and the input matrix is $\bm{B}=\bm{L}_{u}+\bm{I}$. The matrix $\bm{C}$ is a Boolean matrix that selects half of the nodes (the odd labeled nodes from arbitrary labeling). Note that none of the previous methods can be employed to retrieve the network topology since $\bm{B}\neq\bm{I}$ and the network is not fully observed. Even if the covariance matrix is estimated from sampled data, its eigenstructure does not represent the eigenstructure of the state network topology.

We first estimate the system matrices using the system identification framework and then employ the AP method initialized with a random symmetric matrix that has as eigenvalues the estimated state network eigenvalues. The constraint set in (34) is the convex relaxation of the combinatorial Laplacian set [15]. We further enrich this set with the system identification constraints to enforce the feasibility of the realization. Fig. 5 reports the results after $30$ iterations of the AP method.

From Fig. 5(a), we observe that the estimated state graph does not exactly share the eigenbasis with the original one, i.e., the graph mode projections do not form a diagonal matrix. However, we could perfectly match the input graph. This behavior is further seen in the eigenvalues, where those of the input graph are matched by the estimated eigenvalues. For the state graph, a perfect eigenvalue match is possible if a final projection onto $\mathcal{M}$ is performed. These results are also reflected in the estimated topologies in Fig. 5. Perfect reconstruction of the input graph support is achieved, while the state graph presents a different arrangement in the nodes and it is not regular.

Despite the differences in the state graphs, the triple $\{\hat{\bm{S}}_{x},\hat{\bm{S}}_{u},\bm{C}\}$ realises (approximately) the same system as the triple $\{\bm{S}_{x},\bm{S}_{u},\bm{C}\}$. This is because the product of the involved system matrices is preserved, i.e., although the structure of the state graph is different, the observations can be reproduced with high confidence using the estimated system matrices. With the estimated graphs, we can predict the system output with an NRMSE fitness of $\approx 95\%$. In the supplementary material, we compare the response of the true system to an arbitrary excitation with the obtained system response that uses the estimated graphs.

VIII Conclusions

This paper introduced a general framework for graph topology identification through state-space models and subspace techniques. We showed that it is possible to retrieve the matrix representation of the involved graphs from the system matrices by exploiting the geometric structure of the input and output data. In particular, we discussed the challenges of retrieving the network topology under partial observations and proposed an alternating projections method to recover a set of matrices that realizes the system and are conform with user-defined constraints. The proposed theoretical analysis is corroborated with numerical results. Future research is needed in three main directions. First, the focus should be on improving the scalability of the proposed techniques to larger graphs. Second, research is needed in employing subspace models to learn a coarser graph that drives the system dynamics in large data sets. Third, extensions of the current approach to nonlinear systems, such as the ones in chemical reaction networks are also worth being investigated.

IX Appendix

IX-A Proof Theorem 2

We construct the proof of the convergence of the method with arguments similar to  [58].

Note that any matrix $\bm{S}_{k}\in\mathcal{M}$ can be represented by a matrix $\bm{Q}\in{\rm O}(N)$. Then, since ${\rm O}(N)$ is compact, yields the set $\mathcal{M}$ is compact. The latter implies $\mathcal{M}$ contains a convergent sequence and that $\bm{S}_{k}$ has a limit point denoted by $\bm{S}^{*}$. Further, consider the implementation of (38)

[TABLE]

and

[TABLE]

From the compactness of ${\rm O}(N)$, $\bm{Q}_{k}$ converges to a point $\bm{Q}$ that, in turn, makes $\bm{S}_{k+1}$ converge to $\bm{S}$.

Further, since the successive projections between $\mathcal{S}$ and $\mathcal{M}$ is a nonincreasing function over the iterations [58, Thm. 2.3], the limit $\lim_{k\rightarrow\infty}\|\bm{S}_{k}-P_{\mathcal{S}}(\bm{S}_{k})\|$ exists. Then, since $P_{\mathcal{S}}(\cdot)$ is a continuous operation, we have that

[TABLE]

Then, since $P_{\mathcal{S}}(\bm{S}_{k})$ converges to $P_{\mathcal{S}}(\bm{S}^{*})$ and the orthogonal matrix $\bm{Q}_{k}$ converges to $\bm{Q}$, it follows that $P_{\mathcal{S}}(\bm{S}^{*})=\bm{Q}\bm{\Lambda}^{*}\bm{Q}^{T}$ for some $\bm{\Lambda}^{*}$. Further, the eigendecomposition $\bm{S}=\bm{Q}\bm{\Lambda}_{o}\bm{Q}^{T}$ implies that $\bm{S}=P_{\mathcal{M}}(P_{\mathcal{S}}(\bm{S}^{*}))$.

From (48) and since $\bm{S}$ is the projection of $P_{\mathcal{S}}(\bm{S}^{*})$ onto $\mathcal{M}$, we have

[TABLE]

By considering then that the projection of $\bm{S}$ onto $\mathcal{S}$ is unique (Assumption A.2), we have $P_{\mathcal{S}}(\bm{S}^{*})=P_{\mathcal{S}}(\bm{S})$. Further, since $\bm{S}$ is a projection of $P_{\mathcal{S}}(\bm{S}^{*})$ onto $\mathcal{M}$, $\bm{S}$ is a fixed point.

Finally, consider that the limit in (48) is zero. An arbitrary sequence generated by the AP method converges to some point $\tilde{\bm{S}}\in\mathcal{M}$, which is a limit point. From the inequality

[TABLE]

and by taking the limit for $k\rightarrow\infty$, we observe that $\tilde{\bm{S}}$ is also a limit of points in $\mathcal{S}$. As both sets are closed by assumption, the result follows.

IX-B Proof Theorem 3

Consider that $\mathcal{S}$ and $\mathcal{M}$ are a convex set and a smooth manifold, respectively. Therefore, they are both super-regular sets [59], which implies that $\mathcal{S}$ and $\mathcal{M}$ are super-regular sets at $\bar{\bm{S}}$. By the results from [59, Thm. 5.15] and [59, Thm. 5.17] that guarantee the convergences of alternating projections in super-regular sets (under the same condition at the intersection as in this theorem), the convergence guarantee follows.

IX-C Proof Proposition 2

Consider that the set

[TABLE]

is compact. This holds since the set $D_{N,e}$ is equivalent to the set of diagonal matrices whose diagonal entries lie in $[-\epsilon,\epsilon]$. Hence, the compactness.

Consider now the map $\phi:{\rm O}(N)\times\mathcal{D}_{N,\epsilon}\rightarrow\mathcal{M}_{\epsilon}$ given by

[TABLE]

which defines the set $\mathcal{M}_{\epsilon}$ and observe that ${\rm O}(N)$ is compact. The map $\phi$ is continuous and surjective, thus the set $\mathcal{M}_{\epsilon}$ is compact.

IX-D Proof Theorem 4

The proof starts by expressing the set $\mathcal{M}_{\epsilon}$ as

[TABLE]

where $\mathcal{M}(\bm{\Lambda}_{\epsilon})$ is defined similarly to the set $\mathcal{M}$ [cf. Theorem 1] by substituting $\bm{\Lambda}_{o}$ with $\bm{\Lambda}_{o}+\bm{\Lambda}_{\epsilon}$. Then, we can expand the projection problem for $\bm{S}$ as

[TABLE]

To obtain the first equality, we used the union description of $\mathcal{M}_{\epsilon}$ in (49). The second and third equalities follow from Theorem 1. The last equality follows from the invariance of the Frobenius norm under unitary transforms. Finally, observe that the non-increasing ordering of eigenvalues for both $\bm{\Lambda}$ and $\bm{\Lambda}_{o}$ solves the last optimization problem, thus providing the best Frobenius-norm approximant for $\bm{S}$.

IX-E Enforcing Positive semidefinitness under noise

Observe that the set $\mathcal{M}_{\epsilon}$ [cf. Prop. 2]

[TABLE]

specifies a neighborhood where the eigenvalues of $\bm{M}$ must lie. However, if the desired operator $\bm{M}$ is positive semidefinite (PSD), the constraint in the norm $\|\bm{\Lambda}_{\epsilon}\|_{2}$ might lead to matrices that are not PSD. The following proposition introduces an alternative set $\mathcal{M}_{\bm{a},\bm{b}}$ that guarantees that $\bm{M}$ is PSD and shows that this new set is also compact.

Proposition 3**.**

Consider the set of $N\times N$ diagonal matrices $\mathcal{D}_{N}$. Let $\bm{\Lambda}_{o}\in\mathcal{D}_{N}$ be a fixed diagonal matrix with non-increasing eigenvalues $[\bm{\Lambda}_{o}]_{ii}\geq[\bm{\Lambda}_{o}]_{jj}$ for $i<j$. Let also $\bm{\Lambda}_{\epsilon}\in\mathcal{D}_{N}$ be a diagonal matrix accounting for errors on the elements of $\bm{\Lambda}_{o}$ with bounded elements $[\bm{\Lambda}_{\epsilon}]_{ii}\in[a_{i},b_{i}]$.

Then, the set

[TABLE]

is compact for vectors $\bm{a}$ and $\bm{b}$ such that $[\bm{a}]_{i}=a_{i}$ and $[\bm{b}]_{i}=b_{i}$.

Proof.

The proof follows similarly to that of Proposition 2 for the set $\mathcal{M}_{\epsilon}$.

Consider the map

[TABLE]

which defines the set $\mathcal{M}_{\bm{a},\bm{b}}$ and observe that ${\rm O}(N)$ and all intervals $[a_{i},b_{i}]$ are compact.

Further, it can be shown that the map $\phi$ is continuous and surjective, thus the set $\mathcal{M}_{\bm{a},\bm{b}}$ is compact. Alternatively, consider that the Cartesian product of compact sets is compact (under the appropriate topology) leading to a compact set $\mathcal{M}_{\bm{a},\bm{b}}$. ∎

IX-F Proof Theorem 5 (Sketch.)

We first show that the set $\mathcal{M}_{\epsilon}^{m}$ is compact. Consider the map

[TABLE]

which defines the set $\mathcal{M}_{\epsilon}^{m}$ and observe that ${\rm O}(N)$, $\mathcal{D}_{m,\epsilon}$, and $\mathcal{D}_{N-m,\rho}$ are compact. Then, since the Cartesian product of compact sets is compact (under the appropriate topology), the set $\mathcal{M}_{\epsilon}^{m}$ is compact.

To prove the optimality of the provided best approximant, we proceed as follows. First, following the construction of [61, 62], it can be shown that for finding the shortest distance from $\bm{S}$ to $\mathcal{M}_{\epsilon}^{m}$, it suffices to find the shortest distance to a particular substructure of $\mathcal{M}_{\epsilon}^{m}$ defined by the permutation $\bm{\sigma}$. Then, following the same construction of Theorem 4, it can be proven the structure of the optimal $\bm{\Lambda}_{\epsilon}^{*}$.

IX-G Convergence of Alternating Projections

In the following, for completeness, we present further results on the convergence of the proposed alternating projections method for the discussed sets in the manuscript. Here, we fix the set $\mathcal{S}$ as the set of nonnegative matrices, $\mathcal{S}_{\rm NN}$. The selection of this feasible set is because this set encompasses the problem of finding stochastic matrices with fixed eigenvalues.

Similar results as the ones shown for the case of set $\mathcal{M}_{\epsilon}^{m}$ can be observed. The error with respect to the projection is always monotonic non-increasing, while the iterate error is not necessarily guaranteed to be monotonic. Despite this, in both cases (sets $\mathcal{M}$ and $\mathcal{M}_{\epsilon}$) the method converges to the desired accuracy of the iterates error.

For these experiments, the employed graph is a $3$-regular graph with $30$ nodes as shown in Fig. 8.

IX-H Response comparison for partially observed network

Here, a comparison of the predicted output associated with the estimated network structures when compared to the true output of the system for an arbitrary excitation. As seen in Fig. 9, an approximate $95\%$ of NRMSE is obtained by using the estimated network structures. This s one of the advantages of using the proposed framework. In instances where the problem can be represented by a state-space formulation, finding a graph that satisfies the spectral characteristic of the identified system, within the set of feasible matrix network representations, allows for consistent network topology identification.