{'fraction_non_alphanumeric': 0.06506493219176566, 'fraction_numerical': 0.025663149045375103, 'mean_word_length': 3.8793177892918824, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Sparse High-Dimensional Graphical Model Learning For Dependent Time Series

We consider the problem of inferring the conditional independence graph (CIG) of a sparse, high-dimensional stationary multivariate Gaussian time series. A sparse-group lasso-based frequency-domain formulation of the problem based on frequency-domain sufficient statistic for the observed time series is presented. We investigate an alternating direction method of multipliers (ADMM) approach for optimization of the sparse-group lasso penalized log-likelihood. We provide sufficient conditions for convergence in the Frobenius norm of the inverse PSD estimators to the true value, jointly across all frequencies, where the number of frequencies are allowed to increase with sample size. This results also yields a rate of convergence. We also empirically investigate selection of the tuning parameters based on Bayesian information criterion, and illustrate our approach using numerical examples utilizing both synthetic and real data.

central role in our proposed approach, and use of DFT can also be viewed as a way to decorrelate time-domain dependent data by transforming it to frequency-domain where {d x (f m )} n/2 m=0 is an approximately independent sequence (see Sec. II-A). Representing data {x(t)} n−1 t=0 as well as the DFT sequence {d x (f m )} n−1 m=0 as column vectors, the two can be related via a unitary n × n matrix, signifying that {d x (f m )} n−1 m=0 represents any {x(t)} n−1 t=0 via an orthonormal basis in C n [24, p. 71]. One may also view this transformation (linear projection) as feature extraction from raw data, more suitable than raw data for further processing for some intended task. In our case this transformation is invertible. In other applications, not necessarily related to graphical modeling, one may resort to combined projection and dimension reduction to a lower rank subspace; possible examples include [25], [26] and references therein. Low rank bilinear projections of matrices (or higherorder tensors) are proposed in [25] to improve the performance of existing two-dimensional discriminant analysis algorithms for classification tasks. A self-supervised learning method is proposed in [26] to train a deep feature extraction network without the need for numerous manually labeled video samples. It is not yet clear how such ideas and approaches apply to graphical model learning for dependent time series.

A. Related Work
Prior work on graphical modeling for dependent time series in the low-dimensional settings (sample size n p) is concerned with testing whether {i, j} ∈ E for all possible edges in the graph, based on some nonparametric frequencydomain test statistic such as partial coherence [10], [11], [17]- [19], [27], [28] which, in turn, is based on estimates of S x (f ) given time-domain data. These approaches do not scale to high dimensions where p is comparable to or larger than the sample size n. As an alternative to nonparametric modeling of time series, parametric graphical models utilizing a (Gaussian) vector autoregressive (VAR) process model of x(t), have been advocated in [13], [29], [30], but these approaches are suitable for only small values of p. As an alternative to exhaustive search over various edges, a penalized maximum likelihood approach in conjunction with VAR models has been used in [14] where the penalty term incorporates sparsity constraints, making it suitable for high-dimensional setting. For every pair of series components, the corresponding partial coherence is thresholded to decide if it is zero (exclude the edge), or nonzero (include the edge). No systematic (or principled) method is given for threshold selection.
Nonparametric approaches for graphical modeling of real time series in high-dimensional settings have been formulated in frequency-domain in [31], [32] using a neighborhood regression scheme, and in the form of penalized loglikelihood in frequency-domain in [33]- [35], all based on estimates of S x (f ) at various frequencies, given time-domain data. A key model assumption used in these papers is that 1], is sparse, i.e., for any f , the number of off-diagonal nonzero elements are much smaller than the total number p(p−1) of off-diagonal elements. Sparsity is enforced during optimization of the chosen objective function, making the problem well-conditioned. Reference [35] considers latent variable graphical models. References [34], [35] exploit the framework of [33], and therefore, inherit some of its drawbacks, discussed in more detail later in Sec. IV-B (see also Remark 1 in Sec. II-E). Reference [34] considers a Bayesian framework with a focus on multiple time series realizations whereas [33] deals with a single realization of the time series. Sufficient conditions for graph edge recovery are provided in [33] whereas there is no such analysis in [34], [35]. A recent approach [36] considers a vector Gaussian time series with uncorrelated samples but possibly nonstationary covariance matrices. This paper [36] addresses real-valued time series using a neighborhood regression scheme. However, the results in [36] are applied to the discrete Fourier transform (DFT) of stationary time series without analyzing the resultant complexvalues series (complex in the frequency-domain), and without noting the fact the DFT over the entire observation set does not result in (approximately) uncorrelated Gaussian sequence; see Secs. II-A and IV-B later where it is pointed out that only "half" of the DFT sequence is uncorrelated.

