Controllability of Bandlimited Graph Processes Over Random Time Varying Graphs

TL;DR

This paper studies how to control network states over randomly changing graphs by leveraging bandlimited graph signals, proposing methods to ensure expected controllability and minimize deviations in stochastic environments.

Contribution

It introduces the concept of controllability in the mean for bandlimited graph processes on random time-varying graphs and develops control strategies accounting for topological randomness.

Findings

Control strategies effectively drive expected network states towards targets.

Graph filters can remove out-of-band frequency content caused by control.

Proposed methods perform well on synthetic and social network models.

Abstract

Controllability of complex networks arises in many technological problems involving social, financial, road, communication, and smart grid networks. In many practical situations, the underlying topology might change randomly with time, due to link failures such as changing friendships, road blocks or sensor malfunctions. Thus, it leads to poorly controlled dynamics if randomness is not properly accounted for. We consider the problem of controlling the network state when the topology varies randomly with time. Our problem concerns target states that are bandlimited over the graph; these are states that have nonzero frequency content only on a specific graph frequency band. We thus leverage graph signal processing and exploit the bandlimited model to drive the network state from a fixed set of control nodes. When controlling the state from a few nodes, we observe that spurious, out-of-band frequency content is created. Therefore, we focus on controlling the network state over the desired frequency band, and then use a graph filter to get rid of the unwanted frequency content. To account for the topological randomness, we develop the concept of controllability in the mean, which consists of driving the expected network state towards the target state. A detailed mean squared error analysis is performed to quantify the statistical deviation between the final controlled state on a particular graph realization and the actual target state. Finally, we propose different control strategies and evaluate their effectiveness on synthetic network models and social networks.

Equations

$$\mathbf{x}={\mathbf{V}}_K{\tilde{\mathbf{x}}}_K$$

$$\frac{d\mathbf{x}(s)}{ds}=-\mathbf{L}\mathbf{x}(s),\mathbf{x}(0)={\mathbf{x}}_0.$$

$${\mathbf{x}}_t=\mathbf{A}{\mathbf{x}}_{t-1},\mathbf{A}=\mathbf{I}-ϵ\mathbf{L},t\in \mathbb{N}$$

$${\mathbf{x}}_t=\mathbf{A}{\mathbf{x}}_{t-1},\mathbf{A}=\mathbf{W},t\in \mathbb{N}.$$

$${\mathbf{x}}_t=\mathbf{A}{\mathbf{x}}_{t-1}+\mathbf{B}{\mathbf{u}}_{t-1}$$

$$\mathbf{\Omega }=[\mathbf{B},\mathbf{A}\mathbf{B},\dots ,{\mathbf{A}}^{T-1}\mathbf{B}]$$

$${\mathcal{C}}_{M,N}=\left\{\mathbf{C}\in \{0,1{\}}^{M\times N}:\mathbf{C}{1}_N=1_M,{\mathbf{C}}^{\mathsf{T}}{1}_M\le 1_N\right\}$$

$$\begin{array}{rl}{\tilde{\mathbf{x}}}_t& ={\mathbf{V}}^{\mathsf{H}}\mathbf{A}\mathbf{V}{\tilde{\mathbf{x}}}_{t-1}+{\mathbf{V}}^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}\\ & \triangleq \tilde{\mathbf{A}}{\tilde{\mathbf{x}}}_{t-1}+{\mathbf{V}}^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}\end{array}$$

$$\left[\begin{array}{c}{\tilde{\mathbf{x}}}_{t,K}\\ {\tilde{\mathbf{x}}}_{t,N-K}\end{array}\right]=\left[\begin{array}{c}\text{diag}({\mathbf{a}}_K){\tilde{\mathbf{x}}}_{t-1,K}\\ \text{diag}({\mathbf{a}}_{N-K}){\tilde{\mathbf{x}}}_{t-1,N-K}\end{array}\right]+\left[\begin{array}{c}{\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}\\ {\mathbf{V}}_{N-K}^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}\end{array}\right]$$

$${\tilde{\mathbf{x}}}_{t,K}={\tilde{\mathbf{A}}}_K{\tilde{\mathbf{x}}}_{t-1,K}+{\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}.$$

$$\underset{\mathbf{u}}{\mathrm{minimize}}$$

$${\mathbf{x}}_{t+1}=\mathbf{A}{\mathbf{x}}_t+\mathbf{B}{\mathbf{u}}_t,t=0,\dots ,T-1$$

$${\mathbf{x}}_0=0_N,{\mathbf{x}}_T={\mathbf{x}}^{\ast }$$

$${\mathbf{x}}_t={\mathbf{A}}_{t-1}{\mathbf{x}}_{t-1}+{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}$$

$${\mathbf{\mu }}_t$$

$${\tilde{\mathbf{\mu }}}_{t,K}={\check{\mathbf{A}}}_K{\tilde{\mathbf{\mu }}}_{t-1,K}+{\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}$$

$$\begin{array}{rl}& \mathrm{MSE}(T)=\mathbb{E}\left[{\Vert \mathbf{H}{\mathbf{x}}_T-{\mathbf{x}}^{\ast }\Vert }_2^2\right]\\ & =\alpha -2\sum _{\tau =0}^{T-1}{\mathbf{\beta }}_{\tau }^{\mathsf{T}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{\tau }+\sum _{\tau =0}^{T-1}\sum _{{\tau }^{\prime }=0}^{T-1}{\mathbf{u}}_{\tau }^{\mathsf{T}}\mathbf{C}{\mathbf{\Gamma }}_{\tau ,{\tau }^{\prime }}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{{\tau }^{\prime }}\end{array}$$

$${\tilde{\mathbf{\mu }}}_{T,K}=\sum _{\tau =0}^{T-1}{\check{\mathbf{A}}}_K^{T-\tau -1}{\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{\tau }$$

$$\tilde{\mathbf{\Omega }}\mathbf{u}=\left[{\mathbf{I}}_K,{\check{\mathbf{A}}}_K,\cdots ,{\check{\mathbf{A}}}_K^{T-1}\right]\left({\mathbf{I}}_T\otimes {\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}\right)\mathbf{u}={\tilde{\mathbf{x}}}_K^{\ast }$$

$$\underset{\mathbf{C}\in {\mathcal{C}}_{M,N}^{\ast },\mathbf{u}\in {\mathbb{R}}^{TM}}{\min }$$

$$\tilde{\mathbf{\Omega }}\mathbf{u}={\tilde{\mathbf{x}}}_K^{\ast }$$

$$\tilde{\mathbf{\Omega }}=\left[{\mathbf{I}}_K,{\check{\mathbf{A}}}_K,\cdots ,{\check{\mathbf{A}}}_K^{T-1}\right]\left({\mathbf{I}}_T\otimes {\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}\right)$$

$${\mathbf{u}}^{\ast }={\tilde{\mathbf{\Omega }}}^{\mathsf{H}}{\left[\tilde{\mathbf{\Omega }}{\tilde{\mathbf{\Omega }}}^{\mathsf{H}}\right]}^{-1}{\tilde{\mathbf{x}}}_K^{\ast }.$$

$$\underset{\mathbf{C}\in {\mathcal{C}}_{M,N}^{\ast }}{\min }$$

$${\mathbf{u}}^{\ast }={\tilde{\mathbf{\Omega }}}^{\mathsf{H}}{\left[\tilde{\mathbf{\Omega }}{\tilde{\mathbf{\Omega }}}^{\mathsf{H}}\right]}^{-1}{\tilde{\mathbf{x}}}_K^{\ast },$$

$$\tilde{\mathbf{\Omega }}=\left[{\mathbf{I}}_K,{\check{\mathbf{A}}}_K,\cdots ,{\check{\mathbf{A}}}_K^{T-1}\right]\left({\mathbf{I}}_T\otimes {\mathbf{V}}_K^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}\right).$$

$$\frac{\partial \mathrm{MSE}(T)}{\partial {\mathbf{u}}_t}=-2\mathbf{C}{\mathbf{\beta }}_t+2\sum _{\tau =0}^{T-1}\mathbf{C}{\mathbf{\Gamma }}_{t,\tau }{\mathbf{C}}^T{\mathbf{u}}_{\tau }$$

$$\mathbf{\Gamma }=\left[\begin{array}{cccc}{\mathbf{\Gamma }}_{T-1,T-1}& {\mathbf{\Gamma }}_{T-1,T-2}& \cdots & {\mathbf{\Gamma }}_{T-1,0}\\ {\mathbf{\Gamma }}_{T-2,T-1}& {\mathbf{\Gamma }}_{T-2,T-2}& \cdots & {\mathbf{\Gamma }}_{T-2,0}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{\Gamma }}_{0,T-1}& {\mathbf{\Gamma }}_{0,T-2}& \cdots & {\mathbf{\Gamma }}_{0,0}\end{array}\right]\in {\mathbb{R}}^{NT\times NT}$$

$${\mathbf{\Gamma }}_C\mathbf{u}=\left({\mathbf{I}}_T\otimes \mathbf{C}\right)\mathbf{\Gamma }\left({\mathbf{I}}_T\otimes {\mathbf{C}}^{\mathsf{T}}\right)\mathbf{u}=\left({\mathbf{I}}_T\otimes \mathbf{C}\right)\mathbf{\beta }={\mathbf{\beta }}_C$$

$${\mathbf{u}}_C^{\ast }={\mathbf{\Gamma }}_C^{-1}{\mathbf{\beta }}_C.$$

$$\underset{\mathbf{C}\in {\mathcal{C}}_{M,N}^{\ast }}{\min }$$

