{'fraction_non_alphanumeric': 0.05774141681641918, 'fraction_numerical': 0.01853778492386267, 'mean_word_length': 3.909914093336429, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Networked Linear Dynamical Systems under Non-white Excitation from a Single Trajectory

We consider a networked linear dynamical system with p agents/nodes. We study the problem of learning the underlying graph of interactions/dependencies from observations of the nodal trajectories over a time-interval T . We present a regularized non-casual consistent estimator for this problem and analyze its sample complexity over two regimes: (a) where the interval T consists of n i.i.d. observation windows of length T/n (restart and record), and (b) where T is one continuous observation window (consecutive). Using the theory of M -estimators, we show that the estimator recovers the underlying interactions, in either regime, in a time-interval that is logarithmic in the system size p. To the best of our knowledge, this is the first work to analyze the sample complexity of learning linear dynamical systems driven by unobserved not-white wide-sense stationary (WSS) inputs.


Introduction
A networked linear dynamical system (LDS) is a network of agents/nodes, each of whose state evolves over time (in discrete or continuous steps) as a linear function of an external excitation and the states of its neighboring nodes in the network. The framework of LDS has been used to model dynamics in systems biology [1,2], financial markets [3], energy [4], transportation [5] and other critical networks [6,7]. Learning the dependencies, or topology learning, in a networked LDS is crucial for inference of influence pathways and subsequent control for the corresponding networks. As such, strategies to recover the underlying network structure from nodal time-series in LDS have been researched and can be classified into two categories: active and passive. Active learning involves efficient manipulation or interventions of nodal dependencies and injecting exogenous inputs into the LDS to infer the edges in the network by identifying the resulting changes [8]. Passive methods, on the other hand, use historical or streaming time-series of nodal states to infer the underlying topology. Our work falls within the domain of passive structure estimation. Very few works discuss learning such systems but are limited to the asymptotic regime (infinite sample limit). Examples include [9,10,11].
Prior Work: Tractable passive topology learning in networked LDS and Vector Auto-Regressive processes (VAR) has been shown using the framework of l 1 -regularized regression ( [12,13] and references therein), where the focus is on extending the results from the static Lasso or Graphical Lasso [14,15,16] to one with correlated samples, by showing that properties such Restricted strong convexity hold. A similar approach for continuous time stochastic differential equation has been studied in [17]. A graphical model for VAR processes, without performance guarantees, has been proposed in [18]. Least squared regression based identification of unstable dynamical systems using a single trajectory has been studied in [19,20]. However to the best of our knowledge, all existing work with sample complexity analysis of the learning algorithm relies on the assumption that unobserved exogenous inputs to the system are i.i.d. or white Gaussian noise, or that the exogenous inputs are observed [21,22]. On learning networked LDS with temporally correlated inputs (wide-sense stationary and cyclo-stationary) using state observations (not inputs), [11] presents a consistent algorithm using non-causal regression (Wiener filter) for exact recovery. Similarly, [23] recovers the underlying topology in networked LDS using the framework of directed mutual information. However, these works do not provide for guarantees in the finite sample regime, aside from numerical examples. The overarching goal of this work is thus to provide a structure learning algorithm for networked LDS driven by temporally-correlated inputs, with guarantees on its performance at finite samples. We present a regularized Wiener filter estimator for this problem and use the theory of M -estimators to determine the observation window necessary to guarantee correct topology estimation.
The rest of the paper is organized as follows. In Section 2, we describe the mathematical model of networked linear dynamical system and our consistent learning algorithm. The main results are presented in Theorems 2.2 and 2.4. Section 3 contains results on M-estimators used in the proof of our theorems, with sketches of proofs in Section 4. Section 5 contains simulation results, and Section 6 summarizes the paper and includes potential extensions and generalizations.

Main Results
Consider a graph G = (V, E) of p + 1 nodes in set V = {1, ..., p + 1} and undirected edge set E ⊂ V × V . We denote the set of two-hop neighbors in G by set E M , where E M = {(i, j)|(ij) ∈ E or ∃k, s.t (ik), (jk) ∈ E}. Note that E M \ E is the set of 'strict' two-hop neighbors in the graph G, that do not form edges in E. Each node i ∈ V is associated with a real-valued scalar state variable {x i (k), k ∈ Z} that evolves in discrete time 2 according to the following linear dynamical equation: where, {e i (k), k ∈ Z}, is an exogenous input. While samples of x i (k) are correlated in time due to the system dynamics, prior work on guaranteed learning of networked LDS include only temporally uncorrelated or white excitations/inputs e i (k). In this work, we consider e(k) k∈Z = [e 1 (k)...e p+1 (k)] T to be a zero-mean Wide-Sense Stationary (WSS) Gaussian process, uncorrelated across nodes, i.e., ∀k 1 , k 2 , τ ∈ Z, The time series vector x(k) k∈Z = [x 1 (k)...x p+1 (k)] T ∈ R p+1 is thus a zero mean jointly Gaussian WSS processes. We can write the frequency-domain representation of Eq. 1 at frequency f ∈ [0, 2π) as, where, spectral density Φ x , at frequency f as Topology Learning: Consider n state trajectories of all the nodes in V for the graph G = (V, E) excited by unobserved WSS (temporally correlated) inputs, such that the r th state trajectory (x r ) has N samples. Let T = n × N be equal to the total observation window. For the r th state trajectory, define the Discrete Fourier Transform (DFT) 3 is where, r ∈ {1, · · · , n}, i , · · · , X n i ] T ∈ C n×p respectively. We assume that X and Y are column-normalized, that is, For any quantity β ∈ C, we use (β) and (β) to denote its real and imaginary components. We list the following result from [11] that enables consistent estimation of all edges in E (as described in Figure 1), using nodal state trajectories.
The proof of Lemma 2.1 (see [11] for details) follows by showing that . The result then follows from algebraic properties of Φ −1 x (inverse power spectral density) derived from Eq. 3. It is worth noting that, in the time-domain, Eq. 7 is equivalent to a non-causal regression of the time-series, termed as "Wiener filter" [9]. This is effectively a non-causal extension of the connection between the inverse covariance matrix and the neighborhood regression used in learning static Gaussian graphical models [15,16,24]. For the finite sample regime, we study the problem of estimating edgesÊ such that P[Ê = E] ≥ 1 − for any user-defined threshold ∈ (0, .5). Estimating Φ x and then inverting it requires significant amount of data in the high dimensional setting. Instead, we use a regularized version of Eq. 7 as our graph estimator.