B. Our Contributions
In this paper we address the same problem as in [33], namely, first estimate the inverse PSD S −1 x (f ) at distinct frequencies, given time-domain data, and then select the graph edge {i, j} based on whether or not [S −1 x (f )] ij = 0 for every f . As in [33], we use a penalized log-likelihood function in frequency-domain as our objective function to estimate S −1 x (f ), given time-domain data. However, there are significant differences between our choice of frequencies and objective function, and our analysis, and that of [33]. We enumerate these differences below and elaborate on them some more later in the paper after introducing some concepts and notation in Sec. II.
(i) Our log-likelihood function is different from that in [33]. In terms of the DFT d x (f m ), f m = m/n, of real-valued time-domain data {x(t)} n−1 t=0 , defined in (1), we establish that {d x (f m )} n/2 m=0 , (n even), is a sufficient statistic for our problem, a fact not recognized in [33] (see also Remark 1 in Sec. II-E, and Sec. IV-B), who uses some redundant frequencies f m 's from the set ( 1 2 , 1), in addition to those from the set [0, 1 2 ]. Our proposed frequencydomain formulation is based on {d x (f m )} (n/2)−1 m=1 , neglecting d x (f m ) at m = 0 and m = n/2 where the DFT is real-valued Gaussian. We use a sparse-group lasso penalty on off-diagonal elements of inverse PSD estimates at individual frequencies (lasso) as well as jointly across all frequencies (group lasso), whereas [33] uses only a group lasso penalty that is applied to all elements, including diagonal elements, of inverse PSD estimates. The objective of lasso/group lasso penalty is to drive possible zero entries of inverse PSD matrix to zero. However, since PSD matrix at any frequency is assumed to be positive-definite, its inverse is also positive-definite, hence, all diagonal elements of inverse PSD matrix at any frequency are strictly positive, therefore, do not need any sparsity penalty. Our sparse-group lasso penalty is more general than the group lasso penalty of [33] and it reduces to group lasso penalty if we remove the lasso penalty at individual frequencies. In our optimization approach, we allow for the possibility of zero lasso penalty. (ii) In [33], a test statistic based on the estimated inverse PSD's is compared to a nonzero threshold to infer presence/absence of edges in the graph where no method is given regarding how to choose the threshold. In this paper the threshold is set to zero. We also empirically investigate selection of the tuning parameters (lasso and group-lasso penalty parameters) needed for our sparsegroup lasso approach, based on the Bayesian information criterion (BIC). No such results are available in [33]. (iii) We provide sufficient conditions (Theorem 1 in Sec. V) for convergence in the Frobenius norm of the inverse PSD estimators to the true value, jointly across all frequencies, where the number of frequencies are allowed to increase with sample size so that the entire sufficient statistic set is exploited. This results also yields a rate of convergence of the inverse PSD estimators with sample size n. Reference [33] provides sufficient conditions for consistent graph edge recovery in [33,Proposition III.2] when the number of frequencies are fixed independent of sample size, and [34], [35] offer no theoretical analysis. A consequence of using fixed number of frequencies is that in [33], only a subset of sufficient statistic set is exploited; we elaborate on this later in Remark 1 in Sec. II-E. We do not have any counterpart to [33,Proposition III.2] whereas [33] does not have any counterpart to our Theorem 1. (iv) We propose an alternating direction method of multipliers (ADMM) approach for optimization of the sparse-group lasso-based penalized log-likelihood formulation of the problem. Note that [33] also used ADMM but only for group lasso formulation. We discuss practical ADMM issues such as choice of ADMM penalty parameter and stopping rule for termination of the algorithm whereas there is no such discussion in [33].

C. Relationship to Prior Conference Publications
Preliminary versions of parts of this paper appear in conference papers [37], [38]. The sufficient statistic for our problem was first discussed in [37] where a sparse-group lasso problem was also formulated in frequency-domain. But an alternating minimization (AM) based solution to this problem, using the penalty method, is given in [37]. The performance of this method depends upon the penalty parameter, and strictly speaking, convergence of this solution to the desired solution requires the penalty parameter to become large, which can make the problem numerically ill-conditioned. Use of the ADMM approach mitigates the dependence on the penalty parameter used in the AM approach. The theoretical analysis (Theorem 1) presented in Sec. V is partially given in [38] where the proof is incomplete. In this paper we provide a complete proof and also correct some typos. In turn, the proof of Theorem 1 relies on some prior results in [39] which deal with complex-valued Gaussian vectors, not time series. Lemma 1 appears in [37] without a proof. Lemma 3 is from [38] where a complete proof is given, here we simply state and use it. A version of Lemma 4 is stated without proof in [39], here we give a complete proof. Lots of details of proof of Theorem 1 are missing in [38]. The material in Secs. IV-A5, IV-C and VI of this paper does not appear in [37]- [39].

D. Outline and Notation
The rest of the paper is organized as follows. The sufficient statistic in the frequency-domain for our problem and the resulting log-likelihood formulation of the problem are presented in Sec. II where we also provide some background material on Wirtinger calculus needed for optimization w.r.t. complex variables, and Whittle likelihood used in [33]. A sparse-group lasso-based penalized log-likelihood formulation of the problem is introduced in Sec. III. An ADMM algorithm is presented in Sec. IV to optimize the objective function to estimate the inverse PSD and the edges in the graph. Selection of the tuning parameters based on BIC is presented in Sec. IV-C. In Sec. V we analyze consistency (Theorem 1) of the proposed approach. Numerical results based on synthetic as well as real data are presented in Sec. VI to illustrate the proposed approach. We present numerical results for both synthetic and real data. In synthetic data example the ground truth is known and this allows for assessment of the efficacy of various approaches. In real data example where the ground truth is unknown, our goal is visualization and exploration of the dependency structure underlying the data. In both cases we wish to compare our proposed approach with the widely used i.i.d. modeling approach where the underlying time series is either assumed to be i.i.d., or one uses only the covariance of the data. Proofs of Lemma 1 and Theorem 1 are given in Appendices A and B, respectively.
The superscripts * , and H denote the complex conjugate, transpose and Hermitian (conjugate transpose) operations, respectively, and the sets of real and complex numbers are denoted by R and C, respectively. Given A ∈ C p×p , we use φ min (A), φ max (A), |A|, tr(A) and etr(A) to denote the minimum eigenvalue, maximum eigenvalue, determinant, trace, and exponential of trace of A, respectively. We use A 0 and A 0 to denote that Hermitian A is positive semi-definite and positive definite, respectively. The Kronecker product of matrices A and B is denotes by A ⊗ B. For B ∈ C p×q , we define the operator norm, the Frobenius norm and the vectorized 1 norm, respectively, as and we also use θ for θ 2 . Given A ∈ C p×p , A + = diag(A) is a diagonal matrix with the same diagonal as A, and A − = A − A + is A with all its diagonal elements set to zero. We use A − * for (A * ) −1 , the inverse of complex conjugate of A, and A − for (A ) −1 . Given A ∈ C n×p , column vector vec(A) ∈ C np denotes the vectorization of A which stacks the columns of the matrix A. The notation y n = O P (x n ) for random y n , x n ∈ C p means that for any ε > 0, there exists 0 < M < ∞ such that P ( y n ≤ M x n ) ≥ 1 − ε ∀n ≥ 1. The notation x ∼ N c (m, Σ) denotes a complex random vector x that is circularly symmetric (proper), complex Gaussian with mean m and covariance Σ, and x ∼ N r (m, Σ) denotes realvalued Gaussian x with mean m and covariance Σ. The abbreviations PSD, w.r.t., w.h.p., h.o.t., iff and pdf stand for power spectral density, with respect to, with high probability, higher-order terms, if and only if, and probability density function, respectively.