$${\mathbf{u}}_C^{\ast }={\mathbf{\Gamma }}_C^{-1}{\mathbf{\beta }}_C$$

$${\mathbf{\Gamma }}_C=({\mathbf{I}}_T\otimes \mathbf{C})\mathbf{\Gamma }({\mathbf{I}}_T\otimes {\mathbf{C}}^{\mathsf{T}})$$

$${\mathbf{\beta }}_C=({\mathbf{I}}_T\otimes \mathbf{C})\mathbf{\beta }.$$

$${\tilde{\mathbf{x}}}_{t,K}={\tilde{\mathbf{A}}}_K{\tilde{\mathbf{x}}}_{t-1,K}+{\mathbf{V}}_K^{\mathsf{H}}\text{diag}(\mathbf{c}){\hat{\mathbf{u}}}_{t-1}$$

Full text

I Introduction

The controllability of complex networks plays a fundamental role in our understanding of natural and technological systems. Relevant examples involve the control of social, biological, financial, road, communication, and smart grid networks. Different works have highlighted the importance of the network structure when controlling a system evolving on top of that network [2, 3, 4]. Other works controlled said system from a few control or driving nodes [5, 6, 7]. As an illustrative example, consider the Zachary’s Karate club social network in Figure 1 [8]. The network state may be an opinion profile (e.g., members thoughts on a topic) and controlling, or driving, the network amounts to shaping those opinions towards a desired or target state. The social relationships between members affect the ability to control the opinion profile and the objective is to sway all members opinions from a few influencing members (the driving nodes).

While providing seminal contributions on network controllability, the works in [2, 3, 4, 5, 6, 7] ignore the coupling between the underlying topology and the target state (a.k.a. the graph signal). Recent evidence from graph signal processing (GSP) [9, 10, 11] has shown that this coupling can bring substantial benefits in graph signal sampling [12], interpolation [13, 14], adaptive reconstruction [15], and observability of diffusion processes [16]. A common point that unifies [10, 11, 12, 13, 14, 15, 16] is the so-called graph Fourier transform (GFT). The GFT expands the network state onto an orthonormal basis related to the underlying topology —formed by the eigenvectors of a matrix that represents the network, such as the adjacency or the Laplacian matrix— akin to how the discrete Fourier transform expands a time signal on the complex exponential orthonormal basis. The GFT basis vectors can be linked to different variability modes across the graph through the concept of total variation [10, 17]; hence, named the graph oscillating modes.

A particular class of graph signals is that of bandlimited graph signals, i.e., signals that can be expressed by only a few graph oscillating modes. The number of active modes forms the graph signal bandwidth. Likewise for time signals, purely bandlimited graph signals rarely exist. But the bandlimitedness prior poses a powerful and parsimonious model to develop practical tools. This prior is exploited in GSP for sampling [10, 11, 18], where signals arising in economic [19], handwritten digits [12, 14] or brain functional imaging [20], are approximately bandlimited. In the Zachary’s Karate club example, a bandlimited network state corresponds to opinions polarized into a few clusters of like-minded members [12]. Therefore, controlling the system to a bandlimited state implies imposing a similar opinion profile to members that influence each other; this influence is captured by the edge weight.

Network control towards a bandlimited graph signal is considered in [21]. The control signal is fed to a few driving nodes and is percolated through the graph until the target state is reached. The authors determined the trade-off between the control time and the number of driving nodes, provided conditions to reach any bandlimited state, and designed the control signals.

The work in [22] studied the challenge of controlling the network towards a bandlimited state with control signals of limited energy. The main result is the trade-off between the number of driving nodes and the control signal energy. But conditions to drive the network to any bandlimited state were not derived in this limited energy setting.

Along these lines, and in parallel with the shorter version of this paper [1], the work in [23] used GSP to formulate the linear quadratic controller as an autoregressive moving average (ARMA) graph filter [24]. This ARMA formulation is used to control independently each graph mode to the desired state. But these control strategies require all nodes to act as driving nodes.

Altogether, the above works concern network control over time-invariant topologies. But in practice the network structure may change randomly over time due to link losses, or nodes that disappear with a given probability. This might be the case of: $i)$ members that are not present a given day in the Zachary’s Karate club; $ii)$ communication links that are random in sensor networks; $iii)$ power lines and buses that go down in smart grids due to local failures. In these situations, the graph structure is random and can be characterized by an expected graph with a variance around this expectation. Network controllability becomes, therefore, challenging since some connections cannot be exploited; hence, the derivations obtained in the deterministic setting may lead to a completely different state.

Motivated by the above observations, we study the possibilities to perform open loop control over random time varying networks. By exploring GSP tools as in [21, 22, 23], we propose a framework that accounts for the graph randomness in the analysis. Our main contributions are:

We study the problem of controlling dynamics over deterministic networks from a few driving nodes (Section III). We provide conditions on network controllability that relate the minimum number of driving nodes to the target control bandwidth and the control time. This result encompasses the three strategies proposed in [21]. 2. 2.

We formulate the problem of controlling dynamics over random time varying networks from a few driving nodes. (Section IV). We develop the concept of controllability in the mean to drive the expected state towards a graph signal with a target bandwidth content on the expected graph. We also extend the conditions on network controllability from the deterministic to the stochastic setting. 3. 3.

We perform a mean squared error analysis to quantify the statistical deviation of the controlled state from the target state on a particular graph realization (Section IV-B). This analysis illustrates the role of the random graph model, the target state bandwidth, and the control signal. 4. 4.

We propose two control strategies to drive the expected state to the desired bandlimited state with the minimum mean squared error (Section V). 5. 5.

We corroborate the developed framework and study its performance on synthetic (Erdős-Rényi and geometric graphs) and real-world social networks (Facebook sub-graph and Zachary’s Karate club) (Section VI).

To the best of our knowledge, this is the first contribution that approaches network controllability over random time varying topologies. We remark that while this work relies on the bandlimited prior to select the driving nodes, the reader might find of interest other parsimonious models [25, 18, 26] for such a task.

The remaining part of the paper proceeds as follows: Section II sets down the preliminary concepts and Section III contains our formulation of network controllability on deterministic graphs. Section IV formulates the controllability on random graphs. Section V develops the control strategies, while Section VI presents the numerical experiments. Finally, Section VII provides the concluding remarks. The proofs are collected in the appendix.

Notation. Normal letters $a$ (or $A$) denote scalars, bold lowercase letters ${\mathbf{a}}$ vectors, and bold uppercase letters ${\mathbf{A}}$ matrices. The $i$th entry of a vector is $[{\mathbf{a}}]_{i}$ while the $(i,j)$th entry of a matrix is $[{\mathbf{A}}]_{i,j}$. Superscripts $\mathsf{T}$ and $\mathsf{H}$ denote the transpose and the Hermitian, respectively. The $N\times 1$ null vector is ${\mathbf{0}}_{N}$, the $N\times 1$ vector of all ones is ${\mathbf{1}}_{N}$ and the $N\times N$ identity matrix is ${\mathbf{I}}_{N}$. The diagonal operator denoted as $\text{diag}(\cdot)$ is defined such that ${\mathbf{a}}=\text{diag}({\mathbf{A}})$ with $[{\mathbf{a}}]_{i}=[{\mathbf{A}}]_{i,i}$, and ${\mathbf{A}}=\text{diag}({\mathbf{a}})$ is a diagonal matrix with vector ${\mathbf{a}}$ on the main diagonal. The expectation operator is denoted as ${\mathbb{E}}[\cdot]$, the trace operator as $\text{tr}(\cdot)$, the rank of ${\mathbf{A}}$ as $\text{rank}({\mathbf{A}})$, the Kronecker product as $\otimes$, and the element-wise Hadamard product as $\circ$. The $l_{p}$ vector and matrix norm is denoted as $\|\cdot\|_{p}$. The ceiling operator is denoted as $\lceil\cdot\rceil$ and the minimum and maximum operators as $\min\{\cdot\}$ and $\max\{\cdot\}$, respectively. If not otherwise stated, calligraphic letters ${\mathcal{A}}$ indicate sets and the set cardinality is denoted as $|{\mathcal{A}}|$.

II Diffusion Processes on Graphs

In this work, we consider controlling a diffusion process on random time varying graphs towards a desired state. To achieve this, we model diffusion processes through graph signal processing. We introduce the basics of GSP in Section II-A, define the random time varying graph model in Section II-B, and discuss diffusion processes in Section II-C.

II-A Graph signal processing

