Exact Topology Learning in a Network of Cyclostationary Processes

TL;DR

This paper introduces an exact topology learning algorithm for networks of cyclostationary processes, leveraging a lifting technique and Wiener filter properties, with demonstrated accuracy on a resistor-capacitor network.

Contribution

It is the first method to guarantee exact network topology recovery for cyclostationary processes without structural assumptions.

Findings

Accurate reconstruction demonstrated on resistor-capacitor network.

Performance improves with increasing sample size.

First to guarantee exact recovery in this context.

Abstract

Learning the structure of a network from time series data, in particular cyclostationary data, is of significant interest in many disciplines such as power grids, biology and finance. In this article, an algorithm is presented for reconstruction of the topology of a network of cyclostationary processes. To the best of our knowledge, this is the first work to guarantee exact recovery without any assumptions on the underlying structure. The method is based on a lifting technique by which cyclostationary processes are mapped to vector wide sense stationary processes and further on semi-definite properties of matrix Wiener filters for the said processes.We demonstrate the performance of the proposed algorithm on a Resistor-Capacitor network and present the accuracy of reconstruction for varying sample sizes.

Equations

$$\sum _{n=1}^l{a}_{n,i}\frac{d^n{x}_i}{d{t}^n}=\sum _{j=1,j\ne i}^m{b}_{ij}(x_j(t)-x_i(t))+p_i(t),$$

$${\mathsf{S}}_{\mathsf{i}}(z){\mathsf{x}}_i(z)$$

$${\mathsf{x}}_i(z)$$

$$x_i(k)$$

$$X_j(k)$$

$${\mathsf{X}}_j(z)$$

$${\mathsf{H}}_{ji}(z)$$

$$E_i(k)=$$

$${\mathsf{E}}_i(z)=$$

$${\mathsf{X}}_i(z)=$$

$$D(z)=$$

$${\Phi }_{{\mathsf{P}}_i}(z)=$$

$$\Rightarrow {\Phi }_{{\mathsf{E}}_i}^{-1}(z)=$$

$$X(k)=$$

$$\mathsf{X}(z)=$$

$$\displaystyle X(k)=\begin{bmatrix}X_{1}(k)\\ \vdots\\ X_{m}(k)\end{bmatrix};\$$

$$\displaystyle\mathsf{X}(z)=\mathcal{Z}(X(k));\$$

$${\mathbf{B}}_j^{\prime }{\Phi }_X^{-1}{\mathbf{B}}_i={\mathbf{B}}_j^{\prime }(\mathbb{I}-\mathbb{H}{)}^{\ast }{\Phi }_E^{-1}(\mathbb{I}-\mathbb{H}){\mathbf{B}}_i$$

$$=-{\Phi }_{E_j}^{-1}{\mathsf{H}}_{ji}-{{\mathsf{H}}_{ij}}^{\ast }{\Phi }_{E_i}^{-1}+\sum _{k=1}^m({\mathsf{H}}_{kj}{)}^{\ast }{\Phi }_{E_k}^{-1}({\mathsf{H}}_{ki}).$$

$${\mathbf{B}}_j^{\prime }{\Phi }_X^{-1}(\omega ){\mathbf{B}}_i=-{\Phi }_{E_j}^{-1}(\omega ){\mathsf{H}}_{ji}(\omega )-({\mathsf{H}}_{ij}(\omega ){)}^{\ast }{\Phi }_{E_i}^{-1}(\omega )+$$

$$\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}({\mathsf{H}}_{kj}(\omega ){)}^{\ast }{\Phi }_{E_k}^{-1}(z){\mathsf{H}}_{ki}(\omega )$$

$$=-0-0+\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}({\mathsf{H}}_{kj}(\omega ){)}^{\ast }{\Phi }_{E_k}^{-1}(\omega ){\mathsf{H}}_{ki}(\omega )$$

$$=\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}[D(\frac{\omega }{T}){\mathsf{h}}_{kj}(\frac{\omega }{T}){]}^{\ast }{\Phi }_{E_k}^{-1}(\omega )D(\frac{\omega }{T}){\mathsf{h}}_{ki}(\frac{\omega }{T}).$$

$$\text{ Using }(\text{ref}phi{E}_{p}hiP)$$