Regularized Wiener Filter Estimator
We propose a Regularized Wiener Filter EstimatorŴ i for a node i ∈ V as follows: where, λ > 0 is the regularization parameter. As β ∈ C p , β 1 is equal to the 1, 2-group norm over [ (β) (β)]. For thresholds τ 1 , τ 2 , we construct setŝ In the remaining of the paper, we find sufficient conditions on n, N and λ and fix thresholds such that P[Ê = E] ≥ 1 − , for given ∈ (0, 0.5). We consider two settings: (i) Restart & Record: The n trajectories of length N each are independent, i.e., the system is reset or sufficient gap is provided before the next N time-intervals are recorded. Thus, {x r } N r=1 are i.i.d., and (ii) Consecutive: In the second and more realistic setting, we consider the n state trajectories to be consecutive, i.e., {x r } N r=1 correspond to N -length intervals from a larger trajectory of length n × N . See Figure. 2 for the two settings considered in this article.

Main Results
The error in estimating W i by Eq. 8 and consequently in topology learning arises due to the finite N in computing X by DFT, as well as the finite n in the empirical regression. For our analysis, we consider the following non-zero parameters of the LDS over graph G = (V, E).
Note that under persistently exciting inputs, Φ e is a positive definite matrix almost surely at all frequencies [9]. Further, G is a connected network. Hence, under standard well-posedness assumptions, (I − H) and P −1 in Eq. 3 in the main document are full-ranked and L ≥ 0. Using norm bounds for matrix products, L and U can be bounded in terms of maximum and minimum eigen-values of (I − H * )(I − H) and P −1 . C, δ relate to the rate of decay of temporal correlation in the system states. Higher values of C and δ −1 imply greater temporal correlation. d, the maximum degree due to edges in E M , is upper-bounded by the square of the maximum nodal degree in G = (V, E).
For the case where the trajectories are i.i.d., the following theorem bounds the error in estimating W i .
holds with a probability of at least 1 − 1 , where c, c are universal positive constants. U, L, C, δ, m i , d are defined in Eq. 10.
Using the above result, we prove our first result on the correctness of our thresholding procedure listed in Eq. 9 for i.i.d. trajectories. We recall that n is the number of realizations of time-series, each of length N.
construct an undirected edge setÊ M andÊ as per Eq. 9 with thresholds τ 1 = τ 2 = m, then E =Ê holds with a probability of at least 1− . U, L, C, δ, m, d are defined in Eq. 10.
The proofs of Theorems 2.1 and 2.2 are provided in Section 4. Next we consider the case where the n trajectories correspond to consecutive intervals from the same time-series of length T = n × N . The following holds for the error for the regularized estimator in Eq. 8.
holds with a probability of at least 1 − 1 , where c, c are universal positive constants. U, L, C, δ, m i , d are defined in Eq. 10.
This leads to our second result on correct topology recovery in the non-i.i.d. case.
construct an undirected edge setÊ M andÊ as per Eq. 9 with thresholds τ 1 = τ 2 = m, then E =Ê holds with a probability of at least 1 − . U, L, C, δ, m, d are defined in Eq. 10.
The proof of Theorems 2.3 and 2.4 are provided in Section 4. In the next section, we present our approach for proving these theorems based on the theory of M-estimators [25], albeit in the complex domain.