Let ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathcal{W}})$ denote a graph with ${\mathcal{V}}=\{v_{1},\ldots,v_{N}\}$ the set of $N$ vertices, ${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$ the set of edges, and ${\mathcal{W}}:{\mathcal{E}}\to{\mathbb{R}}_{+}$ a function that assigns positive weights to the edges. The graph serves as a mathematical representation of the network and its structure is captured by the graph shift operator (GSO) matrix ${\mathbf{S}}\in{\mathbb{R}}^{N\times N}$. The $(i,j)$th element of ${\mathbf{S}}$, $[{\mathbf{S}}]_{i,j}$, is nonzero only if $i=j$ or if $(v_{j},v_{i})\in{\mathcal{E}}$; so that ${\mathbf{S}}$ respects the sparsity of ${\mathcal{G}}$. Standard choices for ${\mathbf{S}}$ are the weighted graph adjacency matrix ${\mathbf{W}}$ [27, 17], the graph Laplacian matrix ${\mathbf{L}}$ [10], or their respective generalizations [28, Chapter 8]. We consider ${\mathbf{S}}$ admits an eigendecomposition ${\mathbf{S}}={\mathbf{V}}\boldsymbol{\Lambda}{\mathbf{V}}^{\mathsf{H}}$, where ${\mathbf{V}}=[{\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{N}]\in{\mathbb{C}}^{N\times N}$ collects the orthonormal eigenvectors and $\boldsymbol{\Lambda}=\text{diag}(\lambda_{1},\ldots,\lambda_{N})\in{\mathbb{C}}^{N\times N}$ contains the associated eigenvalues. This holds for all undirected graphs on which the graph Laplacian can be defined and also for the adjacency matrix of some directed graphs [17, 19].

A graph signal is a mapping from the vertex set to the field of real numbers, i.e., $x_{i}:v_{i}\to{\mathbb{R}}$, for $v_{i}\in{\mathcal{V}}$. An example of a graph signal is the opinion of members in Zachary’s Karate club. We collect all node signals in the vector ${\mathbf{x}}\in{\mathbb{R}}^{N}$ with $[{\mathbf{x}}]_{i}=x_{i}$ being the value of node $i$ [10].

The graph Fourier transform (GFT) is the projection of the graph signal ${\mathbf{x}}$ on the eigenbasis ${\mathbf{V}}$ and is denoted by ${\tilde{\mathbf{x}}}={\mathbf{V}}^{\mathsf{H}}{\mathbf{x}}$ [9, 10]. The elements $[{\tilde{\mathbf{x}}}]_{k}={\tilde{x}}_{k}$ denote the graph Fourier coefficients of ${\mathbf{x}}$, whereas the eigenvectors ${\mathbf{v}}_{k}$ form the basis of graph oscillating modes. Likewise, the inverse GFT is ${\mathbf{x}}={\mathbf{V}}{\tilde{\mathbf{x}}}$, i.e., it writes ${\mathbf{x}}$ as a linear combination of the graph oscillating modes weighed by the Fourier coefficients.

A graph signal is bandlimited if it has only a few nonzero Fourier coefficients. Without loss of generality, assume the first $K$ elements of ${\tilde{\mathbf{x}}}$ are nonzero; so, we can write ${\tilde{\mathbf{x}}}=[{\tilde{\mathbf{x}}}_{K}^{\mathsf{T}},{\mathbf{0}}_{N-K}^{\mathsf{T}}]^{\mathsf{T}}$ where ${\tilde{\mathbf{x}}}_{K}\in{\mathbb{C}}^{K}$ and $K\leq N$. Then, ${\mathbf{x}}$ is written in the compact form

[TABLE]

where ${\mathbf{V}}_{K}\in{\mathbb{C}}^{N\times K}$ is the respective column-trimmed eigenvector matrix. The GFT ${\tilde{\mathbf{x}}}_{K}$ of ${\tilde{\mathbf{x}}}$ writes as ${\tilde{\mathbf{x}}}_{K}={\mathbf{V}}_{K}^{\mathsf{H}}{\mathbf{x}}$. This representation connects the signal bandwidth with the sampling and reconstruction strategies as shown in [29, 12, 30, 15, 14]. We will also exploit bandlimitedness in Section IV to control the network from a few driving nodes.

II-B Random time varying graphs

We consider the following random graph model.

Definition 1** ($\text{RES}(p)$ graph model [31]).**

Given an underlying graph ${\mathcal{G}}\!=\!({\mathcal{V}},{\mathcal{E}})$, a random edge sampling (RES) graph realization ${\mathcal{G}}_{t}\!=\!({\mathcal{V}},{\mathcal{E}}_{t})$ of ${\mathcal{G}}$ consists of the same set of nodes ${\mathcal{V}}$ and assumes the edge $(v_{i},v_{j})\in{\mathcal{E}}$ is sampled at time $t$ (i.e., $(v_{i},v_{j})\in{\mathcal{E}}_{t}$) with a probability $0<p\leq 1$. The edges are sampled independently over both the graph and the temporal dimension and are mutually independent from the graph signal if the latter has a stochastic nature.

In other words, the $\text{RES}(p)$ model states that the realization ${\mathcal{G}}_{t}=({\mathcal{V}},{\mathcal{E}}_{t})$ is drawn from the underlying graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$, where the instantaneous edge set ${\mathcal{E}}_{t}\subseteq{\mathcal{E}}$ is generated via an independent Bernoulli process with probability $p$. Let us from now on denote with ${\mathbf{W}}$ and ${\mathbf{D}}=\text{diag}({\mathbf{W}}{\mathbf{1}}_{N})$ the adjacency and degree matrix of ${\mathcal{G}}$, respectively. If the graph is undirected, we will also consider the unnormalized graph Laplacian matrix ${\mathbf{L}}={\mathbf{D}}-{\mathbf{W}}$. To ease the exposition, denote with ${\mathbf{W}}_{t}$, ${\mathbf{D}}_{t}$, and ${\mathbf{L}}_{t}$ the respective matrices of the instantaneous graph ${\mathcal{G}}_{t}$ and with ${\bar{\mathbf{W}}}={\mathbb{E}}[{\mathbf{W}}_{t}]$, ${\bar{\mathbf{D}}}={\mathbb{E}}[{\mathbf{D}}_{t}]$, and ${\bar{\mathbf{L}}}={\mathbb{E}}[{\mathbf{L}}_{t}]$ those of the expected graph $\bar{{\mathcal{G}}}$. Under the $\text{RES}(p)$ model, it holds that ${\bar{\mathbf{W}}}=p{\mathbf{W}}$, ${\bar{\mathbf{D}}}=p{\mathbf{D}}$, and ${\bar{\mathbf{L}}}=p{\mathbf{L}}$.

We assume the following.

Assumption 1**.**

The GSO of the underlying graph ${\mathcal{G}}$ has an upper bounded spectral norm $\|{\mathbf{S}}\|_{2}\leq\varrho$ for some $\varrho<\infty$.

This assumption is generally met in practice and implies the graphs of interest have finite dimension and edge weights.

We also remark that more complex models than the $\text{RES}(p)$ can be found in literature. In particular, many of these results can be readily extended to the model in which each edge $(v_{i},v_{j})$ is sampled independently with a different probability $p_{ij}$. But for clarity of exposition, we will focus only on the $\text{RES}(p)$ model.

II-C Diffusion on graphs from a GSP perspective

The continuous-time diffusion of a signal ${\mathbf{x}}_{0}$ on a graph ${\mathcal{G}}$ with Laplacian matrix ${\mathbf{L}}$ is described by the differential equation [32, 33]

[TABLE]

This equation can be discretized as [34]

[TABLE]

which is stable if $0<\epsilon\leq 1/\|{\mathbf{L}}\|_{2}$. Alternatively, a diffusion process on a graph can be interpreted as the discrete-time shift of ${\mathbf{x}}_{0}$ through the graph edges [27, 35]

[TABLE]

Model (3) is used when the process is defined over a continuous space that has been discretized, usually in the form of a mesh, as for heat diffusion processes [33]. Model (4) is employed when the underlying support is naturally a graph, as for sensor network communications [35]. In essence, these are two examples of processes that describe the network state evolution by ${\mathbf{x}}_{t}={\mathbf{A}}{\mathbf{x}}_{t-1}$ and relate this state to the underlying time-invariant topology (the transition matrix ${\mathbf{A}}$ depends on the shift operator). In this paper, we consider the more general case of random time varying topologies, i.e. ${\mathbf{x}}_{t}={\mathbf{A}}_{t-1}{\mathbf{x}}_{t-1}$, and we abstract the relationship between the transition matrix and the underlying topology as follows.

Assumption 2**.**

Let ${\mathbf{A}}_{t}$ be the time varying transition matrix of a diffusion process over a random time varying graph ${\mathcal{G}}_{t}\in\text{RES}(p)$. Then, ${\mathbb{E}}[{\mathbf{A}}_{t}]$ and ${\mathbf{S}}$ share the same eigenvectors.

That is, we consider diffusions on random graphs such that the eigenvectors of the expected transition matrix and the underlying GSO coincide. The following lemma shows this is the case for the diffusion models (3) and (4) on $\text{RES}(p)$ graph realizations.

Lemma 1**.**

Let ${\mathcal{G}}$ be a graph satisfying Assumption 1 and let ${\mathcal{G}}_{t}$ be a $\text{RES}(p)$ realization of it. For the diffusion models

1. (i)

${\mathbf{S}}={\mathbf{L}}$* and ${\mathbf{A}}_{t}={\mathbf{I}}-\epsilon{\mathbf{L}}_{t}$ [cf. (3)],* 2. (ii)

${\mathbf{S}}={\mathbf{W}}$* and ${\mathbf{A}}_{t}={\mathbf{W}}_{t}$ [cf. (4)],*

Assumption 2 and $\|{\mathbf{A}}_{t}\|_{2}\leq\varrho$ hold.

Other models that satisfy these conditions are the wave equation on graphs and graph-based ARMA models, see [16].

III Controllability on Deterministic Graphs

Consider the $N$-state linear system

[TABLE]

where ${\mathbf{x}}_{t}\in{\mathbb{R}}^{N}$ denotes the state value on all nodes at time $t$, ${\mathbf{u}}_{t}\in{\mathbb{R}}^{M}$ is the control signal injected on $M\leq N$ nodes, and ${\mathbf{A}}\in{\mathbb{R}}^{N\times N}$ and ${\mathbf{B}}\in{\mathbb{R}}^{N\times M}$ are the transition and control input matrix, respectively. The relationship between the network state ${\mathbf{x}}_{t}$ and the underlying topology is captured in (5) through the transition matrix ${\mathbf{A}}$; it shares the eigenvectors with ${\mathbf{S}}$ and this is the case for models (3) and (4).

System (5) is controllable if and only if the controllability matrix

[TABLE]

has full rank $N$ [36, Section 6.2.1]. While full rank of $\boldsymbol{\Omega}$ guarantees the convergence of ${\mathbf{x}}_{t}$ to any target signal ${\mathbf{x}}^{*}$, we focus on controlling the network state towards a bandlimited graph signal ${\mathbf{x}}^{\ast}={\mathbf{V}}_{K}{\tilde{\mathbf{x}}}_{K}^{\ast}$. Here, ${\tilde{\mathbf{x}}}_{K}^{\ast}\in{\mathbb{C}}^{K}$ determines the desired frequency response over the $K$ frequencies of interest; the target bandwidth. We thus define the bandwidth controllability as follows.

Definition 2** (Bandwidth controllability).**

An $N$-state system on a graph is bandwidth controllable from $M\leq N$ nodes if, for any initial state ${\mathbf{x}}_{0}$ and some final time $T$, there exists a sequence of control signals $\{{\mathbf{u}}_{t}\in{\mathbb{R}}^{M},t=0,1,\ldots,T-1\}$ acting on a fixed set of $M$ nodes that drive the network state to a value ${\mathbf{x}}^{*}$ with any frequency content ${\tilde{\mathbf{x}}}_{K}^{\ast}={\mathbf{V}}_{K}^{\mathsf{H}}{\mathbf{x}}^{\ast}$ over the $K\leq N$ target bandwidth.

Lead by the promising results of bandlimited graph signal reconstruction from samples on a few nodes [29, 12, 30, 15, 14], we aim to control ${\mathbf{x}}_{t}$ through a fixed, time-invariant, set of nodes ${\mathcal{S}}$ of cardinality $|{\mathcal{S}}|=M\leq N$. Let then ${\mathbf{B}}={\mathbf{C}}^{\mathsf{T}}$ denote a binary matrix that selects these nodes. More formally, ${\mathbf{C}}$ belongs to the combinatorial set

[TABLE]

that selects $M$ out of $N$ different nodes and the ordering relation $\leq$ among vectors stands for the elementwise partial ordering [39, Example 2.23]. Observe that ${\mathbf{C}}{\mathbf{C}}^{\mathsf{T}}={\mathbf{I}}_{M}$ and ${\mathbf{C}}^{\mathsf{T}}{\mathbf{C}}=\text{diag}({\mathbf{c}})$ with ${\mathbf{c}}\in\{0,1\}^{N}$, such that $[{\mathbf{c}}]_{i}=1$ if and only if node $v_{i}$ belongs to ${\mathcal{S}}$.

With this in place, we write the linear system on graphs (5) in the GFT domain as

[TABLE]

where ${\tilde{\mathbf{A}}}={\mathbf{V}}^{\mathsf{H}}{\mathbf{A}}{\mathbf{V}}$ is a diagonal matrix containing the eigenvalues of ${\mathbf{A}}$ [cf. Assumption 2]. For convenience, we write ${\tilde{\mathbf{A}}}=\text{diag}({\mathbf{a}})\in{\mathbb{C}}^{N\times N}$ with ${\mathbf{a}}\in{\mathbb{C}}^{N}$ the vector containing the eigenvalues of ${\tilde{\mathbf{A}}}$, known also as the spectral response of ${\mathbf{A}}$. Then, by splitting (8) into the $K$ frequencies of interest we have

[TABLE]

where ${\tilde{\mathbf{x}}}_{t}=[{\tilde{\mathbf{x}}}_{t,K}^{\mathsf{T}},{\tilde{\mathbf{x}}}_{t,N-K}^{\mathsf{T}}]^{\mathsf{T}}$, ${\mathbf{a}}=[{\mathbf{a}}_{K}^{\mathsf{T}},{\mathbf{a}}_{N-K}^{\mathsf{T}}]^{\mathsf{T}}$ and ${\mathbf{V}}=[{\mathbf{V}}_{K},{\mathbf{V}}_{N-K}]$.

For a system that is bandwidth controllable [cf. Def. 2], recursion (9) can drive the network state ${\mathbf{x}}_{t}$ to any signal ${\mathbf{x}}_{T}$ that has a desired frequency response ${\tilde{\mathbf{x}}}_{K}^{\ast}$ over the target bandwidth $K$. Thus, we can focus on the $K$ frequencies of interest, determine the driving nodes through matrix ${\mathbf{C}}$, and design the control signals $\{{\mathbf{u}}_{t}\}$ such that the system

[TABLE]

reaches ${\tilde{\mathbf{x}}}_{K}^{\ast}$ at time $T$, i.e. ${\tilde{\mathbf{x}}}_{T,K}={\tilde{\mathbf{x}}}_{K}^{\ast}$. But we observe that focusing on these $K$ frequencies leads to a non-zero value also on the $N-K$ remaining frequencies. One way to suppress this undesired content is to force ${\tilde{\mathbf{x}}}_{t,N-K}={\mathbf{0}}_{N-K}$. This requires the design of a sampling matrix ${\mathbf{C}}$ that satisfies ${\mathbf{V}}_{N-K}^{\mathsf{H}}{\mathbf{C}}^{\mathsf{T}}={\mathbf{0}}_{(N-K)\times M}$ [cf. (9)]. The latter might be infeasible or might severely constraint the selection of driving nodes. Therefore, to avoid the out-of-band frequency content, we use a frequency-selective graph filter ${\mathbf{H}}={\mathbf{V}}_{K}{\mathbf{V}}_{K}^{\mathsf{H}}$ to force bandlimitedness on the final network state ${\mathbf{x}}^{\ast}={\mathbf{H}}{\mathbf{x}}_{T}$. This filter can be implemented locally through a polynomial in the shift operator ${\mathbf{S}}$ with degree at most $N$ [37].

Hereinafter, the design variables are the sampling matrix ${\mathbf{C}}$ and the control signals ${\mathbf{u}}_{t}\in{\mathbb{R}}^{M}$ for $t=0,\ldots,T-1$. And, since the initial state ${\mathbf{x}}_{0}$ is considered known we fix, without loss of generality, ${\mathbf{x}}_{0}={\mathbf{0}}_{N}$ as the common practice in control literature [22], [36, Section 2.3.2], [38, Section 2.1]. With this set down, we claim our first contribution, which we will also exploit for controllability on random graphs in Section IV.

Proposition 1**.**

Consider the linear system (10) describing a process over a deterministic graph ${\mathcal{G}}$. A necessary condition to control the system in a finite time $T$ towards a target frequency content ${\tilde{\mathbf{x}}}^{*}_{K}$ over the $K$ frequencies of interest is to select $M\geq\lceil K/T\rceil$ driving nodes.

Proposition 1 provides a necessary condition on the minimum number of nodes to control system (10) in $T$ instants. It shows the trade-off between the cardinality of the sampling set $M$, the signal bandwidth $K$, and the control time $T$. Thus, for $T\geq K$ there is the potential to control the network by acting on a single node. But this is not sufficient since the system controllability is affected by the network topology [2, 3, 4], i.e. the influence of the driving nodes on the frequency content, see also [12, 19, 21]. Hence, the driving nodes should be carefully picked to guarantee the controllability of (10). In what follows, we show the relation of Proposition 1 with [21] and [22].

a) Relation with [21]. Proposition 1 encompasses under a single condition the three graph signal reconstruction strategies of [21]. In fact, $M\geq K$ and $T=1$ covers the multiple node-single time seeding strategy; $M=1$ and $T\geq K$ is a necessary condition to control the signal for the single node-multiple time seeding strategy; and $M\geq\lceil K/T\rceil$ covers the more involved multiple node-multiple time seeding approach. This is expected, since graph signal reconstruction through percolation is a particular case of system (5).

b) Relation with [22]. Differently from our approach, [22] focuses on designing the control signal ${\mathbf{u}}=[{\mathbf{u}}_{T-1}^{\mathsf{T}},\ldots,{\mathbf{u}}_{1}^{\mathsf{T}},{\mathbf{u}}_{0}^{\mathsf{T}}]^{\mathsf{T}}\in{\mathbb{R}}^{MT\times 1}$ as a trade-off between sparsity in the vertex domain and signal energy. This problem writes as