$$=\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}{D}^{\ast }(\frac{\omega }{T})\frac{b_{kj}}{{\mathsf{S}}_{\mathsf{k}}^{\ast }(\frac{\omega }{T})}{\Phi }_{P_k}^{-1}|{\mathsf{S}}_{\mathsf{k}}(\frac{\omega }{T}){|}^{2}D(\frac{\omega }{T})\frac{b_{ki}}{{\mathsf{S}}_{\mathsf{k}}(\frac{\omega }{T})}$$

$$=\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}(\frac{b_{kj}{b}_{ki}}{|{\mathsf{S}}_{\mathsf{k}}(\frac{\omega }{T}){|}^2})|{\mathsf{S}}_{\mathsf{k}}(\frac{\omega }{T}){|}^2[D^{\ast }(\frac{\omega }{T}){\Phi }_{P_k}^{-1}(\omega )D(\frac{\omega }{T})]$$

$$=\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}(b_{kj}{b}_{ki})[D^{\ast }(\frac{\omega }{T}){\Phi }_{P_k}^{-1}(\omega )D(\frac{\omega }{T})]$$

$${\mathbf{B}}_j^{\prime }{\Phi }_X^{-1}(\omega ){\mathbf{B}}_i$$

$$=\sum _{k\in {\mathcal{C}}_G(j)\cap {\mathcal{C}}_G(i)}$$

$$\forall \omega \in [0,2\pi ).$$

$$x_1(k)$$

$$x_2(k)$$

$$x_3(k)$$

$${\mathsf{x}}_i(z)$$

$$x_i(k)$$

$$x_i(kT)$$

$$x_i(kT+p)$$

$$\text{Put }t=p-b$$

$$x_i(kT+p)$$

$$e_i(kT+p)$$

$$X_i(k)=\left[\begin{array}{c}{x}_i(kT)\\ x_i(kT+1)\\ ⋮\\ x_i(kT+T-1)\end{array}\right]$$

Full text

I INTRODUCTION

Networks are a framework used extensively for analysis of the behavior of complex systems like the brain [1], climate [2], gene regulation dynamics [3], disease spread [4], power grid [5] and others. Network representations help identify essential influence pathways in a complex system, thereby enhance the understanding of the dynamics of complex systems. A problem of interest to multiple communities like control theory, machine learning and signal processing is the inference of the influence pathways among the entities of interest from observation of the entities, which is sometimes referred as structure learning or topology learning [6, 7, 8]. For example: given a time series collection of stock prices over a time horizon, it is of interest to infer the influence paths among the collection of stocks from the observed time series of stock prices.

Fundamentally, there exist two different approaches for inference of the influence pathways among the entities from observations. The first is active learning, where a particular entity is perturbed by an external agent and its influence on the rest of the entities is examined [9]. However, the disadvantage of such an approach is that it is not always possible to actively excite a specific entity to glean the network structure. For example: it is not always possible to turn off or change generator set points to infer the structure of a power distribution network. The second approach to structure learning involves passive or non invasive approaches, where, the system is not perturbed actively and topology is inferred solely from the observations of the variables of interest [10],[11], [12],[14]. In this article we focus on the passive approach to topology inference.

Inference of the network structure from observations assuming that the entities are a collection of random variables is studied in [10, 14]. However, such approaches are ineffective when lagged correlations are significant and the past of one entity influences the present of another entity. In such situations, the observed entities are modeled as stochastic processes or time series. In the time series setting, topology inference under the assumption of wide sense stationary time series using multivariate Wiener filtering is described in [15] and using power spectral density in [16]. These approaches are not applicable to the case of non stationary time series. In this article we utilize multivariate Wiener filtering for topology inference with consistency guarantees for non stationary time series, which are cyclostationary. Cyclostationary processes are characterized by a periodic mean and correlation function. Many phenomenon in nature exhibit periodic behaviour [17]; examples include weather patterns [18][2], regulatory processes in human body [19], [20], [21], trends in stock market [22], [23], motion of mechanical systems [24], [25], communication systems [26], [27], and planet movements [28]. In [29], the authors use multivariate Wiener filtering for topology inference in a collection of cyclostationary time series. However, the approach presented in [29] introduces possibly many spurious links and does not provide exact inference of the underlying network topology from the observations. In this article, we extend the work in [29] and present an algorithm with guarantees for exact reconstruction of a network structure in a collection of cyclostationary processes. Our approach does not require any structural assumptions on the underlying graph structure nor any knowledge of the system parameters. Our results are based on the properties of the Wiener filters of collection of wide sense stationary time series for inference of exact topology in a network of wide sense stationary processes [30, 31]. We validate the algorithm presented through implementation on time series of nodal states generated from an interconnected Resistor-Capacitor (RC) network with input cyclostationary process. We compare the performance of our algorithm against the algorithms presented in [29] and show the superiority of the approach presented here over prior work. To the authors’ best knowledge, this is the first work on exact topology reconstruction of cyclostationary processes with provable guarantees.