II. BACKGROUND AND PROBLEM FORMULATION
Given time-domain data {x(t)} n−1 t=0 originating from a p−dimensional stationary Gaussian sequence, our objective is to first estimate the inverse PSD S −1 x (f ) at distinct frequencies, and then select the edge {i, j} in the time series GGM G based on whether or not [S −1 x (f )] ij = 0 for every f . We will follow a maximum likelihood approach, and to this end we need to express the likelihood function of time-domain data in terms its PSD matrix.

A. Sufficient Statistic and Log-Likelihood
, where f m = m/n, m = 0, 1, · · · , n − 1. Since is periodic in m with period n, and is periodic in normalized frequency f m with period 1. Since so d x (f m ) for m = 0, 1, · · · , (n/2), (n even), completely determines d x (f m ) for all integers m. As proved in [40, p. 280, Sec. 6.2], for any statistical inference problem, the complete sample is a sufficient statistic, and so is any one-to-one function of a sufficient statistic. Since the inverse DFT yields (one-to-one transformation) m=0 is a sufficient statistic, which can be further m=0 is a sufficient statistic for any statistical inference problem, including our problem of estimation of inverse PSD.

B. Complex Gaussian Vectors
Here we recall some facts regarding proper and improper complex Gaussian random vectors from [42]. We need these results to clarify different expressions for the pdf of a proper complex Gaussian vector, used later for optimization w.r.t. complex variables using Wirtinger calculus [42, Appendix 2], [43]. Define R vw = E{vw } for (zero-mean) v, w ∈ R p , and define the covariance matrix R vw = E{vw H }, and the complementary covariance matrixR vw = E{vw } [42, Sec. 2.2], for zero-mean v, w ∈ C p . Given u = u r + ju i ∈ C p , with real part u r and imaginary part u i , define the augmented complex vector y and the real vector z as The pdf of an improper complex Gaussian u is defined in terms of that of the augmented z or y [42, Sec. 2.3.1]. Assume E{u} = 0. Then we have z ∼ N r (0, R zz ) where Since z ∼ N r (0, R zz ), its pdf is given by (assuming R zz 0) One can also express (8) For proper complex u,R uu = 0, and (9) reduces to Since R uu = R H uu , |R uu | = |R * uu |, it follows that (u H R −1 uu u) * = u H R −1 uu u, for proper u, and therefore, (10) has the familiar form used in (4).