[TABLE]

where the constant $\gamma$ trades the control signal energy $\|{\mathbf{u}}\|_{2}^{2}$ with sparsity $\|{\mathbf{u}}\|_{0}$. Problem (11) yields a sparse control signal ${\mathbf{u}}$ across time, but the driving nodes are not necessarily fixed. In this regard, Proposition 1 imposes a minimum dimension on ${\mathbf{u}}$ such that controllability is possible from a fixed set ${\mathcal{S}}$ of driving nodes.

IV Controllability on Random Graphs

When the network topology is time varying, system (5) should change to reflect the time dependency in the transition matrix. When the time variation is random, the controllability of the network state should follow a statistical approach. We propose a statistical framework in this section, where in Section IV-A we develop the concept of controllability in the mean and in Section IV-B we perform the mean squared error analysis.

IV-A Mean controllability

The dynamics of a time varying system on random graphs are given by

[TABLE]

where under the $\text{RES}(p)$ model in Definition 1, $\{{\mathbf{A}}_{t}\}$ is a set of i.i.d. random matrices with ${\mathbb{E}}[{\mathbf{A}}_{t}]={\bar{\mathbf{A}}}$. The deterministic design variables are contained in the second term of (12), ${\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{t-1}$. The state ${\mathbf{x}}_{t}$ depends on the random system matrices $\{{\mathbf{A}}_{\tau}\}_{\tau=0}^{t-1}$, which are independent from ${\mathbf{A}}_{t}$ and from the deterministic design variables ${\mathbf{C}}$ and $\{{\mathbf{u}}_{\tau}\}_{\tau=0}^{t-1}$. Note also that ${\bar{\mathbf{A}}}$ and ${\mathbf{S}}$ share the same eigenvectors, i.e. ${\bar{\mathbf{A}}}={\mathbf{V}}\text{diag}({\bar{\mathbf{a}}}){\mathbf{V}}^{\mathsf{H}}$; thus, it captures in statistics the relation between the underlying topology and state ${\mathbf{x}}_{t}$. We can then write the mean evolution of (12) as

[TABLE]

where $\boldsymbol{\mu}_{t}={\mathbb{E}}[{\mathbf{x}}_{t}]$. System (13) is a deterministic and time-invariant system analogous to (5). We develop the following controllability concept.

Definition 3** (Bandwidth controllability in the mean).**

An $N$-state system on a random graph of the form in (12) with mean evolution in (13) is bandwidth controllable in the mean from $M\leq N$ nodes if, for any initial state ${\mathbf{x}}_{0}$ and some final time $T$, there exists a sequence of control signals $\{{\mathbf{u}}_{t},t=0,\ldots,T-1\}$ acting on a fixed set of $M$ nodes that drive the mean network state to a value ${\mathbf{x}}^{*}$ with any frequency content ${\mathbf{x}}_{K}^{\ast}={\mathbf{V}}_{K}^{\mathsf{H}}{\mathbf{x}}^{\ast}$ over the $K\leq N$ target bandwidth.

Our goal is to control the mean system to a desired bandlimited graph signal ${\mathbf{x}}^{\ast}={\mathbf{V}}_{K}{\tilde{\mathbf{x}}}_{K}^{\ast}$ in a finite time $T$ from a few nodes. We do so by designing the input signals $\{{\mathbf{u}}_{t},t=0,\ldots,T-1\}$ and the node driving set ${\mathcal{S}}$ (through matrix ${\mathbf{C}}$). In analogy to Section III, we focus on the $K$ frequencies of interest of the mean system [cf. (10)]

[TABLE]

and drive it to the desired frequency content ${\tilde{\mathbf{x}}}_{K}^{\ast}$ [cf. Def. 3], with $\check{{\mathbf{A}}}_{K}=\widetilde{{\bar{\mathbf{A}}}}_{K}=\text{diag}({\bar{\mathbf{a}}}_{K})\in{\mathbb{C}}^{K\times K}$ containing the eigenvalues of ${\bar{\mathbf{A}}}$ which determine the spectral response of the system evolution on the expected graph. Then, we apply a (deterministic) linear filter ${\mathbf{H}}={\mathbf{V}}_{K}{\mathbf{V}}_{K}^{\mathsf{H}}$ to keep only those $K$ desired frequencies such that the mean network state $\boldsymbol{\mu}^{\ast}={\mathbb{E}}[{\mathbf{H}}{\mathbf{x}}_{T}]={\mathbf{H}}\boldsymbol{\mu}_{T}={\mathbf{x}}^{\ast}$ results in a bandlimited graph signal.

Similarly to Proposition 1, we claim the following.

Proposition 2**.**

Consider the linear system (12) describing a process over a sequence of $\text{RES}(p)$ graphs ${\mathcal{G}}_{t}$ with in-band mean evolution (14). A necessary condition to control the mean system in finite time $T$ towards a target frequency content ${\tilde{\mathbf{x}}}^{*}_{K}$ over the $K$ frequencies of interest is to select $M\geq\lceil K/T\rceil$ driving nodes.

Like Proposition 1, Proposition 2 establishes a necessary condition to control a linear system, now, on random time varying graphs. As such, the same trade-off between the number of driving nodes $M$, the signal bandwidth $K$, and the control time $T$ applies here. The next corollary extends this result to a sufficient condition under some restrictions on the eigenvector basis.

Corollary 1**.**

Under the hypothesis of Proposition 2, if there exists a set of driving nodes ${\mathcal{S}}$ built by $M$ nodes such that the corresponding $M$ rows of ${\mathbf{V}}_{K}$ are linearly independent vectors, then $M\geq K$ is a sufficient condition to control system (12) in the mean.

An algorithm for finding such $M$ nodes is readily available in [12, Algorithm 1].

IV-B Mean squared error analysis

In Section IV-A, we discussed that system (12) can be controlled in the mean. Therefore, it is paramount also to quantify the mean squared error (MSE) of the controlled state to gain statistical insight into how close the filtered final state on a specific graph realization ${\mathbf{H}}{\mathbf{x}}_{T}$ is to the actual desired signal ${\mathbf{x}}^{\ast}$. Towards this end, define $\boldsymbol{\Phi}_{b,a}={\mathbf{A}}_{b}{\mathbf{A}}_{b-1}\cdots{\mathbf{A}}_{a+1}{\mathbf{A}}_{a}$ as the state transition matrix between time instants $b\geq a$. The following theorem determines the MSE.

Theorem 1**.**

Let Assumptions 1 and 2 hold and let ${\mathbf{x}}_{t}$ be a graph process defined over a sequence of $\text{RES}(p)$ graphs ${\mathcal{G}}_{t}$ described by linear system (12). Given also a set of driving nodes ${\mathcal{S}}$ characterized by selection matrix ${\mathbf{C}}$ and a set of control signals $\{{\mathbf{u}}_{t}\}_{t=0}^{T-1}$ with initial state ${\mathbf{x}}_{0}={\mathbf{0}}_{N}$. The MSE between the filtered signal ${\mathbf{H}}{\mathbf{x}}_{T}$ on a particular graph realization and the actual desired signal is

[TABLE]

which is a quadratic form on ${\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{\tau}$ with coefficients $\alpha=\|{\mathbf{x}}^{\ast}\|_{2}^{2}\in{\mathbb{R}}$, $\boldsymbol{\beta}_{\tau}=({\bar{\mathbf{A}}}^{T-\tau-1})^{\mathsf{T}}{\mathbf{H}}^{\mathsf{T}}{\mathbf{x}}^{\ast}\in{\mathbb{R}}^{N\times 1}$ and $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}={\mathbb{E}}[\boldsymbol{\Phi}_{T-1,\tau+1}{\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}\boldsymbol{\Phi}_{T-1,\tau^{\prime}+1}]\in{\mathbb{R}}^{N\times N}$.

Given ${\mathbf{C}}$ and $\{{\mathbf{u}}_{t}\}$, the MSE in (15) holds for any system described by (12), irrespective of their controllability. The MSE in (15) is a quadratic function in the design variables ${\mathbf{C}}^{\mathsf{T}}{\mathbf{u}}_{\tau}$ and the corresponding coefficients $\alpha$, $\boldsymbol{\beta}_{\tau}$, and $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}$ depend on known quantities: the graph filter ${\mathbf{H}}$; the desired state ${\mathbf{x}}^{\ast}={\mathbf{V}}_{K}{\tilde{\mathbf{x}}}_{K}^{\ast}$; and the statistics (first and second order moments) of the underlying support through ${\bar{\mathbf{A}}}$ and $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}$. This result highlights also the impact the driving nodes have on the overall performance and shows their connection with the underlying support and the target bandwidth. More precisely, the coefficient $\alpha$ provides the MSE floor if there is no control signal (i.e., given by the energy of the target state); $\boldsymbol{\beta}_{\tau}$ takes into account the similarity between the target signal and the controlled signal evolution over the mean graph; and $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}$ accounts for the variability of the random graph. Finally, we note that the computation of $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}$ in (15) might be cumbersome for some transition matrices ${\mathbf{A}}_{t}$. We thus provide in the appendix two practical results that address this issue: first, we provide a general upper bound; second, we show how to exactly compute $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}$ for undirected graphs for the diffusion models in Lemma 1.