The rest of the paper is organized as follows. Section II describes the learning problem of a general linear dynamical system with cyclostationary inputs. In section III, the approach for tackling the topology learning problem from cyclostationary time series data is presented. Analytic results as well as an implementable algorithm is presented. The performance of the algorithm is demonstrated on a simulated RC network in Section IV. Conclusions are presented in Section V.

Notation:

The symbol $:=$ denotes a definition

$\boldsymbol{0}$ : zero matrix of appropriate dimension

$A^{*}$ : the conjugate transpose of matrix $A$

$A^{\prime}$ : transpose of a matrix $A$

$\|A\|_{\infty}$ : largest row sum of matrix $A$

$A\succeq 0$ : $A$ is a positive semi-definite matrix

$A\preceq 0$ : $A$ is a negative semi-definite matrix

II Linear Dynamical System with cyclostationary inputs

In this section, we briefly define and describe a networked cyclostationary process.

II-A Introduction to Cyclostationary Process

A random process $x(t)\in L^{2}(\Omega,\mathcal{F},\mathcal{P})$, with its mean function $m(t):=\mathbb{E}[x(t)]$ and correlation function $R_{x}(s,t):=\mathbb{E}[x(s)x(t)]$ is wide sense cyclostationary (WSCS) or periodically correlated with period $T>0$ if for every $s,t\in\mathbb{Z}$ there is no value smaller than $T$ such that $m(t)=\mathbb{E}[x(t)]=\mathbb{E}[x(t+T)]=m(t+T)$ and $R_{x}(s,t)=\mathbb{E}[x(s)x(t)]=\mathbb{E}[x(s+T)x(t+T)]=R_{x}(s+T,t+T)$. The random processes $x(t)$ and $e(t)$ are said to be jointly wide sense cyclostationary (JWSCS) if $x(t)$ and $e(t)$ are cyclostationary with period $T$ and the cross correlation function $R_{x,e}(s,t)$ is periodic with period $T$ that is $R_{x,e}(s,t)$ $=R_{x,e}(s+T,t+T)$. An example of a cyclostationary process is a periodic signal superimposed with a wide sense stationary (WSS) signal. A WSS signal by itself is also a cyclostationary signal with period $T=1$.

II-B Network of Cyclostationary Processes

Consider a system with $m$ sub-systems indexed by $i$ with the dynamics of each sub-system $x_{i}$ given by

[TABLE]

$i\in\{1,2,..,m\}$, where $x_{i}(t)$ is the time series data of the sub-system $i$; $a_{m,i}$ and $b_{ij}$ are system parameters and $p_{i}(t)$ is a forcing function acting on sub-system $i$, which is a zero mean WSCS of period $T$, uncorrelated with $\{p_{j}(t)\}_{i\neq j}$. Since a linear transformation of $\{p_{j}(k)\}_{j=1}^{m}$ results in $\{x_{j}(k)\}_{j=1}^{m}$, the collection of $(\{x_{j}(k)\}_{j=1}^{m},\{p_{j}(k)\}_{j=1}^{m})$ are JWSCS [29]. All the processes and the system parameters of (1) are real valued and $b_{ij}$ are non-negative. Equation (1) is known as Linear Dynamical Model (LDM).

Converting (1) from time domain to z-domain by applying bilinear transform (Tustin’s method [32]) on (1), we obtain

[TABLE]