D. Whittle Likelihood
Define y = [x (0) x (1) · · · x (n − 1)] ∈ R pn . By assumption y ∼ N r (0, Σ y ) where Σ y = E{yy }. Since Σ y 0, the pdf of y is given by Based on some large sample (n → ∞) results of Whittle [46]- [48], Whittle approximation to f y (y) is stated in [ which specifies the joint pdf up to some constants, while [49, Eqn. (5)] specifies up to some constants. As noted in Sec. II-A, the terms in (17) and (18) corresponding to the indices m = n 2 +1 through n−1 are completely specified by terms corresponding to the indices . Therefore, unlike (4) and (5), (17) and (18), respectively, have lots of redundant frequencies. We note that [33] uses a likelihood function based on such Whittle approximation. The likelihood of [33] is examined further in Remark 1 in Sec. II-E.

E. Problem Formulation
Recall that our objective is to first estimate the inverse PSD S −1 x (f ) at distinct frequencies, and then select the edge {i, j} in the time series GGM G based on whether or not Then the problem is separable in m, and we choose Φ m to minimize where the expression after equality above follows by specify- (10) (correct way to handle complex variates [42]). Using Wirtinger calculus (Sec. II-C), at the optimal solution the gradient of h(Φ m ) w.r.t. Φ * m vanishes: leading to the solution (Φ m denotes estimate of Φ m ) Since [24] and it is known to be a poor estimator of the . Under high dimensional case, one also needs sparsity constraints in order to regularize the problem.
We assume that S x (f m ) is locally smooth (a standard assumption in PSD estimation [24,Chapter 2], [41]), so that S x (f m ) is (approximately) constant over K = 2m t + 1 consecutive frequency points. Pick leading to M equally spaced frequenciesf k in the interval (0, 0.5), at intervals of K/n. It is assumed that for eachf k (local smoothness), Using (24) in (4), we have Defině whereŜ k represents PSD estimator at frequencyf k using unweighted frequency-domain smoothing [41]. The joint pdf ofĎ(f k ) is given by Then we can rewrite (26 ) as This is the joint pdf we will use in the rest of the paper. Remark 1. Ref. [33] cites the Whittle approximation (18) given in [49] as a basis for their negative log likelihood function, which can be inferred from [33,Eqn. (6)] after removing the lasso penalty therein. In the notation of this paper, it is given by whereŜ k is the PSD matrix estimator at frequency f k = k−1 F , k = 1, 2, · · · , F , and X k parametrizes the unknown inverse , which are then used to infer the underlying graph. We will now relate (34) with F frequencies, to (18) with n frequencies.
In the analysis of [33], F is kept fixed while the sample size n is allowed to increase. The estimateŜ k is obtained via Blackman-Tukey method which is mathematically equivalent to frequency-domain weighted smoothing of periodgram where the weights w and effective width W depend upon the window function used in the time-domain on estimated correlation function (see [33,Eqn. (8)] and also discussion in [49, Sec. II-D.2] and [24, Sec. 2.5.1]). A standard assumption is that of local smoothness of the true PSD matrix S x (f k ), hence of X k , as in (24). This allows us to rewrite (34) as The total number of DFT frequencies used in (36), hence in (34), is (2W + 1)F whereas that in (18) is n. ForŜ k to be a consistent estimator of S x (f k ), one must have 2W + 1 → ∞ and 2W +1 n → 0 as n → ∞ [41, Secs. 5.6 and 7.4]: 2W +1 n → 0 makesŜ k asymptotically unbiased and 2W + 1 → ∞ makes estimator covariance matrix tend to zero. To have the same number of frequencies in (34) (or (36)) and (18), one must have (2W + 1)F ≈ n, i.e., F ≈ n/(2W + 1) (→ ∞ for consistency), which is not possible for the approach of [33] with fixed F . Note that the analysis of [33] requiresŜ k to be a consistent estimator of S x (f k ), which is possible by letting 2W + 1 → ∞ and 2W +1 n → 0. But with fixed F , the product (2W + 1)F n for large n since F is fixed. Thus, for large n, [33] exploits only a subset of the sufficient statistic set, whereas in our approach the product KM ≈ n 2 (K and M correspond to 2W + 1 and F , respectively, in (36)), thereby using the entire sufficient statistic set. As noted in Sec. V (after (59)), in our case, as n → ∞, we have K → ∞, M → ∞ and K n → 0.

III. PENALIZED LOG-LIKELIHOOD
We wish to estimate inverse PSD matrix Φ k := S −1 x (f k ). In terms of Φ k we rewrite (33) as where the last expression in (37) follows by specifying the pdf ofĎ in terms of joint pdf ofĎ andĎ * as in (10). Then we have the log-likelihood (up to some constants) [37] ln In the high-dimensional case of K < p 2 (number of realvalued unknowns in S −1 x (f k ))), one has to use penalty terms to enforce sparsity and to make the problem well-conditioned. Consider minimization of a convex objective function g(θ) w.r.t. θ ∈ R p . If θ is known to be sparse (only a few nonzero entries), one may choose θ to minimize a lassopenalized cost g(θ) + λ θ 1 where λ > 0 is the lasso penalty or tuning parameter used to control sparsity of the solution [51]. Suppose θ has M groups (subvectors) θ (m) , m = 1, 2, · · · , M , where rather than just sparsity in θ, one would like a solution which uses only a few of the groups θ (m) (group sparsity). To this end, [52] proposed a group lasso penalty where θ is chosen to minimize a group lasso penalized cost g(θ) + λ M m=1 θ (m) 2 and where λ > 0 is the group-lasso penalty parameter. As noted in [53], [54], while the group-lasso gives a sparse set of groups, if it includes a group in the model, then all coefficients in the group will be nonzero. To enforce sparsity of groups and within each group, sparse-group lasso framework has been proposed in [53], [54] where θ is chosen to minimize a sparse-group lasso penalized cost g(θ)+αλ provides a convex combination of lasso and grouplasso penalties. Note that α = 0 gives the group-lasso fit while α = 1 yields the lasso fit, thus sparse-group lasso penalty is more general than either lasso or group-lasso penalties.
Lasso penalty has been used in [22], [55]- [57], group lasso has been used in [58] and sparse-group lasso has been used in [4], [59], all for graphical modeling of real-valued random vectors (i.e., i.i.d. time series) in various contexts. Group lasso has been used in [33] for graphical modeling of dependent time series. Results of [53], [54] for a regression problem and that of [59] for graphical modeling of random vectors in a multi-attribute context (a random vector is associated with each node of a graph instead of just a random variable), both show significant performance improvements over either just lasso or just group lasso penalties. For our problem in this paper we will use sparse-group lasso penalty.
Imposing a sparse-group sparsity constraint, we propose to minimize a penalized version of negative log-likelihood w.r.t. where λ > 0 and α ∈ [0, 1]. In (41), an 1 penalty term is (lasso), and to the off-block-diagonal group of M terms via . To optimize L SGL ({Φ}), using variable splitting, one may reformulate as in [37]: Using the penalty method, [37] considers the relaxed problem (ρ > 0 is "large") where it is solved via an alternating minimization (AM) based method [60]. The final result depends on ρ and strictly speaking, one must have ρ → ∞ which can make the problem numerically ill-conditioned. In the numerical example of [37], a fixed value ρ = 10 was considered which does not necessarily achieve the constraint Φ k − W k = 0. On the other hand, in the ADMM approach which also uses a penalty parameter ρ, the final solution to the optimization problem does not depend upon the chosen ρ > 0, although the convergence speed depends on it [61]. In this paper we consider ADMM. Some authors [62], [63] have suggested that when both AM and ADMM approaches are applicable, the ADMM which requires dual variables, is inferior to the AM method which is a primal-only method, in terms of computational complexity and accuracy, but their claims do not account for the need to solve the AM problem multiple times, each time with increased value of ρ, and moreover, they do not consider graphical modeling problems.