V Control Strategies

In this section, we propose two control strategies (i.e., find ${\mathbf{C}}$ and $\{{\mathbf{u}}_{t}\}$) for graph processes over random time varying graphs, where, depending on the scenario, one can be preferred over the other. In Section V-A we propose an unbiased control strategy, while in Section V-B we introduce a control strategy that leverages the bias-variance trade-off to minimize the MSE.

V-A Unbiased controller

The mean state in (14) for $t=T$ can be expanded as

[TABLE]

where $\check{{\mathbf{A}}}_{K}=\widetilde{{\bar{\mathbf{A}}}}_{K}=\text{diag}({\bar{\mathbf{a}}}_{K})\in{\mathbb{C}}^{K\times K}$ and $\tilde{\boldsymbol{\mu}}_{0,K}={\mathbf{0}}_{K}$. For an unbiased controller, it must hold at final time $T$ $\tilde{\boldsymbol{\mu}}_{T,K}={\tilde{\mathbf{x}}}_{K}^{\ast}$. Combining (16) and $\tilde{\boldsymbol{\mu}}_{T,K}={\tilde{\mathbf{x}}}_{K}^{\ast}$, we obtain

[TABLE]

with in-band controllability matrix $\tilde{\boldsymbol{\Omega}}\in{\mathbb{C}}^{K\times TM}$ and input vector is ${\mathbf{u}}=[{\mathbf{u}}_{T-1}^{\mathsf{T}},{\mathbf{u}}_{T-2}^{\mathsf{T}},\ldots,{\mathbf{u}}_{1}^{\mathsf{T}},{\mathbf{u}}_{0}^{\mathsf{T}}]^{\mathsf{T}}{\mathbb{R}}^{MT\times 1}$. For $\tilde{\boldsymbol{\Omega}}$ being of full rank $K$ (i.e. a controllable system), system (17) has infinite solutions on ${\mathbf{u}}$. Also often exists more than one set of nodes that guarantees controllability. We then select the set of nodes and design the control signals to minimize the $\mathrm{MSE}(T)$, while guaranteeing the solution is unbiased.