Analysis of Regularized Wiener Filter Estimator
Our regularized regression involves working with complex-valued random variables X r i , X r i defined in Eq. 5 for a node i ∈ V . Their probability distribution is as follows, Thus, X r i ∼ N (0,Φ i ) and X r i ∼ N (0,Φ i ). The following result bounds the difference betweenΦ x and Φ x (see Eq. 4) for a N -length trajectory, and is used in our analysis.
The proof is provided in the Appendix. The regularized Wiener filter estimator Eq. 8 belongs to a class of regularized M -estimators. Note that the regularizer ( . 1 ) in Eq. 8 satisfies decomposability property with respect to the following complex-valued subspaces: We follow the approach in [25] to bound the error∆ [25] states that only two conditions are sufficient to control the error ∆ 2 .
As shown in [25], Eq. 13 ensures that∆, defined in Eq. 12, belongs to the set Second condition (restricted eigenvalue property): The following proposition, similar to Theorem 1 in [25], bounds the error ∆ 2 .
For completion, we outline a proof of Eq. 14 and Proposition 1 for complex-valued variables in Appendix 7.2, following the real-valued analysis in [25]. We now show that Eq. 13 The proofs for Lemmas 3.2 and 3.3 are provided in Section 4.1.
Consecutive (non i.i.d.) trajectories: For the non i.i.d. scenario, the guarantee on lambda and restricted eigenvalue condition are stated next with proofs provided in Section 4.2.
, with a probability of at least 1 − 2 .

Proof of Main Theorems
2 holds with a probability of at least 1 − 1 .
Proof of Theorem 2.2. Choose 1 = p . It follows from definition of m, that 1 d would satisfy the conditions on n and λ specified in Theorem 2.1 for a i ∈ V . Therefore, Ŵ i − W i 2 ≤ m 2 holds with a probability of at least 1 − p . Using a union bound for all the p + 1 nodes, we have Next we prove the lemmas on λ and restricted eigenvalue property, introduced in Section 3, that are used in the Theorems' proofs.