A. ADMM Algorithm
Given the results {Φ (m) }, {W (m) }, {U (m) } of the mth iteration, in the (m + 1)st iteration, an ADMM algorithm executes the following three updates: ) up to some terms not dependent upon Φ k 's, where As in [61, Sec. 6.5] but accounting for complex-valued vectors/matrices in this paper compared to real-valued vectors/matrices in [4], and therefore using the Wirtinger calculus, the solution to arg min ) is as follows. A necessary and sufficient condition for a global optimum is that the gradient of The solution to (48) Let V DV H denote the eigen-decomposition of the Hermitian matrixŜ k − ρA where D is diagonal with real values on the diagonal, and V V H = V H V = I. Then we have whereD := V H Φ k V . Assume thatD is a diagonal matrix and solve (50) for diagonalD. That is,D should satisfy The solutionD satisfies (51) and yieldsD > 0 for any ρ > 0. Therefore, so constructedD 0, and hence, (48) and (49).
Here we use Lemma 1 first stated (but not proved) in [37], and is based on the real-valued results of [53]. Lemma 1 is proved in Appendix A.
where (b) + := max(0, b), soft-thresholding operator S(b, β) := (1 − β/|b|) + b (for complex scalar b = 0), and vector operator [S(a, β)] j = S(a j , β), . Invoking Lemma 1, the solution to minimization of (53) is For the scaled Lagrangian formulation of ADMM [61], for k = 1, 2, · · · , M , update U  (43) is done under the equality constraints Φ k − W k = 0, k = 1, 2, · · · , M . The error in this equality during ADMM iterations is called primal residual as it measures the primal feasibility [61,Sec. 3.3]. The primal residual matrix is given by and the primal residual vector is r p = vec(R p ) ∈ C p 2 M , vectorization of R p . At the (m + 1)st iteration, the primal residual matrix will be denoted by ). As m → ∞, one must have r (m+1) p → 0. Based on some dual feasibility conditions for optimization of the ADMM problem, a dual residual at the (m + 1)st iteration is defined in [61,Sec. 3.3]. For our problem, the dual residual matrix at the (m + 1)st iteration is given , and the dual residual vector is r The convergence criterion is met when the norms of these residuals are below primary and dual tolerances τ pri and τ dual , respectively: for some µ > 1. As stated in [61, Sec. 3.4.1], "The idea behind this penalty parameter update is to try to keep the primal and dual residual norms within a factor of µ of one another as they both converge to zero." For all numerical results presented later, we used ρ (0) = 2 (initial value of ρ), µ = 10, and τ abs = τ rel = 10 −4 . Convergence: The objective function L SGL ({Φ}), given by (40), is strictly convex in {Φ} for Φ k 0, k = 1, 2, · · · , M . It is also closed, proper and lower semi-continuous. Hence, for any fixed ρ > 0, the ADMM algorithm is guaranteed to converge [61,Sec. 3.2], in the sense that we have primal residual convergence to 0, dual residual convergence to 0, and objective function convergence to the optimal value. For varying ρ, the convergence of ADMM has not yet been proven

B. Further on Comparison with Existing Works
Here we briefly summarize comparisons with [33]- [35]. As noted earlier in Sec. II-A, {d x (f m )} n/2 m=0 is a frequencydomian sufficient statistic for this problem. However, in [33] (and [34], [35]), one uses some estimate of PSD S x (f m ), or d x (f m ), for m = 0, 1, · · · , n − 1 appealing to Whittle approximation which, as discussed in Sec. II-D, has lots of redundant frequencies. Also, in [33], frequencies (f m orf k in our notation) are fixed a priori as (f − 1)/F ∈ [0, 1), f = 1, 2, · · · , F , for some even integer F . For instance, in the simulation example of [33], F = 4, leading to four frequencies {0, 0.25, 0.5, 0.75} for any data size n. Note that S x (0.75) = S * x (0.25) (and so are their estimates), so there is no new information in it. Also, S x (0) and S x (0.5), and their estimates, are real-valued, not complex, but as their estimates are based on d x (f k )'s for k in a neighborhood of m = 0 or 0.5, any information in the imaginary part of d x (f k )'s is not exploited. Furthermore, as noted in item (i) in Sec. I-B, [33] considers only group-lasso penalty which is subsumed by our more general sparse-group lasso. In our analysis presented later in Sec. V, the number of frequencies are allowed to increase with sample size, and as discussed in Remark 1, Sec. II-E, increasing the number of frequencies allows one to exploit the entire sufficient statistic set, unlike the analysis in [33] where with fixed F , one uses only a subset of sufficient statistics. where 2KM are total number of real-valued measurements in frequency-domain and 2K are number of real-valued measurements per frequency point, with total M frequencies in (0, π). Each nonzero off-diagonal element ofΦ k consists of two real variables, but sinceΦ k is Hermitian, the number of (nonzero) real unknowns inΦ k equal the number of nonzero elements ofΦ k . Pick α and λ to minimize BIC. We use BIC to first select λ over a grid of values with fixed α, and then with selected λ, we search over α values in [0, 0.3]. This sequential search is computationally less demanding than a two-dimensional search.