Let ${\mathcal{C}}^{\ast}_{M,N}=\{{\mathbf{C}}\in{\mathcal{C}}_{M,N}:\text{rank}(\tilde{\boldsymbol{\Omega}})=K\}$ be the set of selection matrices that satisfy controllability. The optimal unbiased control strategy can be posed as

[TABLE]

where $\mathrm{MSE}(T)$ is given in (15). Oftentimes, we are interested in controlling the system with minimum energy [7, 22]. In such cases, the minimum energy control signal is [39, Section 6.2]

[TABLE]

Then, within the minimum energy framework, we select the nodes that minimize the $\mathrm{MSE}(T)$ as follows

[TABLE]

Problem (20) is non-convex due to the binary nature of the optimization variable ${\mathbf{C}}$. A heuristic solution is to follow a constrained greedy approach as described in Algorithm 1. The objective is to greedily select the nodes that improve the $\mathrm{MSE}$ while satisfying controllability. Specifically, for each candidate driving node, we check $\text{rank}(\tilde{\boldsymbol{\Omega}})$ increases until we reach controllability as indicated in line 11 [cf. Proposition 2]. Since we are looking for the minimum energy controller, then line 12 entails computing ${\mathbf{u}}^{\ast}=(\tilde{\boldsymbol{\Omega}}^{\mathsf{H}}\tilde{\boldsymbol{\Omega}})^{-1}\tilde{\boldsymbol{\Omega}}^{\mathsf{H}}{\tilde{\mathbf{x}}}_{K}^{\ast}$ if $\text{rank}(\tilde{\boldsymbol{\Omega}})=Tm$, and (19) if $\text{rank}(\tilde{\boldsymbol{\Omega}})=K$. While this constrained greedy approach has no theoretical guarantees [40], our numerical results in Section VI show that Algorithm 1 exhibits a performance close to the optimal solution.

V-B Biased controller

When the requirement for an unbiased controller is not strict, we can leverage the bias-variance trade-off to further reduce the mean squared error of the controlled state. Given a fixed sampling set ${\mathbf{C}}$, the $\mathrm{MSE}(T)$ (15) is a quadratic function on the control signals $\{{\mathbf{u}}_{t},t=0,1,\ldots,T-1\}$. Therefore, we express ${\mathbf{u}}_{t}$ as a function of ${\mathbf{C}}$ as follows.

First, the derivative of $\mathrm{MSE}(T)$ in (15) w.r.t. ${\mathbf{u}}_{t}$ is

[TABLE]

for $t=0,\ldots,T-1$ with $\boldsymbol{\beta}_{t}$ and $\boldsymbol{\Gamma}_{t,\tau}$ given in (15). By defining $\boldsymbol{\beta}=[\boldsymbol{\beta}_{T-1}^{\mathsf{T}},\ldots,\boldsymbol{\beta}_{0}^{\mathsf{T}}]^{\mathsf{T}}\in{\mathbb{R}}^{NT}$,

[TABLE]

and by setting (21) to zero, we can write

[TABLE]

with $\boldsymbol{\Gamma}_{C}\in{\mathbb{R}}^{MT\times MT}$ and $\boldsymbol{\beta}_{C}\in{\mathbb{R}}^{MT\times 1}$. By construction $\boldsymbol{\Gamma}_{C}$ has rank $MT$ since ${\mathbf{C}}\in{\mathcal{C}}_{M,N}$. Then, for a sampling matrix ${\mathbf{C}}$ such that $\text{rank}(\boldsymbol{\Gamma}_{C})\geq MT$, $\boldsymbol{\Gamma}_{C}$ is nonsingular and leads to the (parameterized) minimum $\mathrm{MSE}(T)$ control signals

[TABLE]

From this relation between the control signal and the sampling matrix, we consider a two-stage optimization approach [39, Section 4.1.3] to find ${\mathbf{C}}$ that minimizes the $\mathrm{MSE}(T)$. This optimization problem writes as

[TABLE]

To deal with the non-convexity of (25), similarly to (20), we rely on a constrained greedy approach analogous to Algorithm 1. Specifically, we replace line 10 by the computation of $\boldsymbol{\Gamma}_{C}$ and $\boldsymbol{\beta}_{C}$, line 11 by $\text{rank}(\boldsymbol{\Gamma}_{C})\geq mT$, and line 12 by (24).

VI Numerical Experiments

We evaluate the proposed control strategies on different scenarios to analyze the different trade-off when controlling the network. We compare the unbiased minimal energy controller (20) and the biased controller (25) with the Percolation control strategy of [21] and with the Min. Energy approach of [22]. Next, we consider synthetic network models, namely Erdős-Rényi (ER) graphs [41] and geometric graphs, while in Section VI-B we test the methods on real-world social networks, namely on the Zachary’s Karate Club [8] and on a Facebook subnet [42].

VI-A Synthetic network models

The ER graph forms the edges between any two nodes randomly and independently with probability $p_{\textrm{ER}}$ and has an average degree of $p_{\textrm{ER}}N$. The geometric graph draws nodes uniformly at random in the $[0,1]^{2}$ plane and computes the Euclidean distance $d_{ij}$ between any pair of nodes. We assigned a Gaussian kernel edge weights $w_{ij}={\mathcal{W}}(v_{i},v_{j})=e^{-d_{ij}^{2}}$ and kept only the $k_{\textrm{NN}}$ nearest neighbors per node; the parameter $k_{\textrm{NN}}$ controls the average degree. For both models, we considered realizations that result in connected graphs. To account for the randomness in the generative models and in the edge loss, we averaged the performance over $500$ different underlying graphs where for each of them we accounted also for $5000$ RES realizations.

Unless otherwise specified, we set $N=100$, $p_{\textrm{ER}}=0.5$, $k_{\textrm{NN}}=5$, and the RES link loss probability to $p_{\text{RES}}=0.95$. The control time is $T=8$, the number of driving nodes is $M=8$, and the initial state is ${\mathbf{x}}_{0}={\mathbf{0}}_{N}$. The target state ${\mathbf{x}}^{\ast}$ has a bandwidth of $K=10$ with GFT coefficients ${\tilde{\mathbf{x}}}_{K}^{\ast}$ decaying linearly as $[{\tilde{\mathbf{x}}}_{K}^{\ast}]_{k}=1-(k-1)/K$ for $k=1,\ldots,K$. We measured the controllability performance between the bandlimited controlled state and target one through the normalized MSE: $\mathrm{MSE}(T)={\mathbb{E}}[\|{\mathbf{H}}{\mathbf{x}}_{T}-{\mathbf{x}}^{\ast}\|^{2}]/\|{\mathbf{x}}^{\ast}\|^{2}$.

Time-invariant network. To set a baseline, we first compared the biased (25) and the unbiased (20) controllers with the Percolation [21] and Min. Energy [22] strategies on a fixed time-invariant network. This is the same control scenario that is considered in [21, 22] and is equivalent to setting $p_{\text{RES}}=1$. We report the results for geometric graphs in Figure 2. Figure 2(a) is a parametric simulation as a function of time horizon $T$ and Figure 2(b) is a parametric simulation as a function of the number of samples. In general, we observe that the biased estimator has a performance similar to the Min. Energy and slightly better than Percolation. The unbiased controller lags behind in terms of $\mathrm{MSE}$.

Graph connectivity. In the first random time varying experiment, we studied the impact of the graph connectivity on the controllability performance. We accounted for the graph connectivity by changing the average degrees, i.e., $p_{\textrm{ER}}N$ for the ER graph and $k_{\textrm{NN}}$ in the geometric graph. Figure 3 shows the MSE as the connectivity increases. From Figure 3(a), we observe that the control on geometric networks improves with the average degree. This is intuitively satisfying since larger degrees lead to a higher connectivity between nodes; hence, it renders them more robust to the RES model. Contrarily, for the ER model in Figure 3(b) this behavior is not as much emphasized. We attribute this phenomenon to the large average degree of the ER graphs (above $30$) and to the relatively high value of $p_{\textrm{RES}}$. That is, the loss of a few edges does not impact the overall ability to control the network.

Link loss. In the second experiment, we analyzed the impact of $p_{\text{RES}}$ for a fixed average degree. From Figure 4, we note that as $p_{\textrm{RES}}$ increases (fewer links are lost) the MSE reduces and leads to an easier to control network. This is because a higher $p_{\text{RES}}$ yields realizations with fewer edge losses, thus more similar to the underlying (mean) graph.

Control time. In the third and last experiment, we analyzed the impact of the control time horizon $T$. From Figure 5, we observe that the proposed strategies are not significantly affected by changes in $T$ as they only improve slightly.

From this set of experiments, we make three key observations. First, the proposed strategies offer the best performance. Second, the biased controller achieves the lowest MSE. This is expected since it levers the bias-variance trade-off to minimize the overall MSE at expenses of a bias in the controlled state. Third, not accounting for the graph randomness affects seriously the performance, even for $p_{\text{RES}}=0.999$. In fact, the deterministic alternatives of Percolation and Min. Energy have a worse performance by orders of magnitude compared with the proposed techniques. This contrast is particularly evident when comparing with the simulations for a time-invariant network in Figure 2. This could be explained by the fact that losing a link has a huge impact in the topology of the graph and severely affects the eigenbasis, thereby, changing the subspace of signals that are bandlimited on a given graph.