where, $\mathsf{h}_{ii}(z)=0$, $\mathsf{S_{i}}(z)=\sum_{n=1}^{l}a_{n,i}(\frac{2(1-z^{-1})}{\Delta t(1+z^{-1})})^{n}+\sum_{j=1,j\neq i}^{m}b_{ij}$ is a scalar function, $\mathsf{h}_{ij}(z)=\frac{b_{ij}}{\mathsf{S_{i}(z)}}$, $\mathsf{e}_{i}(z)=$ $\frac{\mathsf{p}_{i}(z)}{\mathsf{S_{i}(z)}}$. Note that $\mathsf{h}_{ij}(z)$ is the Z-transform of $h_{ij}(n)$, that is $\mathsf{h}_{ij}(z)=\mathcal{Z}[h_{ij}(n)]$. Here, $\mathsf{x}_{i}(z)$, $\mathsf{p}_{i}(z)$, $\mathsf{e}_{i}(z)$ are the z-transforms of $x_{i}(k),p_{i}(k)\text{ and }e_{i}(k)$ respectively. $e_{j}(k)$ is the input noise which is zero mean WSCS and uncorrelated with $\{e_{i}\}_{i=1,i\neq j}^{m}$.

The Linear Dynamic Graph (LDG) given by $\mathcal{G}(V,E)$ associated with the LDM described by (1) is obtained by setting each vertex (node) $i\in V$ to represent every subsystem $i$ and by placing a direct edge in $E$ from vertex $i$ to vertex $j$ if $b_{ji}\not=0$. The set of children of node $x_{j}$ is defined as $\mathcal{C}_{G}(x_{j}):=\{x_{i}|(x_{j},x_{i})\in E\}$, its set of parents as $\mathcal{P}_{G}(x_{j}):=\{x_{i}|(x_{i},x_{j})\in E\}$ and its set of kins as $\mathcal{K}_{G}(x_{j}):=\{x_{i}|x_{i}\neq x_{j}\text{ and }x_{i}\in\mathcal{C}_{G}(x_{j})~{}\cup\mathcal{P}_{G}(x_{j})~{}\cup~{}\mathcal{P}_{G}(\mathcal{C}_{G}(x_{j}))\}.$ The kin topology/moral graph of the directed graph $\mathcal{G}(V,E)$ is defined as an undirected graph $\mathcal{G}_{M}(V,\tilde{E})$ with a vertex set $V$ and an edge set $\tilde{E}=\{(x_{i},x_{j})|x_{i}\in\mathcal{K}_{G}(x_{j})\forall j\}$. The topology of the directed graph $\mathcal{G}(V,E)$ is defined as an undirected graph $\mathcal{G}_{T}(V,\tilde{E})$ with a vertex set $V$ and an edge set $\tilde{E}=\{(x_{i},x_{j})|(x_{i},x_{j})\text{ or }(x_{j},x_{i})\in E\}$ and is denoted by $top(\mathcal{G})=\mathcal{G}_{T}$. Example of topology and kin topology are shown in Fig. 1(c) and (d) respectively. In the topology $\mathcal{G}_{T}(V,\tilde{E})$ of the directed graph $\mathcal{G}$, a neighbor is defined as $\mathcal{N}_{\mathcal{G}_{T}}(x_{j}):=\{x_{i}|(x_{i},x_{j})\in\tilde{E}\}$ and a two hop neighbor is defined as $\mathcal{N}_{\mathcal{G}_{T}}(x_{j},2):=\{x_{i}\in V|(x_{j},x_{k}),(x_{i},x_{k})\in\tilde{E}\text{ for some }x_{k}\in V\}$

II-C Lifting to vector wide sense stationary processes

The vector random process $X_{i}(k)$ is said to be (weakly) stationary if $E([X_{i}(k)]_{j})=m_{i}(j)$, where $[X_{i}(k)]_{j}$ is the $j^{th}$ element of $X_{i}(k)$ and $R_{X_{i}}^{jk}(s,t)=E([X_{i}(s)]_{j}[X_{i}(t)]_{k})=E(x_{i}(s+j-1)x_{i}(t+k-1))=E(x_{i}(s-t+j)x_{i}(k))=R_{X_{i}}^{jk}(s-t)$ for all $s,t\in\mathbb{Z},\text{ and }j,k\in\{1,2,\cdots,T\}$

Any cyclostationary process with a period $T$ can be represented as multivariate, vector, wide sense stationary process through the lifting process [33] described next. Processes $x_{i}(k)$ and $e_{i}(k)$ are lifted to $T$-variate WSS processes $X_{i}(k)=[x_{i}(kT)\ ...\ x_{i}(kT+T-1)]^{\prime}$ and $E_{i}(k)=[e_{i}(kT)\ ...\ e_{i}(kT+T-1)]^{\prime}$ respectively and their z-transform’s are denoted as $\mathsf{X}_{i}(z)$ and $\mathsf{E}_{i}(z)$ respectively. The vector process $E_{i}(k)$ is uncorrelated with $\{E_{j}(k)\}_{j\neq i}$.