C. BIC for Tuning Parameter Selection
We search over λ in the range [λ , λ u ] selected via the following heuristic (as in [59]). The heuristic lies in the fact that we limit the range of α and λ values over which the search is performed. For α = α 0 (=0.1), we first find the smallest λ, labeled λ sm , for which we get a no-edge model (i.e., |Ê| = 0, whereÊ denotes the estimated edge set based on (56)). Then we set λ u = λ sm /2 and λ = λ u /10. The given choice of λ u precludes "extremely" sparse models while that of λ precludes "very" dense models. The choice α ∈ [0, 0.3] reflects the fact that group-lasso penalty across all frequencies is more important ((i, j)th element of inverse PSD at all frequencies must be zero for edge {i, j} ∈ E) than lasso penalty at individual frequencies.

V. CONSISTENCY
In this section we analyze consistency of the proposed approach. In Theorem 1 we provide sufficient conditions for convergence in the Frobenius norm of the inverse PSD estimators to the true value, jointly across all frequencies, where the number of frequencies are allowed to increase with sample size. As discussed in Remark 1, Sec. II-E, increasing the number of frequencies allows one to exploit the entire sufficient statistic set, unlike the analysis in [33]. Theorem 1 also yields a rate of convergence w.r.t. sample size n. We follow proof technique of [22] which deals with i.i.d. time series models and lasso penalty, to establish our main result, Theorem 1.
Define p × (pM ) matrix Ω as With 0 ≤ α ≤ 1, re-express the objective function (40) as where we now allow p, M , K (see (22), (23)), and λ to be functions of sample size n, denoted as p n , M n , K n and λ n , respectively. We take p n to be a non-decreasing function of n, as is typical in high-dimensional settings. Note that K n M n ≈ n/2. Pick K n = a 1 n γ and M n = a 2 n 1−γ for some 0.5 < γ < 1, 0 < a 1 , a 2 < ∞, so that both M n /K n → 0 and K n /n → 0 as n → ∞ (cf. Remark 1). As discussed in Remark 2 later, Theorem 1 clarifies how to choose M n and K n (or put additional restrictions on it) so that for given {p n }, the estimate of Ω converges to its true value in the Frobenius norm. Assume (A2) Define the true edge set of the graph by E 0 , implying (22), and use Ω 0 to denote the true value of Ω). Assume that card(E 0 ) = |(E 0 )| ≤ s n0 . (A3) The minimum and maximum eigenvalues of p n × p n PSD S 0 (f ) 0 satisfy Therefore, by (61), for 1 ≤ i ≤ p n , Finally, by [41, p. 279], LetΩ λ = arg min Ω : Φ k 0 L SGL (Ω). Theorem 1 whose proof is given in Appendix B, establishes consistency ofΩ λ . Theorem 1 (Consistency). For τ > 2, let where Given any real numbers δ 1 ∈ (0, 1), δ 2 > 0 and C 1 > 0, let Suppose the regularization parameter λ n and α ∈ [0, 1] satisfy Then if the sample size is such that n > max{N 2 , N 3 } and assumptions (A1)-(A3) hold true,Ω λ satisfies with probability greater than 1 − 1/p τ −2 n . In terms of rate of convergence (i.e., for large n), A sufficient condition for the lower bound in (69) to be less than the upper bound for every α ∈ [0, 1] is C 1 = 2(1 + α( √ M n − 1)). • Remark 2. Theorem 1 helps determine how to choose M n and K n so that for given {p n }, lim n→∞ Ω λ − Ω 0 F = 0. This behavior is governed by (71), therefore we have to examine ln(M n ) √ M n r n . As noted before, since K n M n ≈ n/2, if one picks K n = a 1 n γ , then M n = a 2 n 1−γ for some 0 < γ < 1, 0 < a 1 , a 2 < ∞. Suppose that the maximum number of nonzero elements in S −1 0 (f ), p n + s n0 , satisfy p n + s n0 = a 3 n θ for some 0 ≤ θ < 1, 0 < a 3 < ∞. Then we have If θ = 0 (fixed graph size and fixed number of connected edges w.r.t. sample size n), then we need 2 3 < γ < 1. By (72), we must have 1 > γ > 2 3 + θ 3 . If θ > 0, γ has to be increased beyond what is needed for θ = 0, implying more smoothing of periodogram d (29)), leading to fewer frequency test points M n . Clearly, we cannot have θ ≥ 1 because p n + s n0 = O(n θ ) will require γ > 1 which is impossible.

VI. NUMERICAL EXAMPLES
We now present numerical results for both synthetic and real data to illustrate the proposed approach. In synthetic data example the ground truth is known and this allows for assessment of the efficacy of various approaches. In real data example where the ground truth is unknown, our goal is visualization and exploration of the dependency structure underlying the data.