Proofs of M -estimator conditions for (restart & record) trajectories
Proof of Lemma 3.2. Let E := Y − X W i . We show that 1 n X H E ∞ is bounded with a high probability and choose λ greater than that bound. Separating X = X R + ιX I , E = E R + ιE I , into real and imaginary parts (specified by subscripts R and I respectively), we have Let some symmetric permutation matrix P , such that its covariance matrix C 2 = P C 1 P . Rewriting To bound the right side of Eq. 17, we first show that either function is Lipschitz. Consider first Thus, f is a Lipschitz function with Lipschitz constant ) nL and the union bound for all j ∈ {1, · · · , p}, Using a similar analysis, Using the Union bound on Eq. 17, we have Proof of Lemma 3.3. Note that X = [X 1 i , · · · , X n i ] T , where X r i is computed from r th trajectory, defined in Eq. 5. Separating into real and imaginary parts, we have: where X 1 := [X R − X I ], X 2 := [X I X R ] and v = (∆ T R ∆ T I ) T . For simplicity, in this proof we drop the superscript r in X r i and X r i . Note that the rows of X 1 and X 2 are i.i.d. samples of the real random vectors, holds for all ∆ ∈ D(W i ) with high probability, we prove the restricted eigenvalue property for group structured norms on both terms in Eq. 19. LetΣ be the covariance of random vector [(X i ) T R (X i ) T I ] T . Then, Thus, [(X i ) T R − (X i ) T I ] T and [(X i ) T I (X i ) T R ] T have means 0 and covariance Σ 1 and Σ 2 , respectively, where, From Eq. 19, Consider the following definitions: Clearly, if ∆ ∈ D(W i ), defined in Eq. 14, then v ∈ D 2 (W i ) and vice versa. For v ∈ D 2 , the group norm v 1,2 ≤ 4 v M2 1,2 . Using Cauchy Schwartz inequality and the definition of bounded degree Below is a result, derived from Section 5 of [25] and [27] for Gaussian random matrices.
We use Eq. 22 with X = X 1 and X 2 to obtain a lower bound on Note that, forΣ, Σ 1 and Σ 2 defined in Eqs. 20, 21, we have ρ for two non-negative numbers a, b in Eq. 19, we get that the following holds with a probability of 1 − 2c exp(−cn), where c, c are universal positive constants: Using Eq. 21, the definition ofΦ i : Thus, λ min ((Σ 1 . Following Lemma gives the result.
Proof of Lemma 3.5. Let Z := X ∆ ∈ C n . Its real and imaginary components are We find the lower bounds on 1 n Z R 2 2 , 1 n Z I 2 2 and then combine them to obtain a lower bound of 1 n Z 2 2 . Applying Lemma I.2 from the Supplementary material of [28] on Z R and Z I , we have, with individual probability atleast 1 − [2 exp(− (25) Note that the diagonal values of E(Z R Z T R ) are all equal to v T Σ 1 v, and those of E(Z I Z T I ) are equal to v T Σ 2 v with Σ 1 , Σ 2 defined in Eq. 21. Using this with Lemma 7.1 in the Appendix, we get, holds with probability of atleast 1 − [4 exp(− n for a large n). Then for n ≥ 33 2 log p[ U L + 0.5 + 4 , the statement holds whenever n ≥ 2 log( 4p 2 p 2 2 −4 ).

Numerical Results
We demonstrate the numerical implementation of recovering topology on a Desktop PC with Intel Xeon E5-1620 Processor (8x 3.7 GHz) and 32 GB RAM. We considered a two-dimensional square grid G with p+1 nodes and generate samples for Eq. 1 with exogenous input e(k) = s p 5[w(k)−0.3w(k− 1)], where w(k), w(k − 1) are sampled from a standard Normal distribution, s p ∈ R (p+1)×(p+1) is a diagonal matrix containing constants; h ∈ R (p+1)×(p+1) is a weighted adjacency matrix. The scaling s p is chosen such that X and Y in the estimator are column-normalized.
We reconstruct the topology with a probability of at least 1 − , where = 0.05. The numerical experiments are conducted in MATLAB R2020b.
For a choice of n, we generate n trajectories, either independently (restart & record) or taken as consecutive intervals of a larger trajectory (consecutive). Each trajectory is of length N = 4CU δ −1 (1−δ −1 ) 2 , rounded to the nearest integer. From the trajectories, we compute the samples {X r i , X r i } n r=1 for all the nodes i ∈ V at frequency f = π 2 . The Regularized Wiener Filter Estimator is solved using CVXR [29], with λ  Figure. 3. In Figure 3(b), the correlation strength δ −1 of the trajectories is high, and consequently N is large. For small δ −1 , the length of each trajectory can be reduced significantly for reconstructing the topology. For example, in Figure 3(c), N is much smaller. Further, n min is of the order of ≈ 10 7 rather than a more conservative estimate of ≈ 10 16 as provided by the main theorems.

Extensions and Path Forward
In this article, we presented a regularized Wiener filter estimator to learn the structure of a discretetime networked LDS. We analyzed the sample complexity of our estimator and showed that it linearly depends the logarithm of the number of nodes p in two cases, one where trajectories of nodal states are collected in independent observation windows of equal length, and another where the trajectories pertain to a single continuous observation window.
While we discuss our method for first-order discrete-time LDS, our estimator can be extended to learning related networks as highlighted next.
VAR(p) models with correlated inputs: Lemma 2.1 and our subsequent analysis follows directly if higher-order delays (at the same node) are included in the LDS Eq. 1. L, U and the sample complexity will need to be changed accordingly.

Continuous time LDS:
Considering a fixed sampling time ∆T and a time-discretization function, the continuous time LDS can be converted to a discrete-time LDS with related frequency domainrepresentation [11]. The analysis will involve merging the error due to discretization with the finite sample analysis.
LDS under cyclo-stationary processes: Cyclo-stationary processes represent a generalization of WSS processes where the statistics such as mean, correlation function are periodic functions of time. As shown in [30], a lifting operation can be used to represent time-evolution of a cyclo-stationary process as a WSS process with vector-valued states. The remaining analysis of the sample complexity will be similar but involve vector-valued states.
Directed graphs under correlated inputs: Note that our estimator uses properties of the inverse power spectral density (see Lemma 2.1 and discussion). Under a strict causality assumption on the linear filters [9,23], it has been shown that directed edges can be recovered using inverse power spectral density. The framework presented here can thus be extended to efficiently learn the topology for strictly causal linear filters in the LDS setting.
Finally, we plan to analyze the restrictions of our algorithm in learning networked LDS with spatially correlated inputs, and estimating directed networks with non-causal dependencies, where only approximate reconstruction may be possible using passive methods.
Using the same approach,

M-estimator in complex-domain
Our estimator (Eq. 8) is a M-estimator in the complex domain. For completion, we present the proof of Eq. 14 and Proposition 1 by following the real-valued analysis in [25].
Eq. 14 states that if λ ≥ 2 To prove that, note that where, Re(x) denotes the real part of the complex number x. Moreover, . By optimality ofŴ i (λ) = W i +∆ in the Regularized Wiener Filter Estimator, , whenever Eq. 13 and Eq. 15 hold.
Proof of Proposition 1. Let K(δ) := {∆ ∈ C p | ∆ ∈ D(W f i ) and ∆ 2 = δ}. Let F(∆) be the difference between the objective of the Regularized Wiener Filter Estimator evaluated at W i + ∆ and W i . For a ∆ ∈ K(δ), the following holds: Note that F(0) = 0. Using Lemma 4 from the Supplementary material of [25] (uses convexity of F(∆)), it then follows that

Results on Covariance Matrices
Here we present few results on norms involving state covariance and error matrices that are used in the proofs in the main document.
Proof. For r ∈ {1, · · · , n}, X r i (r th entry of Y) and (X r i ) T (r th row of X ) are computed using N consecutive samples from {x((r − 1)N ), · · · , x((r − 1)N + N − 1)}, using Eq. 5. For r, c ∈ {1, · · · , n} we have, Consider In the remaining, we find an upper bound for max r c,c =r B rc . From Eq. 30, Thus, max r n c=1,c =r ≤ C(1 + Sa) 2δ −N 1 − δ −N (∵ δ −N − 1 is negative and can be ignored), Proof. As the n trajectories are i.i.d., using the Schur complement lemma, we get, From the definition of L, U in Eq. 10, we have, 1 Proof. Writing real and imaginary parts of E 1 and using inequality of matrix 2 and ∞-norms, we have,