Lemma II.1

Given the dynamics of a network of cyclostationary processes $\{x_{i}(k)\}_{i=1}^{m}$ in (II-B), the dynamics of the lifted vector WSS processes $\{X_{i}(k)\}_{i=1}^{m}$ is given by $\mathsf{X}_{j}(z)=\sum_{i=1}^{m}\mathsf{H}_{ji}(z)\mathsf{X}_{i}(z)+\mathsf{E}_{j}(z)$, where $\mathsf{H}_{ji}(z)=D(z^{\frac{1}{T}})\mathsf{h}_{ji}(z^{\frac{1}{T}})$

Proof:

See the appendix ∎

The dynamics of $\{X_{j}\}_{j=1}^{m}$ is given by

[TABLE]

Here, $\mathsf{H}_{ji}(z)$ is a $T\times T$ transfer matrix, $\mathsf{E}_{i}(z)=\mathcal{Z}(E_{i}(k))$ and $\mathsf{X}_{i}(z)=\mathcal{Z}(X_{i}(k))$. The collection $\{X_{i}(k)\}_{i=1}^{m}$ are jointly vector WSS processes. The LDG for the vector processes $\{X_{i}(k)\}_{i=1}^{m}$ is defined as follows: a directed edge from $X_{i}$ to $X_{j}$ exists if and only if $\mathsf{h}_{ji}\neq 0$.

Lemma II.2

A directed edge from a cyclostationary process $x_{i}$ to $x_{j}$ in the LDG of (1) exists if and only if there exists a directed edge from $X_{i}$ to $X_{j}$

Proof:

Given that a directed edge exists from cyclostationary process $x_{i}$ to $x_{j}$ then $\mathsf{h}_{ji}(z)\neq 0$ in (II-B). It follows from (4) that $\mathsf{H}_{ji}(z)\neq\boldsymbol{0}$ which implies that there is a directed edge from vector stationary process $X_{i}$ to $X_{j}$. To prove the converse, if suppose there is a directed edge from $X_{i}$ to $X_{j}$ then $\mathsf{H}_{ji}(z)\neq\boldsymbol{0}$ in (4). This implies that $\mathsf{h}_{ji}(z)\neq 0$ from (4) and there is a directed edge from $x_{i}$ to $x_{j}$. ∎

Remark 1

The above lemma concludes that kin relationship’s in the LDG of $\{x_{j}\}_{j=1}^{m}$ identical to the kin relationship’s in the LDG of $\{X_{j}\}_{j=1}^{m}$. This enables us to reconstruct the topology of LDG for the cyclostationary processes $\{x_{j}\}_{j=1}^{m}$ by reconstructing the topology for their equivalent vector stationary processes $\{X_{j}\}_{j=1}^{m}$.

The Power Spectral Density (PSD) of the vector WSS process $X_{i}(k)$ is related to that of scalar WSCS process $x_{i}(k)$. The PSD matrix ${\Phi}_{P_{i}}(z),{\Phi}_{E_{i}}(z)$ are block matrices. Note that $\mathsf{P}_{i}(z)=\mathsf{S_{i}}(z^{\frac{1}{T}})\mathsf{E}_{i}(z)$. Therefore,

[TABLE]

The rest of the article focuses on the reconstruction of topology of a dynamical network of $T\times 1$ vector stationary processes that are jointly wide sense stationary (in contrast to $x_{i}$ which are jointly wide sense cyclostationary).

III Learning the structure from time series data

Here we assume that $x_{i}(k)$ is available as a measured time series for all subsystems $i$. Given the time series data of the dynamical system with noise modeled as cyclostationary, the main aim of this article is to reconstruct the topology of the LDG of the system. The time series $x_{i}(k)$ is lifted to vector WSS process $X_{i}(k)$. The output processes $\{X_{j}(k)\}_{j=1}^{m}$ of (4) are defined in a compact way as:

[TABLE]