VI-B Real world graphs

We consider the formation of opinion profiles on two social networks, namely the Zachary’s Karate Club [8] in Figure 1 of $N=34$ nodes and a Facebook subnetwork [42] of $N=234$ nodes.

We set the control time to $T=8$, the number of selected nodes to $M=\text{round}(0.08N)$, $p_{\text{RES}}=0.95$, and ${\mathbf{x}}_{0}={\mathbf{0}}_{N}$. For simplicity of presentation, we focus only on the biased controller strategy which has consistently yielded the best performance. Likewise, we do not compare it with the methods in [21, 22] given their poor performance on random time varying graphs. We again averaged the performance over $5000$ different RES realizations.

Bandwidth and spectrum of the desired state. In this experiment, we analyzed the impact that the desired state bandwidth and its GFT have on the controllability performance. We considered four different GFTs for the desired state, namely: (i) a step low-pass $[{\tilde{\mathbf{x}}}_{K}^{\ast}]_{k}=1$ for $k=1,\ldots,K$; (ii) a step high-pass, where the active frequencies correspond to the $K$ eigenvectors with highest total variation; (iii) a linear decay response given by $[{\tilde{\mathbf{x}}}_{K}^{\ast}]_{k}=1-(k-1)/K$ for $k=1,\ldots,K$; and (iv) an exponential decay response with $[{\tilde{\mathbf{x}}}_{K}^{\ast}]_{k}=e^{1-k}$ for $k=1,\ldots,K$. For a fair comparison, we normalized all desired states to unit energy and analyzed different values of $K$; a fraction of $N$ between $0.15$ and $0.27$.

The results are depicted in Figure 6. First, we observe that controlling the system to a higher bandwidth state is harder since the set of graph frequencies to guarantee controllability increases. Second, we observe that the high-pass response is harder to achieve and responses that decay to zero (like the linear decay and the exponential decay) yield lower MSE. This is because high-pass responses are translated in the vertex domain as states having dissimilar values in adjacent nodes. They render the control of network dynamics to such states more challenging. Hence, we conclude that it is easier to drive the network state to a signal that varies smoothly on the nodes compared with a signal that has highly different values in connecting vertices; e.g., it is easier to convince someone to vote a conservative candidate, if she is surrounded by members that have the same political inclination.

Sampling heuristics. In the last experiment, we focused on the impact of the control nodes. We compared the proposed constrained greedy selecting heuristic in Algorithm 1 with the optimal combinatorial solution and a uniformly random sampling scheme. We fixed $K=10$ and considered the linear decay desired state ${\tilde{\mathbf{x}}}_{K}^{\ast}$ such that $[{\tilde{\mathbf{x}}}_{K}^{\ast}]_{k}=1-(k-1)/K$ for $k=1,\ldots,K$. The obtained results are shown in Figure 7. We observe that the MSE decreases as more control nodes are selected. We also observe that the greedy heuristic yields a performance similar to the optimal solution and represents a considerable improvement over random selection.

VII Conclusions

In this paper, we studied the problem of controlling network states. We considered a random time varying network to be driven to a desired bandlimited state. To cope with the randomness in the underlying support, we introduced the concept of controllability in the mean, where we postulated to control the system as if it were running on the expected graph. We then carried out a detailed mean squared error analysis to quantify the deviation of the target signal, when the control is designed for the expected graph but ran on any given random network realization. We used this analysis to propose two different control strategies and evaluated their performance on both synthetic graph models and real-world social networks. We concluded that it is of paramount importance to take into account the random nature of the underlying topology. We leave as future work the analysis of more complex random network models, other parsimonious graph signal models, and other control strategies that involve spectral or energetic constraints. Another direction worth investigating is the proposal of other heuristic solutions to the respective optimization problems.

Appendix A Proof of Lemma 1.

Proof.

For model $(i)$, ${\mathbf{S}}={\mathbf{L}}$ and the system transition matrix is ${\mathbf{A}}_{t}={\mathbf{I}}-\epsilon{\mathbf{L}}_{t}$ for $0<\epsilon\leq 1/\|{\mathbf{L}}\|_{2}$. First, we prove that $\|{\mathbf{L}}_{t}\|_{2}\leq\|{\mathbf{L}}\|_{2}\leq\varrho$. Note that ${\mathcal{G}}_{t}\subseteq{\mathcal{G}}$ for every $t$ and, therefore, from the Laplacian interlacing property [43] this condition always holds. The proof of Assumption 2 is straightforward, i.e., from ${\mathbb{E}}[{\mathbf{A}}_{t}]={\mathbf{I}}-\epsilon{\mathbb{E}}[{\mathbf{L}}_{t}]={\mathbf{I}}-\epsilon p{\mathbf{L}}$, which means that ${\mathbb{E}}[{\mathbf{A}}_{t}]$ and ${\mathbf{L}}$ share the same eigenvectors. For the last condition, note that $\|{\mathbf{A}}_{t}\|_{2}=\|{\mathbf{I}}-\epsilon{\mathbf{L}}_{t}\|_{2}\leq 1$ since $\epsilon\leq 1/\|{\mathbf{L}}\|_{2}\leq 1/\|{\mathbf{L}}_{t}\|_{2}$. Therefore, $\|{\mathbf{A}}_{t}\|_{2}$ is upper bounded by some finite $\varrho$.

For model $(ii)$, ${\mathbf{S}}={\mathbf{W}}$ and the system transition matrix is ${\mathbf{A}}_{t}={\mathbf{W}}_{t}$. To prove that $\|{\mathbf{W}}_{t}\|_{2}\leq\|{\mathbf{W}}\|_{2}$, recall that for connected graphs, the largest eigenvalue is positive and real [44, Theorem 0.2]. Then, since ${\mathbf{W}}$ is considered to be normal and Assumption 1 holds, $\|{\mathbf{W}}\|_{2}=\lambda_{\max}({\mathbf{W}})\leq\max\deg({\mathcal{G}})\leq\varrho<\infty$. Likewise, since ${\mathcal{G}}_{t}\subseteq{\mathcal{G}}$, then $\max\deg({\mathcal{G}}_{t})\leq\max\deg({\mathcal{G}})<\infty$ and therefore $\|{\mathbf{W}}_{t}\|_{2}\leq\varrho<\infty$ for all $t$. The proofs of the last two conditions are straightforward since ${\mathbb{E}}[{\mathbf{A}}_{t}]=p{\mathbf{W}}$ and $\|{\mathbf{A}}_{t}\|_{2}=\|{\mathbf{W}}_{t}\|_{2}\leq\varrho<\infty$. This completes the proof. ∎

Appendix B Proof of Proposition 1 and Corollary 1

Proof of Proposition 1.

Recall that ${\mathcal{S}}$ is the set of the selected $M$ nodes and that ${\mathbf{C}}^{\mathsf{T}}{\mathbf{C}}=\text{diag}({\mathbf{c}})$, where ${\mathbf{c}}\in\{0,1\}^{N}$ with $[{\mathbf{c}}]_{i}=1$ if $v_{i}\in{\mathcal{S}}$ and $[{\mathbf{c}}]_{i}=0$, otherwise. System (10) is equivalent to

[TABLE]

where ${\hat{\mathbf{u}}}_{t}\in{\mathbb{R}}^{N\times 1}$ denotes the zero-extended control signal such that $[{\hat{\mathbf{u}}}_{t}]_{i}=[{\mathbf{u}}_{t}]_{i}$ if $v_{i}\in{\mathcal{S}}$ and $[{\hat{\mathbf{u}}}_{t}]_{i}=0$, otherwise. Then, system (26) is controllable iff the $K\times TN$ matrix

[TABLE]

is full rank. Observe that

[TABLE]

holds from $\text{rank}({\mathbf{A}}{\mathbf{B}})\leq\min\{\text{rank}({\mathbf{A}}),\text{rank}({\mathbf{B}})\}$. Therefore, to ensure the full rank $K$ of $\tilde{\boldsymbol{\Omega}}$, $M\geq\lceil K/T\rceil$ must hold, for some $T\geq 1$. This concludes the proof. ∎

Proof of Corollary 1.

Recall that, for two matrices ${\mathbf{X}}\in{\mathbb{R}}^{M\times N}$ and ${\mathbf{Y}}\in{\mathbb{R}}^{N\times K}$ [45, Section 0.4.6]

[TABLE]

The mean system (14) is controllable, iff

[TABLE]

has rank $K$ with ${\tilde{\mathbf{A}}}_{K}=\text{diag}({\bar{\mathbf{a}}}_{K})$.

The first term in (30) has rank

[TABLE]

while the second term has rank

[TABLE]

since ${\mathbf{Y}}$ consists of the Kronecker product of ${\mathbf{V}}_{K}^{\mathsf{H}}\text{diag}({\mathbf{c}})$ with an identity matrix. Note that ${\mathbf{V}}_{K}^{\mathsf{H}}\text{diag}({\mathbf{c}})$ selects indeed rows of ${\mathbf{V}}_{K}$.

Now, if $M\geq K$ and the node set ${\mathcal{S}}$ are such that the selected $M$ rows of ${\mathbf{V}}_{K}$ form a set of $K$ linearly independent vectors, then $\text{rank}({\mathbf{V}}_{K}^{\mathsf{H}}\text{diag}({\mathbf{c}}))=K$. This implies that $\text{rank}({\mathbf{Y}})=TK$ and in virtue of (29), we obtain