A. Synthetic Data
Consider p = 128, 16 clusters (communities) of 8 nodes each, where nodes within a community are not connected to any nodes in other communities. Within any community of 8 nodes, the data are generated using a vector autoregressive (VAR) model of order 3. Consider community q, q = 1, 2, · · · , 16. Then x (q) (t) ∈ R 8 is generated as with w (q) (t) as i.i.d. zero-mean Gaussian with identity covariance matrix. Only 10% of entries of A if not, we re-draw randomly till this condition is fulfilled. The overall data x(t) is given by First 100 samples are discarded to eliminate transients, and generate stationary Gaussian data. This set-up leads to approximately 3.5% connected edges. In each run, we calculated the true PSD S(f ) for f ∈ [0, 0.5] at intervals of 0.01, and then take ii (f ))/p, averaged over 100 runs, turns out to be 75.10 (±1.85). Therefore, the threshold of 10 −6 in f |S −1 ij (f )| > 10 −6 is quite low, resulting in some very "weak" edges in the graph.
Simulation results are shown in Fig. 1 where the performance measure is F 1 -score for efficacy in edge detection. The F 1 -score is defined as F 1 = 2 × precision × recall/(precision + recall) where precision = |Ê ∩ E 0 |/|Ê|, recall = |Ê ∩ E 0 |/|E 0 |, and E 0 andÊ denote the true and estimated edge sets, respectively. For our proposed approach, we consider three different values of M ∈ {2, 4, 6} for five samples sizes n ∈ {128, 256, 512, 1024, 2048}. For M = 2, we used K = 31, 63, 127, 255, 511 for n = 128, 256, 512, 1024, 2048, respectively, for M = 4, we used K = 15, 31, 63, 127, 255 for n = 128, 256, 512, 1024, 2048, respectively, and for M = 6, we used K = 9, 21, 41, 85, 169 for n = 128, 256, 512, 1024, 2048, respectively. These approaches are labeled as "proposed: M=2," "proposed: M=4," and "proposed: M=6," in Fig. 1. The tuning parameters λ and α were selected by searching over a grid of values to maximize the F 1 -score (over 100 runs). The search for α was confined to [0, 0.3]. For a fixed α = 0.1, we first picked the best λ value, and then with fixed best λ value, search over α ∈ [0, 0.3]. Fig. 1 shows the results for thus optimized (λ, α). In practice, one cannot calculate the F 1 -score since ground truth is unknown. For M = 4 we selected (λ, α) in each run via BIC as discussed in Sec. IV-C (where knowledge of the ground truth is not needed). The obtained results based on 50 runs are shown in Fig. 1, labeled as "proposed: M=4,BIC." The conventional i.i.d. modeling approach exploits only the sample covariance 1 n n−1 t=0 x(t)x (t) whereas the proposed approaches exploit the entire correlation function (equivalently PSD), and thus, can deliver better performance. In Fig. 1, the label "IID model" stands for the ADMM lasso approach ( [61, Sec. 6.4]) that models data as i.i.d., and the corresponding results are based on 100 runs with lasso parameter λ selected by exhaustive search over a grid of values to maximize F 1 score. We also show the results of the ADMM approach of [33], labeled "GMS" (graphical model selection), which was applied with F = 4 (four frequency points, corresponds to our M = 4) and all other default settings of [33] to compute the PSDs (see also Sec. IV-B). The lasso parameter λ for [33] was selected by exhaustive search over a grid of values to maximize F 1 score. The F 1 -scores are shown in Fig. 1 and average timings per run are shown in Fig. 2 for sample sizes n = 128, 256, 512, 1024, 2048. All ADMM algorithms were implemented in MATLAB R2020b, and run on a Window Home 10 operating system with processor Intel(R) Core(TM) i5-6400T CPU @2.20 GHz with 12 GB RAM. It is seen from Fig. 1 that with F 1 -score as the performance metric, our proposed approach for all three values of M (number of normalized frequency points in (0, 0.5)) significantly outperforms the conventional IID model approach. The BIC-based approach for M = 4 yields performance that is close to that based on optimized parameter selection for n ≥ 512. It is also seen that [33] performs better than IID modeling but much worse than our proposed sparse group lasso approach, while also taking more time to convergence. The stopping rule, variable penalty and thresholds selected for all approaches were the same since all approaches are ADMM-based; these values have been specified as ρ (0) = 2, µ = 10, and τ abs = τ rel = 10 −4 in Sec. IV-A5.
The conventional i.i.d. modeling approach estimates the (sparse) precision matrix Ω = E{x(t)x (t)} −1 = R −1 xx (0): there is an edge {i, j} in CIG iff Ω ij = 0. Since this approach ignores R xx (m) for m = 0 for dependent data, its performance is the worst for all sample sizes, although the performance does improve with n since {i, j} ∈ E ⇒ Ω ij = 0 and accuracy of the estimate of Ω ij improves with increasing n. The method of [33] performs better than IID modeling since it does use R xx (m) for m = 0 in estimating S −1 xx (f ). Also, the performance of [33] improves with n as accuracy of the estimates of the PSD improves with n. However, as summarized in Sec. IV-B, [33] makes some peculiar choices which are likely reasons why its performance is inferior to our proposed approach.
Note that for our example, there is no explicit mathematical expression for calculating the true edge set E 0 . In each run, we calculated the true PSD S(f ) for f ∈ [0, 0.5] at intervals of 0.01, and then took {i, j} ∈ E if f |S −1 ij (f )| > 10 −6 , where the average value of diagonal elements ( ii (f ))/p, averaged over 100 runs, turns out to be 75.10 (±1.85). That is, we have some very "weak" edges in the graph which are not easy to detect with relatively "short" sample sizes, resulting in relatively low F 1 scores. Nevertheless, when comparing different approaches, our proposed approach performs much better.