The transfer block matrix $\mathbb{H}(z)$ is $mT\times mT$ matrix with size of each block $T\times T$ and its diagonal entries are $\boldsymbol{0_{T\times T}}$. The matrix $\mathbb{H}(z)$ in (7) characterizes the interconnection between the processes $X_{i}(k)$ and $X_{j}(k)$ with the $(i,j)^{th}$ block of $\mathbb{H}(z)$ being $\mathsf{H}_{ij}(z)$. The relation (7) defines a map from a vector processes $E(k)$ to a vector processes $X(k)$. The LDM described by (7) is well-posed if the operator $(\mathbb{I}-\mathbb{H}(z))$ is invertible almost surely. Therefore for any vector processes ${E}(k)$ there exists a vector ${X}(k)$ of processes and the LDM described by (7) is said to be topologically detectable if $\Phi_{E_{j}}\succ 0$ for any $\omega\in[-\pi,\pi]$ and for any $j=1,...,m$.

The proposed algorithm in this article for exact reconstruction of the true topology is based on inverse PSD of the time series described above.

III-A Reconstruction of Moral Graph

In this section, we present an algorithm for reconstruction of the network structure of LDG $\mathcal{G}$ of (4), which is well-posed and topologically detectable.

From (7) $X=(\mathbb{I}-\mathbb{H}(z))^{-1}E$ which implies that $\Phi_{X}^{-1}=(\mathbb{I}-\mathbb{H})^{*}\Phi_{E}^{-1}(\mathbb{I}-\mathbb{H})$. We now consider the $(i,j)$ block matrix of $\Phi_{X}^{-1}$, for $i\neq j$. Using the block diagonal structure of $\Phi_{E}$, we have

[TABLE]

where $\boldsymbol{B}_{j}=[\boldsymbol{0}\ \boldsymbol{0}\ ..\ \boldsymbol{0}\ \boldsymbol{I}_{T\times T}\ \boldsymbol{0}\ ..\ \boldsymbol{0}]^{\prime}$ is a matrix $\in\mathbb{R}^{mT\times T}$ that has an identity matrix $\boldsymbol{I}_{T\times T}$ as the $j^{th}$ block and other blocks as $\boldsymbol{0}$.

If $(i,j)$ are not kins then $j$ cannot be the parent of $i$ ($\mathsf{H}_{ij}=\boldsymbol{0}$) or child of $i$ ($\mathsf{H}_{ji}=\boldsymbol{0}$) or spouse of $i$ (which implies that there does not exist any $k\in\{1,...,m\}\text{ with }\mathsf{H}_{kj}\neq\boldsymbol{0}\ \text{and}\ \mathsf{H}_{ki}\neq\boldsymbol{0}$).

This gives a sufficient criterion to identify kins which include parents, children, spouses and the moral graph is obtained (see Fig. 1(d)). The limitation with the moral graph is that it may include substantial spurious links between two-hop neighbors as shown in Fig. 1(d). Thus, it is important to devise a method for the detection and removal of spurious links from the moral graph and to deduce the topology of the LDG; which forms the focus of the next section.

III-B Pruning the Spurious Strict Two Hop Neighbor Links

The following theorem will provide a sufficient condition for elimination of all strict two-hop neighbors in the moral graph $\mathcal{G}_{M}$ which results in true topology ${\mathcal{G}_{T}}$

Theorem III.1

Consider a LDM $(\mathbb{H}(z),{E})$ which is well-posed and topologically detectable, with its associated LDG $\mathcal{G}$ and topology $\mathcal{G_{T}}$. Let the output of the LDM be given by $X(k)$ according to (7). Let the vertices $i$ and $j$ in $\mathcal{G_{T}}$ be such that, $i\in N_{\mathcal{G_{T}}}(j,2)$ but $i\notin N_{\mathcal{G_{T}}}(j)$, that is, $i,j$ are strict two-hop neighbors. Then, $\boldsymbol{\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}\succeq 0$, for all $\omega\in[0,2\pi)$.

Proof:

Given that $i\in N_{\mathcal{G_{T}}}(j,2)$ and $i\notin N_{\mathcal{G_{T}}}(j)$. Using (III-A)

[TABLE]

It follows from (6) that ${{\Phi}_{P_{k}}}(\omega)$ is also positive definite, thus ${{\Phi}_{P_{k}}^{-1}}(\omega)\succ 0$ and $[D^{*}(\frac{\omega}{T}){\Phi}_{P_{k}}^{-1}(\omega)D(\frac{\omega}{T})]\succeq 0$.

For a common child $k\in\mathcal{C}_{G}(j)\cap\mathcal{C}_{G}(i)$, the values of $b_{kj},b_{ki}$ strictly positive. Thus,