[TABLE]

yielding that the mean system (14) is controllable. ∎

Appendix C Proof of Theorem 1

Proof.

The MSE can be rewritten as

[TABLE]

where each term is computed next.

First, to compute ${\mathbb{E}}[{\mathbf{x}}_{T}]$ and ${\mathbb{E}}[{\mathbf{x}}_{T}^{\mathsf{T}}{\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}{\mathbf{x}}_{t}]$, note that ${\mathbf{x}}_{T}$ can be written as

[TABLE]

where $\boldsymbol{\Phi}_{b,a}={\mathbf{A}}_{b}{\mathbf{A}}_{b-1}\cdots{\mathbf{A}}_{a+1}{\mathbf{A}}_{a}$ is the state transition matrix in the interval $[a,b]$ for $b>a$. Since under the $\text{RES}(p)$ model ${\mathbf{A}}_{t}$ are i.i.d. matrices, ${\mathbb{E}}[\boldsymbol{\Phi}_{b,a}]={\bar{\mathbf{A}}}^{b-a+1}$. Thus, the expectation of (35) is

[TABLE]

For the second order moment ${\mathbb{E}}[{\mathbf{x}}_{T}^{\mathsf{T}}{\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}{\mathbf{x}}_{T}]$, denote by ${\mathbf{Q}}={\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}$ and by substituting (35) we have

[TABLE]

Define $\boldsymbol{\Gamma}_{\tau,\tau^{\prime}}={\mathbb{E}}[\boldsymbol{\Phi}_{t-1,\tau+1}{\mathbf{Q}}\boldsymbol{\Phi}_{t-1,\tau^{\prime}+1}]\in{\mathbb{R}}^{N\times N}$ [cf. (15)], so that (37) can be compactly written as

[TABLE]

Finally, by substituting (38) in the first term of the MSE (34) and (36) in the second term, we obtain the claimed expressions. This completes the proof.∎

Appendix D Special Cases: Useful Computations of the Quadratic Term $\Gamma_{\tau,\tau^{\prime}}$ in Theorem 1

Computation of the quadratic term $\Gamma_{\tau,\tau^{\prime}}$ in Theorem 1 can turn out to be quite cumbersome for arbitrary graph shift operators ${\mathbf{S}}_{t}$ or transition matrices ${\mathbf{A}}_{t}$. In what follows, we offer two corollaries of Theorem 1 that address this issue. In particular, Corollary 2 gives an upper bound on the MSE that does not entail computation of second-order moments, while in Corollary 3 we show that, for the usually found case of undirected graphs, diffusion models in Lemma 1 admit an exact computation. Proofs follow after the statement of the corollaries.

Corollary 2**.**

Under the same conditions of Theorem 1 and from Lemma 1, the MSE (15) can be upper bounded by

[TABLE]

The result of Corollary 2 can be interpreted as the worst case scenario to account for the variability in the topology. In fact, (39) shows only first order dependence from the $\text{RES}(p)$ model, but it does not show dependence from the second order moment. For the models in Lemma 1, we can consider for model $(i)$ $\varrho=1$ (${\bar{\mathbf{A}}}={\mathbf{I}}-\epsilon p{\mathbf{L}}$), while for model $(ii)$ $\varrho=\max\deg({\mathcal{G}})$ (${\bar{\mathbf{A}}}=p{\mathbf{W}}$). Further insight then on the role of the graph variability is given by Corollary 3, which shows the explicit dependence on the link activation probability $p$.

Corollary 3**.**

Under the same conditions of Theorem 1 and additionally given that ${\mathcal{G}}$ is an undirected graph, for the diffusion models in Lemma 1, the following holds:

[TABLE]

for $\tau\leq\tau^{\prime}$, ${\mathbf{Q}}_{0}={\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}$, and ${\bar{\mathbf{A}}}^{\tau^{\prime}-\tau}=({\mathbb{E}}[{\mathbf{A}}_{t}])^{\tau^{\prime}-\tau}$.

For model $(i)$ in Lemma 1 (${\bar{\mathbf{A}}}={\mathbf{I}}-\epsilon p{\mathbf{L}}$), ${\mathbf{Q}}_{a}$ is

[TABLE]

For model $(ii)$ in Lemma 1 (${\bar{\mathbf{A}}}=p{\mathbf{W}}$), ${\mathbf{Q}}_{a}$ is

[TABLE]

Proof of Corollary 2.

From the first term in (34), we have

[TABLE]

The trace argument in (43) can be expanded as

[TABLE]

From Lemma 1, we have $\|{\mathbf{A}}_{t}\|_{2}\leq\varrho$; so that, by the submultiplicativity of the spectral norm, we can write

[TABLE]

Also, observe that for any square matrix ${\mathbf{X}}$ and positive semidefinite matrix ${\mathbf{Y}}$, it holds that $\text{tr}[{\mathbf{X}}{\mathbf{Y}}]\leq\|{\mathbf{X}}\|_{2}\text{tr}[{\mathbf{Y}}]$ [46].

Given that the filter ${\mathbf{H}}$ does not amplify any frequency (i.e. $\|{\mathbf{H}}\|_{2}=1$), then plugging back (44) and (45) into (43) yields

[TABLE]

Using (46) to bound (43) and replacing it in (15) yields (39). ∎

Proof of Corollary 3.

Consider $\tau\leq\tau^{\prime}$. From (15), we have

[TABLE]

Then, since for two random matrices ${\mathbf{X}},{\mathbf{Y}}$ it holds that ${\mathbb{E}}[{\mathbf{X}}]={\mathbb{E}}[{\mathbb{E}}[{\mathbf{X}}|{\mathbf{Y}}]]$ [47, Theorem 34.4], (D) becomes

[TABLE]

which under the $\text{RES}(p)$ model (i.e., matrices ${\mathbf{A}}_{a}$ are i.i.d.) can be written as

[TABLE]

Further, for $a\geq 1$ and assuming for now (to be proven later on) that ${\mathbf{Q}}_{a-1}$ is symmetric and positive semidefinite, we proceed to compute the $(i,j)$ entry of matrix ${\mathbf{Q}}_{a}={\mathbb{E}}[{\mathbf{A}}_{T-a}^{\mathsf{T}}{\mathbf{Q}}_{a-1}{\mathbf{A}}_{T-a}]$. Towards this end, denote simply by $[{\mathbf{Q}}_{a-1}]_{ij}=q_{ij}$ and $[{\mathbf{A}}_{T-a}]_{ij}=a_{ij}$ for $i,j=1,\ldots,N$.

For $i\neq j$ and since ${\mathbf{A}}_{T-a}$ is symmetric, the $(i,j)$ element of ${\mathbf{Q}}_{a}$ becomes

[TABLE]

where we have used the independence of the distinct elements in ${\mathbf{A}}_{T-a}$. The second term of (50) groups the element $(i,j)$ together with $(j,i)$ due to symmetry of ${\mathbf{A}}_{T-a}$. Analogously, for the diagonal elements $i=j$, we get

[TABLE]

With this in place, let us fist consider the simpler model $(ii)$ in Lemma 1, where ${\mathbf{S}}_{t}={\mathbf{W}}_{t}$ and ${\mathbf{A}}_{t}={\mathbf{W}}_{t}$. For this case, we have $[{\mathbf{A}}_{t}]_{ij}=Bw_{ij}$, where $B$ is a Bernoulli random variable of parameter $p$ and $w_{ij}=[{\mathbf{W}}]_{ij}$. Then, by substituting ${\mathbb{E}}[a_{ij}]=pw_{ij}$ and ${\mathbb{E}}[a_{ij}^{2}]=(p^{2}+p(1-p))w_{ij}^{2}$ in (50), we get

[TABLE]

which can be written in the compact form

[TABLE]

Likewise, for $i=j$ (51) becomes

[TABLE]

which can also be written compactly as

[TABLE]

By combining (53) and (55) yields (42). Finally, note that if ${\mathbf{Q}}_{a-1}$ is symmetric and positive semidefinite, then so is ${\mathbf{Q}}_{a}$. To complete the proof, observe that ${\mathbf{Q}}_{0}={\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}$ is symmetric and positive semidefinite, thus (42) holds for all $a\geq 1$.

For model $(i)$ in Lemma 1, we proceed in an analogous way. In this case, ${\mathbf{S}}_{t}={\mathbf{L}}_{t}$ and ${\mathbf{A}}_{t}={\mathbf{I}}-\epsilon{\mathbf{L}}_{t}=({\mathbf{I}}-\epsilon{\mathbf{D}}_{t})+\epsilon{\mathbf{W}}_{t}$, where ${\mathbf{D}}_{t}=\text{diag}({\mathbf{W}}_{t}{\mathbf{1}})$ is the degree matrix. This means that $[{\mathbf{A}}_{t}]_{ij}=a_{ij}=\epsilon Bw_{ij}$ if $i\neq j$ and $[{\mathbf{A}}_{t}]_{ii}=a_{ii}=1-\epsilon\sum_{k=1}^{N}B_{k}w_{ik}$ with $B_{k}$ being i.i.d. Bernoulli random variables with probability $p$.

Then, for $i\neq j$, from (50) we have

[TABLE]

Now, recalling that $w_{ii}=0$ (i.e., no self-loops), (56) becomes

[TABLE]

which can be further written in the compact form

[TABLE]

For $i=j$, we start with (51)

[TABLE]

Then, recalling that $w_{ii}=0$, we replace the first and second order moments for each $a_{ij}$ and obtain

[TABLE]

Finally, this can be rewritten as

[TABLE]

completing the proof. ∎