VII. CONCLUSIONS
Graphical modeling of dependent Gaussian time series was considered. A sparse-group lasso-based frequency-domain formulation of the problem was proposed and analyzed where the objective was to estimate the inverse PSD of the data via optimization of a sparse-group lasso penalized log-likelihood cost function. The graphical model is then inferred from the estimated inverse PSD. We investigated an ADMM approach for optimization. We established sufficient conditions for convergence in the Frobenius norm of the inverse PSD estimators to the true value, jointly across all frequencies, where the number of frequencies were allowed to increase with sample  size. We also empirically investigated selection of the tuning parameters based on the Bayesian information criterion, and illustrated our approach using numerical examples utilizing both synthetic and real data. The synthetic data results show that for graph edge detection, the proposed approach significantly outperformed the widely used i.i.d. modeling approach where the underlying time series is either assumed to be i.i.d., or one uses only the covariance of the data. The proposed approach also outperformed the approach of [33] for synthetic data.
Lemma 2 below is specialization of [57, Lemma 1] to Gaussian random vectors. It follows from [57, Lemma 1] after setting the sub-Gaussian parameter σ in [57, Lemma 1] to 1. Lemma 2. Consider a zero-mean Gaussian random vector z ∈ R p with covariance R 0. Given n i.i.d. samples z(t), t = 0, 1, · · · , n − 1, of z, letR = (1/n) n−1 t=0 zz denote the sample covariance matrix. ThenR satisfies the tail bound for all δ ∈ (0, 40 max i (R ii )) • Exploiting Lemma 2, we have Lemma 3 regardingŜ k . We denote S 0 (f k ) as S 0k in this section. A proof of Lemma 3 is in [38, Lemma 2]. The statement of Lemma 3 below is a corrected version of [38,Lemma 2] with no changes to its proof given therein. Lemma 3. Under Assumption (A2),Ŝ k satisfies the tail bound for τ > 2, if the sample size n and choice of K n is such that K n > N 1 = 2 ln(16p τ n M n ), where C 0 is defined in (63). • Lemma 4 deals with a Taylor series expansion using Wirtinger calculus.
0 and Γ k = Γ H k . Then using Wirtinger calculus, the Taylor series expansion of c(Φ k , Φ * k ) is given by where h.o.t. stands for higher-order terms in Γ k and Γ * k . • Proof: Only for the proof of this lemma, we will drop the subscript k, and donate Φ k , Φ 0k and Γ k as Φ, Φ 0 and Γ, respectively. Treating Φ and Φ * as independent variables, the Taylor series expansion of c(Φ, Φ * ) is where Consider the following facts [65], [66] ∂ ln |Φ| Using the above partial derivatives in the first derivative terms in the Taylor series, we have The quadratic terms in the Taylor series yield q,r s,t Given matrices A and B for which product AB is defined, and additionally given matrix Y such that product AY B is defined, we have vec(AY B) = (B ⊗ A)vec(Y ) and tr(AB) = (vec(A)) vec(B). Using these results we have Using (104), we rewrite terms in (103) as In (105) Lemma 4 regarding Taylor series expansion immediately leads to Lemma 5 regarding Taylor series with integral remainder, needed to follow the proof of [22] pertaining to the real-valued case.
Therefore, if we can show that the minimizer∆ must be inside the sphere defined by Θ n (R), and hence ∆ F ≤ Rr n .
We first bound A 1k 's and A 1 . Note that H(Φ 0k , Γ k , v) is a Hermitian matrix and its eigenvalues consist of the eigenvalues of H 11k and H 22k . Since the eigenvalues of A ⊗ B are the product of the eigenvalues of A and eigenvalues of B for Hermitian A and B, the eigenvalues of H 11k are the product of the eigenvalues of (Φ 0k + vΓ k ) − * and that of (Φ 0k + vΓ k ) −1 . But these two matrices have the same set of eigenvalues since one matrix is the complex conjugate of the other, and both have real eigenvalues since they are Hermitian. Since H 11k = H * 22k , it follows that the eigenvalues of H 11k are the same as the eigenvalues of H 22k . Thus where we have used the facts that we have Thus, where we have used the fact that Φ 0k = S −1 0k = φ max (S −1 0k ) = (φ min (S 0k )) −1 ≤ β −1 min and Γ k ≤ Γ k F ≤ ∆ F = Rr n = O(r n ). Therefore, where (134) To bound L 21k , using Lemma 3, with probability > 1 − 1/p τ −2 n , Using Cauchy-Schwartz inequality and Lemma 3, with probability > 1 − 1/p τ −2 n , where s n0 is the cardinality of the true edge set E 0 (see Assumption (A1)). Thus, with probability > 1 − 1/p τ −2 n , Hence with A 2 = Mn k=1 (1/2)A 2k , We now derive an alternative bound on A 2 . We have w.h.p.
where∆ ∈ R pn×pn has its (i, j)th element∆ ij = ∆ (ij) F . We now bound A 3k . Let E c 0 denote the complement of E 0 , given by E c 0 = {{i, j} : [S −1 0 (f )] ij ≡ 0, i = j, 0 ≤ f ≤ 0.5}. For an index set B and a matrix C ∈ C p×p , we write C B to denote a matrix in C p×p such that [ leading to (A 3 = Mn k=1 A 3k ) Similarly, By Cauchy-Schwartz inequality, Set Γ − k 1 = Γ − kE0 1 + Γ − kE c 0 1 in A 2 of (138) to deduce that w.h.p.