[TABLE]

∎

Remark 2

The consequence of Theorem III.1 is that if $i$ is a strict two hop neighbor of $j$, then all the eigenvalues of the matrix ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}(\omega)}{\boldsymbol{B}_{i}}$ are non-negative. This theorem does not guarantee on its converse, but such cases are pathological as it will be evident in next theorem.

Theorem III.2

Given a well-posed and topologically detectable LDM $(\mathbb{H}(z),{E})$ with associated graph $\mathcal{G}$ and its topology $\mathcal{G_{T}}$, the following holds:

Suppose nodes $i$ and $j$ are such that $i\in N_{\mathcal{G_{T}}}(j)$, $i\not\in N_{\mathcal{G_{T}}}(j,2)$, ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}\succeq 0$ for all $\omega\in[0,2\pi)$. Then

$[-{\frac{b_{ji}}{{\mathsf{S}_{j}^{*}}(\frac{\omega}{T})}{\Phi}_{P_{j}}^{-1}(\omega)}D(\frac{\omega}{T})-\frac{b_{ij}}{\mathsf{S}_{i}(\frac{\omega}{T})}D^{*}(\frac{\omega}{T}){{\Phi}_{P_{i}}^{-1}(\omega)}]\succeq 0$** 2. 2.

Suppose $i$ and $j$ are such that $i\in N_{\mathcal{G_{T}}}(j)$ and $i\in N_{\mathcal{G_{T}}}(j,2)$ (both one hop and two hop neighbor), with ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}\succeq 0$ for all $\omega\in[0,2\pi)$. Then $[-{\frac{b_{ji}}{{\mathsf{S}_{j}^{*}}(\frac{\omega}{T})}{\Phi}_{P_{j}}^{-1}(\omega)}D(\frac{\omega}{T})-\frac{b_{ij}}{\mathsf{S}_{i}(\frac{\omega}{T})}D^{*}(\frac{\omega}{T}){{\Phi}_{P_{i}}^{-1}(\omega)}+{\sum_{k\in\mathcal{C}_{G}(j)\cap\mathcal{C}_{G}(i)}}{(b_{kj}b_{ki})}[D^{*}(\frac{\omega}{T}){\Phi}_{P_{k}}^{-1}(\omega)D(\frac{\omega}{T})]]\succeq 0$

Remark 3

The consequence of theorem III.2 is that for nodes $i$ and $j$ that are neighbors but not two hop neighbors, or, nodes $i$ and $j$ that are neighbors and two hop neighbors, the inequalities in theorem III.2 holds for all $\omega\in[0,2\pi)$ when the system parameters satisfy a restrictive and specific set of conditions. In other words, aside for pathological cases, the converse of Theorem III.2 holds. So if the ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}\succeq 0$ then $i\in N_{\mathcal{G_{T}}}(j,2)$. We use this condition to prune out spurious two hop neighbor edges using the algorithm discussed next.

III-C Topology Learning Algorithm

The steps involved to unravel the structure of a LDG of cyclostationary process are summarized in Algorithm table $1$.

IV Results

IV-A Data generation

In this section, a generative graph is considered as shown in Fig.(3) along with its moral graph and topology. Each node $i$ is driven by exogenous input $e_{i}(k)$. The processes $\{e_{2}(k),e_{3}(k)\}$ are wide sense stationary and $e_{1}(k)$ is cyclostationary with $T=2$. There are two directed edges in the generative graph where nodes $\{1,3\}$ are the parents of node $2$. The dynamics of the directed edge from node $1$ to node $2$ is given by $\mathsf{h}_{21}(z)=0.1\{1+0.9z^{-1}+0.5z^{-2}+0.3z^{-3}\}$. Similarly, $\mathsf{h}_{23}(z)=0.2\{1+0.9z^{-1}+0.5z^{-2}+0.3z^{-3}\}$. The data is generated with the following generative model:

[TABLE]

The sample size of each time series $\{x_{1}(k),x_{2}(k),x_{3}(k)\}$ is $3\times 10^{5}$.

IV-B Topology Reconstruction

With the data generated in (10), we can apply the presented algorithm to reconstruct the Topology of the Generative graph. The steps involved in reconstructing the Topology are as follows:

Identify the time period T: Performing the periodogram analysis of the time series $\{x_{1}(k),x_{2}(k),x_{3}(k)\}$, we found the value of $T$ as $2$. For example, the periodogram of $x_{1}(k)$ is shown in Fig.(2). 2. 2.

Lifting the time series: Each time series $x_{i}(k)$ is lifted to a vector time series $X_{i}(k)$ of length $T$. That is, $X_{i}(k)=[x_{i}(2k),x_{i}(2k+1)]^{T}$ for $i=\{1,2,3\}$. 3. 3.

Compute ${\hat{\Phi}_{X}^{-1}(\omega)}$: Using the Welch method [13], an estimate of ${{\Phi}_{X}^{-1}(\omega)}$ is computed. This can be done in MATLAB using cpsd(). For example, ${\boldsymbol{B}_{1}^{\prime}}{\hat{\Phi}_{X}^{-1}}{\boldsymbol{B}_{3}}(f)$ is $2\times 2$ matrix and is shown in Fig.4. 4. 4.

Moral graph reconstruction: We place an undirected edge between two distinct nodes $j$ and $i$ if ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}\neq 0$. Repeating this process for all pairs of nodes, the graph that is constructed matches with the moral graph shown in Fig.3(b). 5. 5.

Topology reconstruction: For every edge $(j,i)$ in the reconstructed moral graph, we remove this edge if ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}\succeq 0$. The eigenvalues of ${\boldsymbol{B}_{j}^{\prime}}{{\Phi}_{X}^{-1}}{\boldsymbol{B}_{i}}$ are shown in Fig. 5. The graph obtained after iterating over all the edges in the moral graph is identical to the topology of the generative graph in Fig.3(c).

This completes the topology reconstruction of the times series $\{x_{1}(k),x_{2}(k),x_{3}(k)\}$. Since we generated data, we calculated the exact expression of ${{\Phi}_{X}^{-1}(\omega)}$ and compared with the estimates in Fig.4 and Fig.5.

V Conclusions

In summary, the topology of the Linear Dynamic Graph of a Linear Dynamic Model with cyclostationary inputs is constructed from the time series data with provable guarantees. It is a data driven approach useful for exact identification of network structure with applications to power grids, thermal networks, network of rotating mechanical systems amongst many others. Our approach doesn’t impose any structural restrictions on the network topology and doesn’t use any knowledge of the system parameters.

In our future work, we will apply this algorithm on experimental data with a wide test cases for the parameters. Regularization methods will be used to improve the performance with less samples and to promote the sparseness. This algorithm can also be used to learn the topology of a network of WSS processes more efficiently than the current state of the art and will be the focus of our future work.

VI Appendix

Lemma 2.1

Proof:

The dynamics of a network of cyclostationary processes are given by

[TABLE]

where $k$ is the time index and $p\in\{0,1,\dots,T-1\}$. Substitute $n=aT+b$, where $b$ takes the values $\{p,p-1,p-2,\dots,p-T+1\}$

[TABLE]

Lets define $H_{ij,pt}[a]=h_{ij}(aT+p-t)$. The value $p$ is the row index starting from [math] to $T-1$ and $t$ is the column index starting from [math] to $T-1$ of the filter matrix $H_{ij}$. Note that for $p=t$ (diagonal entries) $H_{ij,pp}[n]=h_{ij}(nT)$ for $p\in\{0,1,\dots,T-1\}$.

[TABLE]

Lets lift the scalar process $x_{i}(k)$ to a vector process $X_{i}(k)$ by varying $p$ from [math] to $T-1$.

[TABLE]

Iterate the $p$ in equation (3) from [math] to $T-1$ to get the following relation

[TABLE]

[TABLE]

where $H_{ij,pt}*x_{j}(kT+t)=\sum_{a=-\infty}^{\infty}H_{ij,pt}[a]x_{j}(kT+t-aT)$ for $t\in\{0,1,\dots,T-1\}$ with $H_{ij,pt}[n]=h_{ij}(nT+p-t)$. This implies $\mathsf{H}_{ji}(z)=\mathcal{Z}[H_{ji}(k)]=D(z^{\frac{1}{T}})\mathsf{h}_{ji}(z^{\frac{1}{T}})$, where

$D(z)=\begin{bmatrix}z^{0}&z^{1}\dots&z^{T-1}\\ \vdots&\ddots&\vdots\\ z^{-(T-1)}&\dots&z^{0}\\ \end{bmatrix}$. ∎