Low Rank and Structured Modeling of High-dimensional Vector Autoregressions

TL;DR

This paper introduces a novel regularized approach for estimating high-dimensional vector autoregressive models with low-rank and structured sparse components, improving network inference in complex time series data.

Contribution

It develops a regularized framework combining nuclear norm and lasso penalties, along with a fast algorithm, to better capture low-rank and structured sparsity in VAR models.

Findings

Improved estimation accuracy over standard sparse VAR methods.

Theoretical bounds on estimation error rates.

Effective performance demonstrated on synthetic and real datasets.

Abstract

Network modeling of high-dimensional time series data is a key learning task due to its widespread use in a number of application areas, including macroeconomics, finance and neuroscience. While the problem of sparse modeling based on vector autoregressive models (VAR) has been investigated in depth in the literature, more complex network structures that involve low rank and group sparse components have received considerably less attention, despite their presence in data. Failure to account for low-rank structures results in spurious connectivity among the observed time series, which may lead practitioners to draw incorrect conclusions about pertinent scientific or policy questions. In order to accurately estimate a network of Granger causal interactions after accounting for latent effects, we introduce a novel approach for estimating low-rank and structured sparse high-dimensional VAR models. We introduce a regularized framework involving a combination of nuclear norm and lasso (or group lasso) penalty. Further, and subsequently establish non-asymptotic upper bounds on the estimation error rates of the low-rank and the structured sparse components. We also introduce a fast estimation algorithm and finally demonstrate the performance of the proposed modeling framework over standard sparse VAR estimates through numerical experiments on synthetic and real data.

Tables

Table 1. Table 1 : Parameter settings in the proposed algorithms for all the experiments.

Parameters	$\sigma$	${\eta }_{\min }$	$C$	${\beta }_k$
Values	2	${\Vert {𝒳}^{\prime }𝒳\Vert }_2/10$	100	1/k

Table 2. Table 2 : Performance comparison of FNSL with a variant of FISTA on large-scale sparse network structure learning problem.

p	N	method	(TPR, FAR)(%)	EE	T
800	1000	FISTA	(88.5, 14.4)	0.22	11.9
		FNSL	(88.6, 13.9)	0.22	$6.1$
	1500	FISTA	(90.8, 12.3)	0.17	14.7
		FNSL	(90.7, 12.0)	0.17	$8.2$
	2000	FISTA	(91.3, 11.0)	0.14	17.6
		FNSL	(91.3, 11.3)	0.14	$10.8$
900	1000	FISTA	(87.8, 15.0)	0.24	13.1
		FNSL	(87.8, 14.7)	0.24	$7.2$
	1500	FISTA	(89.1, 10.7)	0.19	16.3
		FNSL	(89.1, 10.5)	0.19	$8.7$
	2000	FISTA	(90.9, 12.0)	0.16	20.7
		FNSL	(90.9, 11.8)	0.16	$13.5$
1000	1000	FISTA	(88.5, 15.2)	0.25	18.6
		FNSL	(88.4, 14.8)	0.25	$11.2$
	1500	FISTA	(90.3, 14.1)	0.20	21.5
		FNSL	(90.2, 13.8)	0.20	$12.8$
	2000	FISTA	(91.2, 12.1)	0.16	24.7
		FNSL	(91.2, 12.4)	0.16	$15.1$

Table 3. Table 3 : Comparison of objective function value, CPU time in seconds, and the number of matrix products (AX) for a variant of FISTA and FNSL on sparse network structure with different sample size.

Algorithms	Objective value	CPU	AX
	$p=1000,N=1000$
FISTA	4.140e+5	18.6	66
FNSL	4.089e+5	11.2	44
	$p=1000,N=1500$
FISTA	6.137e+5	21.5	62
FNSL	6.071e+5	12.8	38
	$p=1000,N=2000$
FISTA	8.239e+5	24.7	66
FNSL	8.169e+5	15.1	41

Table 4. Table 4 : Performance comparison of a variant of FISTA and FNSL on estimation of low-rank transition matrices problems.

p	N	method	$\hat{r}$	EE	T
200	400	FISTA	8	0.80	4.9
		FNSL	8	0.80	$3.4$
	1200	FISTA	8	0.63	9.7
		FNSL	8	0.63	$6.9$
	2000	FISTA	8	0.57	14.2
		FNSL	8	0.57	$10.5$
300	400	FISTA	12	0.84	7.8
		FNSL	12	0.84	$5.7$
	1200	FISTA	12	0.72	16.1
		FNSL	12	0.72	$11.7$
	2000	FISTA	12	0.68	18.2
		FNSL	12	0.68	$13.8$
400	400	FISTA	16	0.87	8.5
		FNSL	16	0.87	$5.9$
	1200	FISTA	16	0.82	20.5
		FNSL	16	0.82	$17.4$
	2000	FISTA	16	0.75	31.4
		FNSL	16	0.75	$25.1$

Table 5. Table 5 : True positive rate and false alarm rate of the L+S model on identifying the sparse component S with different α 𝛼 \alpha .

$\alpha$	(TPR, FAR)(%)
p/8	(84.5, 17.4)
p/4	(82.4, 17.2)
p/2	(80.4, 17.5)
p	(73.2, 17.1)
2p	(54.7, 17.0)
4p	(40.2, 17.8)
8p	(22.7, 17.4)

Table 6. Table 6 : Performance comparison of L+S with OLS and Lasso.

p	N	model	(TPR, FAR)(%)	EE	PE
50	100	OLS	(-, -)	0.84	0.72
		Lasso	(73.2, 30.0)	0.69	0.53
		L+S	($76.3$, $18.9$)	$0.48$	$0.47$
	200	OLS	(-, -)	0.52	0.41
		Lasso	(77.3, 35.0)	0.57	0.45
		L+S	($80.4$, $17.5$)	$0.31$	$0.36$
75	100	OLS	(-, -)	0.75	0.37
		Lasso	(71.0, 24.7)	0.75	0.37
		L+S	($79.0$, $18.0$)	$0.51$	$0.29$
	200	OLS	(-, -)	0.53	0.18
		Lasso	(77.0, 28.6)	0.67	0.22
		L+S	($83.8$, $18.3$)	$0.36$	$0.16$
100	100	OLS	(-, -)	3.7	4.0
		Lasso	(57.3, 29.0)	1.06	1.05
		L+S	($52.3$, $20.1$)	$0.92$	$1.0$
	200	OLS	(-, -)	2.07	1.73
		Lasso	(59.4, 25.5)	0.86	0.95
		L+S	($60.4$, $20.5$)	$0.72$	$0.90$

Table 7. Table 7 : Performance comparison of L+S+G with S+G and L+S.

p	N	method	(TPR, FAR)(%)	EE	PE
50	200	S+G	(85.5, 34.1)	0.46	0.53
		L+S	(82.6, 26.4)	0.41	0.51
		L+S+G	($83.3$, $26.9$)	$0.41$	$0.51$
	300	S+G	(91.7, 47.0)	0.37	0.58
		L+S	(88.6, 24.4)	0.31	0.56
		L+S+G	($90.6$, $24.9$)	$0.30$	$0.56$
100	200	S+G	(92.3, 49.9)	0.48	0.73
		L+S	(84.3, 28.4)	0.44	0.72
		L+S+G	($85.3$, $27.3$)	$0.44$	$0.72$
	300	S+G	(94.8, 49.0)	0.44	0.72
		L+S	(89.6, 25.4)	0.37	0.70
		L+S+G	($90.0$, $25.1$)	$0.36$	$0.70$
150	200	S+G	(92.0, 50.5)	0.64	0.73
		L+S	(83.3, 28.8)	0.57	0.71
		L+S+G	($84.0$, $28.1$)	$0.55$	$0.70$
	300	S+G	(93.6, 50.2)	0.55	0.71
		L+S	(85.6, 27.4)	0.46	0.68
		L+S+G	($86.4$, $27.6$)	$0.46$	$0.68$

Table 8. Table 8 : Performance comparison of S+G with Lasso and SGL on sparse plus group-sparse network identification problem.

p	N	model	(TPR, FAR)(%)	EE
50	100	Lasso	(67.9, 16.0)	0.40
		SGL	(71.1, 18.4)	0.41
		S+G	($86.7$, $16.5$)	$0.39$
	200	Lasso	(69.3, 10.4)	0.32
		SGL	(71.1, 11.6)	0.33
		S+G	($90.4$, $10.4$)	$0.28$
100	100	Lasso	(70.2, 21.2)	0.44
		SGL	(71.6, 23.5)	0.46
		S+G	($84.7$, $21.2$)	$0.43$
	200	Lasso	(87.8, 17.4)	0.34
		SGL	(77.7, 19.7)	0.35
		S+G	($87.8$, $17.4$)	$0.32$
200	100	Lasso	(75.2, 43.4)	0.76
		SGL	(76.7, 43.3)	0.76
		S+G	($84.3$, $43.7$)	$0.74$
	200	Lasso	(70.1, 23.1)	0.55
		SGL	(71.4, 23.5)	0.55
		S+G	($78.5$, $23.5$)	$0.54$

Equations

$$X^t=B^{\prime }{X}^{t-1}+{ϵ}^t,t=1,\cdots ,T,$$

$$X^t=B^{{}^{\prime }}{X}^{t-1}+{ϵ}^t,{ϵ}^t\stackrel{i.i.d.}{\sim }N(0,{\Sigma }_{ϵ}),$$

$$B=L^{\ast }+R^{\ast },rank(L^{\ast })=r,$$

$$f_X(\theta )=\frac{1}{2\pi }\left({\mathcal{B}}^{-1}\left(e^{i\theta }\right)\right){\Sigma }_{ϵ}{\left({\mathcal{B}}^{-1}\left(e^{i\theta }\right)\right)}^{†},\theta \in [-\pi ,\pi ],$$

$$\mathcal{M}(f_X)=\underset{\theta \in [-\pi ,\pi ]}{\sup }{\Lambda }_{\max }(f_X(\theta )),$$

$$\text{EuFrak}m(f_X)=\underset{\theta \in [-\pi ,\pi ]}{\sup }{\Lambda }_{\min }(f_X(\theta )),$$

$${\mu }_{\max }(\mathcal{B})=\underset{|z|=1}{\max }{\Lambda }_{\max }\left({\mathcal{B}}^{†}(z)\mathcal{B}(z)\right),$$

$${\mu }_{\min }(\mathcal{B})=\underset{|z|=1}{\min }{\Lambda }_{\min }\left({\mathcal{B}}^{†}(z)\mathcal{B}(z)\right).$$

$$\text{EuFrak}m(f_X)\ge \frac{1}{2\pi }\frac{{\Lambda }_{\min }({\Sigma }_{ϵ})}{{\mu }_{\max }(\mathcal{B})},\mathcal{M}(f_X)\le \frac{1}{2\pi }\frac{{\Lambda }_{\max }({\Sigma }_{ϵ})}{{\mu }_{\min }(\mathcal{B})}.$$

$${\mu }_{\max }(\mathcal{B})\le {\left[1+l+(v_{in}+v_{out})/2\right]}^2$$

$$\displaystyle\underbrace{\left[\begin{array}[]{c}(X^{T})^{\prime}\\ \vdots\\ (X^{1})^{\prime}\end{array}\right]}_{\mathcal{Y}}$$

$$(\hat{L},\hat{R})=\underset{\begin{array}{c}L,R\in {\mathbb{R}}^{p\times p}\\ L\in \Omega \end{array}}{\mathrm{argmin}}l(L,R),$$

$$l(L,R):=\frac{1}{2}{\Vert \mathcal{Y}-\mathcal{X}(L+R)\Vert }_F^2+{\lambda }_N\Vert L{\Vert }_{\ast }+{\mu }_N\Vert R{\Vert }_{\diamond }.$$

$$\frac{1}{2}\Vert \mathcal{X}\Delta {\Vert }_{\text{F}}^2\ge \frac{\zeta }{2}\Vert \Delta {\Vert }_{\text{F}}^2-{\tau }_N{\Phi }^2(\Delta ),\text{for all}\Delta \in {\Re }^{p\times p}$$

$$\Vert {\mathcal{X}}^{\prime }E/N{\Vert }_2\le \varphi (B,{\Sigma }_{ϵ})\sqrt{p/N},\text{and}$$

$$\Vert {\mathcal{X}}^{\prime }E/N{\Vert }_{\max }\le \varphi (B,{\Sigma }_{ϵ})\sqrt{\frac{2\log p}{N}},\text{and}$$

$$\Vert {\mathcal{X}}^{\prime }E/N{\Vert }_{2,\max }\le \varphi (B,{\Sigma }_{ϵ})\frac{\sqrt{m\log K}}{\sqrt{N}},$$

$$\Vert \hat{L}-L^{\ast }{\Vert }_F^2+\Vert \hat{S}-S^{\ast }{\Vert }_F^2\le \frac{4}{{\zeta }^2}(\frac{9}{2}{\lambda }_N^{2}r+4{\mu }_N^{2}s).$$

$$\Vert \hat{L}-L^{\ast }{\Vert }_F^2+\Vert \hat{G}-G^{\ast }{\Vert }_F^2\le \frac{4}{{\zeta }^2}(\frac{9}{2}{\lambda }_N^{2}r+4{\mu }_N^{2}g)$$

$$\Vert \hat{S}+\hat{G}-S^{\ast }-G^{\ast }{\Vert }_F^2\le \frac{4}{{\zeta }^2}(8{\mu }_N^{2}s+9{\nu }_N^{2}g).$$

$$\Vert \hat{L}-L^{\ast }{\Vert }_F^2+\Vert \hat{S}+\hat{G}-S^{\ast }-G^{\ast }{\Vert }_F^2\le \frac{4}{{\zeta }^2}(9{\lambda }_N^{2}r+\frac{25}{2}{\mu }_N^{2}s+8{\nu }_N^{2}g).$$

$$\varphi (B,{\Sigma }_{ϵ})={\Lambda }_{\max }({\Sigma }_{ϵ})\left[1+\frac{1+{\mu }_{\max }(\mathcal{B})}{{\mu }_{\min }(\mathcal{B})}\right]$$

$$\mathbb{P}\left[\Vert \frac{{\mathcal{X}}^{\prime }E}{N}{\Vert }_2>c_0\varphi (B,{\Sigma }_{ϵ})\sqrt{\frac{p}{N}}\right]\le c_1\exp \left[-c_2\log p\right]$$

$$\mathbb{P}\left[\Vert \frac{{\mathcal{X}}^{\prime }E}{N}{\Vert }_{\max }>c_0\varphi (B,{\Sigma }_{ϵ})\sqrt{\frac{\log p}{N}}\right]\le c_1\exp \left[-c_2\log p\right]$$

$$\mathbb{P}\left[\Vert \frac{{\mathcal{X}}^{\prime }E}{N}{\Vert }_{2,\max }>c_0\varphi (B,{\Sigma }_{ϵ})\frac{\sqrt{m\log p}}{\sqrt{N}}\right]\le c_1\exp \left[-c_2\log p\right]$$

$$\mathbb{P}\left[{\Lambda }_{\min }(\frac{{\mathcal{X}}^{\prime }\mathcal{X}}{N})>\frac{{\Lambda }_{\min }({\Sigma }_{ϵ})}{2{\mu }_{\max }(\mathcal{B})}\right]\le c_1\exp \left[-c_2\log p\right]$$

$$\Vert \hat{S}-S^{\ast }{\Vert }_F^2+\Vert \hat{L}-L^{\ast }{\Vert }_F^2\le \frac{c_0{\varphi }^2(B,{\Sigma }_{ϵ}){\mu }_{\max }^2(\mathcal{B})}{{\Lambda }_{\min }^2({\Sigma }_{ϵ})}\frac{(rp+s\log p)}{N}+\frac{32{\Lambda }_{\min }^2({\Sigma }_{ϵ})}{{\mu }_{\max }^2(\mathcal{B})}\frac{s{\alpha }^2}{p^2}.$$

$$\Vert \hat{G}-G^{\ast }{\Vert }_F^2+\Vert \hat{L}-L^{\ast }{\Vert }_F^2\le \frac{c_0{\varphi }^2(B,{\Sigma }_{ϵ}){\mu }_{\max }^2(\mathcal{B})}{{\Lambda }_{\min }^2({\Sigma }_{ϵ})}\frac{(rp+g(m\log p))}{N}+\frac{32{\Lambda }_{\min }^2({\Sigma }_{ϵ})}{{\mu }_{\max }^2(\mathcal{B})}\frac{g{\beta }^2}{K}.$$

$${\eta }_{0,i}=\max \left\{{\eta }_{\min },\frac{\Vert \mathcal{X}(B_i-B_{i-1}){\Vert }_F^2}{\Vert B_i-B_{i-1}{\Vert }_F^2}\right\}\text{for}i>1.$$

$$B_i^{md}=$$

$$B_{i+1}=$$

$${\Gamma }_i=$$

$$Q_{i+1}=$$

$$B_{i+1}^{ag}=(1-{\alpha }_i)B_i^{ag}+{\alpha }_i{B}_{i+1}.$$

Full text

Keywords: lasso, group lasso, nuclear norm, low rank, vector autoregression, probabilistic bounds, identifiability, fast algorithm.

1 Introduction

The problem of learning the network structure among a large set of time series arises in many signal processing, economic, finance and biomedical applications. Examples include processing signals obtained from radars [1, 2], macroeconomic policy making and forecasting [3], assessing connectivity among financial firms [4], reconstructing gene regulatory interactions from time-course genomic data [5] and understanding connectivity between brain regions from fMRI measurements [6]. Vector Autoregressive (VAR) models provide a principled framework for these tasks.

Formally, a VAR model for $p$-dimensional time series $X_{t}$ is defined in its simplest possible form involving a single time-lag as

[TABLE]

where $B$ is a $p\times p$ transition matrix specifying the lead-lag cross dependencies among the $p$ time series and $\{\epsilon^{t}\}$ is a zero mean error process. VAR models for small number of time series (low-dimensional) have been thoroughly studied in the literature [7]. However, the above mentioned applications, where dozens to hundreds of time series are involved, created the need for the study of VAR models under high dimensional scaling and the assumption that their interactions are * sparse* to compensate for the possible lack of adequate number of time points (samples; see [8] and references therein). There has been a growing body of literature on sparse estimation of large scale VAR models [39], including alternative penalties beyond the popular $\ell_{1}$ penalty (lasso), such as the Berhu regularization introduced in [9], group lasso type penalties employed in [22, 21], as well as non-convex penalties akin to a square-lasso one investigated in [11]. Further, [10] examine estimation of the transition matrix and the inverse covariance matrix of the error process through a joint sparse penalty. Note that the problem of sparse estimation of these two model parameters separately from a least squares and maximum likelihood viewpoints is addressed in [8, 3], respectively, where in addition probabilistic finite sample error bounds for the obtained estimates are obtained.

Nevertheless, there are occasions where the sparsity assumption may not be sufficient. For example, during financial crisis periods, returns on assets move together in a more concerted manner [4, 12], while transcription factors regulate a large number of genes that may lead to hub-node network structures [13]. Similarly, in brain connectivity networks, particular tasks activate a number of regions that cross talk in a collaborative manner [14]. Hence, it is of interest to study VAR models under high dimensional scaling where the transition matrix governing the temporal dynamics exhibits a more complex structure; e.g. it is * low rank* and/or * (group) sparse*.

In a low-dimensional regime, where the number of time points scales to infinity, but the number of time series under study remains fixed, [15] examined asymptotic properties of VAR models, where the parameters exhibit reduced rank structure and also discussed connections with canonical correlation analysis of such models presented in [16]. Specifically, the transition matrix $B$ in (2) can be written as the product of two rank-$k$ matrices $\Phi,\Psi$, i.e. $B=\Phi\Psi^{\prime}$, so that in the resulting model specification of the original $p$ time series is expressed as linear combinations $Z^{t}=\Psi X^{t}$ of the original ones, and $\Phi$ specifies the dependence between $X^{t}$ and $Z^{t}$; namely $X^{t}=\Phi^{\prime}Z^{t-1}+\epsilon^{t}$. Hence, to obtain $\Phi$ and $\Psi$ [15] suggest to estimate the parameters of the original model in (1) under the constraint that $B=\Phi\Psi$ and that rank$(B)=k$. Other works include low rank approximations of Hankel matrices that represent the input-output structure of a linear time invariant systems and were studied in [17, 18]. Finally, a brief mention to the possibility that the VAR transition matrix may exhibit such a structure appeared as a motivating example in [19].

On the other hand, there is a mature literature on imposing low rank plus sparse, or pure group sparse structure for many learning tasks for independent and identically distributed (i.i.d.) data. Examples include group sparsity in regression and graphical modeling [20], low rank and sparse matrix approximations for dimension reduction [18], etc. However, as shown in [8], the presence of temporal dependence across observations induces intricate dependencies between both rows and columns of the design matrix of the corresponding least squares estimation problem, as well as between the design matrix and the error term, that require careful handling to establish consistency properties for the model parameters under sparsity and high dimensional scaling. These issues are further compounded when more complex regularizing norms are involved, as discussed in [21]. In this paper, the authors model grouping structures within each column of $B$, but do not consider a low-rank component. In contrast, we focus on groups potentially spanning across different columns and allow a low-rank component in $B$.

More recently, [37] and [38] extended the framework of [18] beyond decomposition of an observable matrix to Gaussian process identification by assuming a low-rank plus sparse structure on the inverse spectral density and the transfer function of a general VAR($d$) system, respectively. Our work is complementary to this recent developments. We directly model the transition matrix of a VAR(1) process that enables us to identify a directed network of (group) sparse Granger causal relationships that are of interest in a number of applications; e.g. in financial economics where firms with higher out-degree are of particular interest in measuring systemic risk [12, 4]. Further, our approach explicitly address identifiability issues for extracting the respective low-rank and sparse components, which in turn are leveraged to obtain probabilistic error bounds that characterize the quality of their estimates. The latter provide insights to the practitioner on sample size requirements and tuning parameter selection for real data applications. Finally, note that our approach to the issue of identifiability builds on [19], wherein we characterize the degree of unidentifiability which guides in an explicit manner the selection of the tuning parameters used in the proposed optimization algorithm.

Further, to estimate the posited model in (1) with $B$ being low-rank and structured sparse (henceforth indicating that it could be either * pure sparse or group sparse or both*), we also introduce a fast accelerated proximal gradient algorithm, inspired by [30, 31], for the corresponding optimization problems. The key idea is that instead of searching for the local Lipschitz constant of the gradient of the smooth component of the objective function, the proposed algorithm utilizes a safeguarded Barzilai-Borwein (BB) initial stepsize [25] and employs relaxed line search conditions to achieve better performance in practice. The latter enables the selection of more “aggressive” stepsizes, while preserving the accelerated convergence rate of $\mathcal{O}(\frac{1}{k^{2}})$, where $k$ denotes the number of iterations required until convergence. Finally, the performance of the model parameters under different structures together with the associated estimation procedure based on the accelerated proximal gradient algorithm are calibrated on synthetic data, and illustrated on three data sets examining realized volatilities of stock prices of $75$ large financial firms before, during and after the $2007-09$ US financial crisis.

Notation: Throughout the paper, we employ the following notation: $\|.\|$, $\|.\|_{2}$ and $\|.\|_{F}$ denote the $\ell_{2}$-norm of a vector, the spectral norm and the Frobenius norm of a matrix, respectively. For a $p\times p$ matrix $B$, the symbol $\|B\|_{*}$ is used to denote the nuclear norm, i.e. $\sum_{j=1}^{p}\sigma_{j}(B)$, the sum of the singular values of a matrix, while $B^{\dagger}$ denotes the conjugate transpose of a matrix $B$. For any matrix $B$, we use $\|B\|_{0}$ to denote $card(vec(B))$, $\|B\|_{1}$ for $\|vec(B)\|_{1}$ and $\|B\|_{\max}$ to denote $\|vec(B)\|_{\infty}$. Further, if $\{G_{1},G_{2},\ldots,G_{K}\}$ denote a partition of $\{1,2,\cdots,p^{2}\}$ into $K$ non-overlapping groups, then we use $\|B\|_{2,1}$ to denote $\sum_{k=1}^{K}\|(B)_{G_{k}}\|_{F}$, $\|B\|_{2,\max}$ for $\max_{k=1,2,...K}\|(B)_{G_{k}}\|_{F}$, while $\|B\|_{2,0}$ denotes the number of nonzero groups in $B$. Here, with a little abuse of notation, we use $B_{G_{k}}$ to denote $vec(B)_{G_{k}}$. In addition, $\Lambda_{\max}(.)$, $\Lambda_{\min}(.)$ denote the maximum and minimum eigenvalues of a symmetric or Hermitian matrix. For any integer $p\geq 1$, we use $\mathbb{S}^{p-1}$ to denote the unit ball $\{v\in\mathbb{R}^{p}:\|v\|=1\}$. We also use $\{e_{1},e_{2},\ldots\}$ generically to denote unit vectors in $\mathbb{R}^{p}$, when $p$ is clear from the context. Finally, for positive real numbers $A,B$, we write $B\succsim A$ if there exists an absolute positive constant $c$, independent of the model parameters, such that $B\geq cA$.

2 Model Formulation and Estimation Procedure

Consider a VAR(1) model where the transition matrix $B$ is low-rank plus structured sparse given by

[TABLE]

where $L^{*}$ corresponds to the low rank component and $R^{*}$ represents either a sparse $S^{*}$, or group-sparse component $G^{*}$. It is further assumed that the number of non-zero elements in the sparse case is $\|S^{*}\|_{0}=s$, while in the group sparse case the number of non-zero groups is $\|G^{*}\|_{2,0}=g$, with $r\ll p,s\ll p^{2}$ and $g\ll p^{2}$. The matrix $L^{*}$ captures persistence structure across * all* $p$ time series and enables the model to be applicable in settings where there are strong cross-autocorrelations, a feature that the standard sparse VAR model is not equipped to handle. The sparse or group sparse component captures additional cross-sectional autocorrelation structure among the time series. Finally, it is assumed that the error terms are serially uncorrelated. Our objective is to estimate $L^{*}$ and $R^{*}$ accurately based on a relatively small number of samples $N\ll p^{2}$.

Stability. In order to ensure consistent estimation, we assume that the posited VAR model in (2) is stable; i.e. its characteristic polynomial $\mathcal{B}(z):=I_{p}-B^{{}^{\prime}}z$ satisfies $\det(\mathcal{B}(z))\neq 0$ on the unit circle of the complex plane $\{z\in\mathcal{C}:|z|=1\}$. This is a common assumption in the literature of multivariate time series [7], required for consistency and asymptotic normality of low-dimensional VAR models. This assumption also ensures that the spectral density of the VAR model

[TABLE]

is bounded above in spectral norm.

It was shown in [8] that this condition is sufficient to establish consistency of some regularized VAR estimates of a sparse transition matrix. Further, the following quantities play a central role in the error bounds of the regularized estimates:

[TABLE]

As shown in [8], $\mathcal{M}(f_{X})$ and $\EuFrak{m}(f_{X})$ together capture the narrowness of the spectral density of a time series. Processes with stronger temporal and cross-sectional dependence have narrower spectra that in turn lead to slower convergence rates for the regularized estimates. For VAR models, $\mathcal{M}(f_{X})$ and $\EuFrak{m}(f_{X})$ are related to $\mu_{\max}(\mathcal{B})$ and $\mu_{\min}(\mathcal{B})$ as follows:

[TABLE]

Proposition 2.2 in [8] provides a lower bound on $\mu_{\min}(\mathcal{B})$. For the special structure of the models considered here, we can get an improved upper bound on $\mu_{\max}(\mathcal{B})$, as shown in the following lemma:

Lemma 2.1**.**

For a stable VAR(1) model of the class (2), we have

[TABLE]

where $l$ is the largest singular value of $L^{*}$, $v_{in}=\max_{1\leq j\leq p}\sum_{i=1}^{p}|R^{*}_{ij}|$ and $v_{out}=\max_{1\leq i\leq p}\sum_{j=1}^{p}|R^{*}_{ij}|$.

Proof.

$\|\mathcal{B}(z)\|=\|I-(L^{*}+R^{*})z\|\leq\|I\|+\|L^{*}\|+\|R^{*}\|$ for any $z\in\mathbb{C}$ with $|z|=1$. The result follows from the fact that $\mu_{\max}(\mathcal{B})=\max_{|z|=1}\|\mathcal{B}(z)\|^{2}$. ∎

2.1 Estimation Procedure

The estimation of VAR model parameters is based on the following regression formulation (see [7]). Given $T+1$ consecutive observations $\{X^{0},\cdots,X^{T}\}$ from the VAR model, we work with the autoregressive design as follows:

[TABLE]

This is a standard regression problem with $N\equiv T$ samples and $q=p^{2}$ variables. Our goal is to estimate $L^{*}$, $R^{*}$ with high accuracy when $N\ll p^{2}$.

There is an inherent identifiability issue in the estimation of the components $L^{*}$ and $R^{*}$. Suppose the low-rank component $L^{*}$ itself is $s$-sparse or $g$-group sparse and the sparse or group-sparse component $R^{*}$ is of rank $r$. In that scenario, we can not hope for any method to estimate $L^{*}$ and $R^{*}$ separately without imposing any further constraints. So, a minimal condition for low-rank and sparse or group-sparse recovery is that the low rank part should not be too sparse and the sparse or group-sparse part should not be low-rank.

This issue has been rigorously addressed in the literature (e.g. [18]) for independent and identically distributed data and resolved by imposing an * incoherence* condition. Such a condition is sufficient for * exact* recovery of the low rank and the sparse or group-sparse component by solving a convex program. In a recent paper, [19] showed that in a noisy setting where exact recovery of the two components is impossible, it is still possible to achieve good estimation error under a comparatively mild assumption. In particular, they formulated a general measure for the radius of nonidentifiability of the problem under consideration and established a non-asymptotic upper bound on the estimation error $\|\hat{L}-L^{*}\|^{2}_{F}+\|\hat{R}-R^{*}\|^{2}_{F}$, which depend on this radius. The key idea is to allow for simultaneously sparse (or group-sparse) and low-rank matrices in the model, and control for the error introduced. We refer the readers to the above paper for a more detailed discussion on this notion of nonidentifiability. In this work, the low-rank plus sparse or group-sparse decomposition problem under restrictions on the radius of nonidentifiability takes the form

[TABLE]

[TABLE]

Here $\Omega=\left\{L\in\mathbb{R}^{p\times p}:\|L\|_{\max}\leq\frac{\alpha}{p}\right\}$ (for sparse) or $\left\{L\in\mathbb{R}^{p\times p}:\|L\|_{2,\max}\leq\frac{\beta}{\sqrt{K}}\right\}$ (for group sparse), $\|\cdot\|_{\diamond}$ represents $\|\cdot\|_{1}$ or $\|\cdot\|_{2,1}$ depending on sparsity or group sparsity of $R$, and $\lambda_{N}$ and $\mu_{N}$ are non-negative tuning parameters controlling the regularizations of low-rank and sparse/group-sparse parts. The parameters $\alpha$ and $\beta$ control for the degree of nonidentifiability of the matrices allowed in the model class. For instance, larger values of $\alpha$ provide sparser estimates of $S$ and allow simultaneously sparse and low-rank components to be absorbed in $\hat{L}$. A smaller value of $\alpha$, on the other hand, tends to produce a matrix $L$ with smaller and pushes the simultaneously low-rank and sparse components to be absorbed in $\hat{S}$. In practice, we recommend choosing $\alpha$ and $\beta$ in the range $[1,p]$ and $[1,K]$, respectively. The issue of selecting them robustly in practice is discussed in Section 5.

Remark. On certain occasions, it may be useful to have both sparse m and group-sparse structures in the model, in addition to the low rank structure. We then have $R^{*}=S^{*}+G^{*}$ in (18) with $\|R\|_{\diamond}=\|S\|_{1}+\frac{\nu_{N}}{\mu_{N}}\|G\|_{2,1}$. However, to guarantee the * simultaneous identifiability* of the sparse and group-sparse components from the low-rank component, stronger conditions need to be imposed on $L$; namely, $\Omega=\left\{L\in\mathbb{R}^{p\times p}:\|L\|_{\max}\leq\frac{\alpha}{p}~{}~{}\&~{}~{}\|L\|_{2,\max}\leq\frac{\beta}{\sqrt{K}}\right\}$.

Remark. Note that the estimated VAR model is not guaranteed to be stable, although the error bound analysis in section 3 ensures its stability with high probability as long as the sample size is large enough and the true generative model is stable. For network reconstruction and visualization purposes, stability of the estimated VAR is not strictly required. However, enforcing stability is essential for forecasting purposes. When there is a small deviation of the estimated model from stability (e.g the spectral radius of the estimated $\hat{B}$ is a little over $1$), stability can be ensured through a post-processing step of shrinking the moduli of eigenvalues of $\hat{B}$ below $1$ while keeping its eigenvectors unchanged. This type of projection argument is common in covariance and correlation matrix estimation with missing data for ensuring positive definiteness of the estimates [42]. However, in case of moderate to large deviation from stability, a closer look at the individual time series is recommended to re-assess the validity of the VAR formulation. For example, in macroeconomics, it is customary to use suitable transformations of the component time series to ensure that each of the individual time series and the resulting VAR model is stable, as opposed to modeling the individual and the joint time series (without transformation) as unit root and co-integrating processes. For instance, see [41, 40] for specific recommendations on useful transformations for macroeconomic time series.

3 Theoretical Properties

Next, we derive non-asymptotic upper bounds on the estimation errors of the low-rank plus structured sparse components of the transition matrix $B$. The main result shows that consistent estimation is possible with a sample size of the order $N\asymp p\,\mathcal{M}^{2}(f_{X})/\EuFrak{m}^{2}(f_{X})$, as long as the process $\{X^{t}\}$ is stable and the radius of nonidentifiability, as measured by $\|L^{*}\|_{\max}$ and/or $\|L^{*}\|_{2,\max}$, is small in an appropriate sense detailed next.

To establish the results, we first consider fixed realizations of $\mathcal{X}$ and $E$ and impose the following assumptions:

1. Restricted Strong Convexity (RSC): There exist $\zeta>0$ and $\tau_{N}>0$ such that

[TABLE]

where $\Phi(\Delta)=\underset{L+R=\Delta}{\text{inf}}\{\lambda_{N}\|L\|_{*}+\mu_{N}\|R\|_{\diamond}\}$, and

2. Deviation Conditions: There exist a deterministic function $\phi(B,\Sigma_{\epsilon})$ of the model parameters $B$ and $\Sigma_{\epsilon}$ such that

[TABLE]

[TABLE]

[TABLE]

where $m$ is the size of the largest group $\max_{1\leq k\leq K}card(G_{k})$.

Later on, we show that assumptions 1) and 2) are indeed satisfied with high probability when the data are generated from the model (2).

Next, we present the non-asymptotic upper bounds on the estimation errors of the low-rank plus structured sparse components, respectively.

Proposition 1**.**

(a) Suppose that the matrix $L^{*}$ has rank at most $r$, while the matrix $S^{*}$ has at most $s$ nonzero entries. Then, for any $\lambda_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{2}$ and $\mu_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{\max}+\frac{4\zeta\alpha}{p}$, any solution $(\hat{L},\hat{S})$ of (18) satisfies

[TABLE]

(b) Suppose that the matrix $L^{*}$ has rank at most $r$, while the matrix $G^{*}$ has at most $g$ non-zero groups. Then, for any $\lambda_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{2}$ and $\mu_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{2,\max}+\frac{4\zeta\beta}{\sqrt{K}}$, any solution $(\hat{L},\hat{G})$ of (18) satisfies

[TABLE]

Remark. It should be noted that if each group in $G^{*}$ has only one element, then we have $K=p^{2}$ and $g$ non-zero entries. For such cases, part (b) of Proposition 1 becomes identical to part (a).

As a byproduct, we also give the estimation error bound of the transition matrix which can be characterized by the sparse plus group-sparse and the low-rank plus sparse and group-sparse components, respectively, under the assumption that the strength of the connections in the group-sparse component $G$ is weak; i.e. $G\in\Psi$ with $\Psi=\left\{G\in\mathbb{R}^{p\times p}:\|G\|_{\max}\leq\frac{\gamma}{p}\right\}$, where $\gamma\in[1,p]$.

Proposition 2**.**

(a) Suppose that the matrix $S^{*}$ has at most $s$ nonzero entries, while the matrix $G^{*}$ has at most $g$ non-zero groups. Then, for any $\mu_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{\max}+\frac{4\zeta\gamma}{p}$ and $\nu_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{2,\max}$, any solution $(\hat{S},\hat{G})$ of (18) satisfies

[TABLE]

(b) Suppose that the matrix $L$ has rank at most $r$, while the matrix $S$ has at most $s$ nonzero entries and the matrix $G$ has at most $g$ non-zero groups. Then, for any $\lambda_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{2}$, $\mu_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{\max}+\frac{4\zeta\alpha}{p}+\frac{4\zeta\gamma}{p}$, and $\nu_{N}\geq 4\|\mathcal{X}^{{}^{\prime}}E\|_{2,\max}+\frac{4\zeta\beta}{\sqrt{K}}$, any solution $(\hat{L},\hat{S},\hat{G})$ of (18) satisfies

[TABLE]

See Appendix A for the detailed proof of Propositions 1 and 2.

Note that the objective in Proposition 2 is not the accurate recovery of the $S^{*}$ and $G^{*}$ components separately. The latter can be in principle achieved, if one sets $\gamma$ to a very small value.

In order to obtain meaningful results in the context of our problem, we need upper bounds on $\|\mathcal{X}^{\prime}E\|_{2}$, $\|\mathcal{X}^{\prime}E\|_{\max}$ and $\|\mathcal{X}^{\prime}E\|_{2,\max}$ and a lower bound on $\Lambda_{\min}(\mathcal{X}^{\prime}\mathcal{X})$ that hold with high probability. For the case of independent and identically distributed data, such high-probability deviation bounds are established in [19]. However, for time series data all entries of the $\mathcal{X}$ matrix are dependent on each other, and hence it is a non-trivial technical task to establish such deviation bounds. A key technical contribution of this work is to derive these deviation bounds, which lead to meaningful analysis for VAR models. The results rely on the measure of stability defined in (2) and an analysis of the joint spectrum of $\{X^{t-1}\}$ and $\{\epsilon^{t}\}$ undertaken next.

Proposition 3**.**

Consider a random realization of $\{X^{0},\ldots,X^{T}\}$ generated according to a stable VAR(1) process (2) and form the autoregressive design (17). Define

[TABLE]

Then, there exist universal positive constants $c_{i}>0$ such that

for $N\succsim p$,

[TABLE]

and for any $N\succsim\log p$,

[TABLE]

and for $N\succsim m\log p$,

[TABLE] 2. 2.

for $N\succsim p\mathcal{M}^{2}(f_{X})/\EuFrak{m}^{2}(f_{X})$,

[TABLE]

See Appendix B for the detailed proof of Proposition 3.

Using the above deviation bounds in the non-asymptotic errors of Propositions 1, we obtain the final result for approximate recovery of the low-rank and the structured sparse components using nuclear and $\ell_{1}/\ell_{2,1}$ norm relaxations, as we show next.

Proposition 4**.**

Consider the setup of Proposition 3. There exist universal positive constants $c_{i}>0$ such that for $N\succsim p\mathcal{M}^{2}(f_{X})/\EuFrak{m}^{2}(f_{X})$, and $\|L^{*}\|_{\max}\leq\alpha/p$, any solution $(\hat{L},\hat{S})$ of the program (18) satisfies, with probability at least $1-c_{1}\exp[-c_{2}\log p]$,

[TABLE]

Remark: The error bound presented in the above proposition consists of two key terms. The first term is the error of estimation emanating from randomness in the data and limited sample capacity. For a given model, this error goes to zero as the sample size increases. The second term represents the error due to the unidentifiability of the problem. This is more fundamental to the structure of the true low-rank and structured sparse components, and depends only on the model parameters and does not vanish, even with infinite sample size.

Further, the estimation error is a product of two terms - the second term $(rp+s\log p)/N$ involves the dimensionality parameters and matches the parametric convergence rate for independent observations. The effect of dependence in the data is captured through the first part of the term: $\frac{c_{0}\phi^{2}(B,\Sigma_{\epsilon})\mu_{\max}^{2}(\mathcal{B})}{\Lambda_{\min}^{2}(\Sigma_{\epsilon})}$. As discussed in [8], this term is larger when the spectral density is more spiky, indicating a stronger temporal and cross-sectional dependence in the data.

Proposition 5**.**

Consider the setup of Proposition 3. There exist universal positive constants $c_{i}>0$ such that for $N\succsim p\mathcal{M}^{2}(f_{X})/\EuFrak{m}^{2}(f_{X})$, for any $G^{0}$ with $\|L^{*}\|_{2,\max}\leq\beta/\sqrt{K}$, any solution $(\hat{L},\hat{G})$ of the program (18) satisfies, with probability at least $1-c_{1}\exp[-c_{2}\log p]$,

[TABLE]

See Appendix B for the detailed proof of Proposition 4 and 5.

Remark. Based on Proposition 5, similar conclusions can be obtained as that for the low rank plus sparse case.

4 Computational Algorithm and its Convergence Properties

Next, we introduce a fast algorithm for estimating the transition matrix $B$ from data. For ease of presentation and to convey the key ideas clearly, we first present the algorithm for $B$ representing a single structure (e.g. only low rank, or only group sparse, or only sparse), and in addition establish its convergence properties. Subsequently, we modify the algorithm to handle the composite structures considered in this paper and also establish its convergence.

The fast network structure learning (FNSL) Algorithm 1 is described next. A safeguarded BB initial value is selected, as the initial choice of the nominal step $\eta_{i}$, i.e.

[TABLE]

For notational convenience, the penalty term for estimating the transition matrix $B$ is denoted by $P_{B}(B,\lambda)$, where $\lambda>0$ represents the tuning parameter.

The specific $B_{i+1}$ update depends on the employed penalty term; for an $\ell_{1}$ penalty inducing sparsity, it corresponds to soft-thresholding [27], for a group sparse penalty to group soft-thresholding [20], while for a nuclear norm penalty to singular value thresholding [28].

It can also been seen in Algorithm 1, that for $\alpha_{i}\equiv 1$ for $\forall i\geq 1$, then $B_{i}^{md}=B_{i}$ and $B_{i+1}^{ag}=B_{i+1}$, which leads to the traditional gradient descent algorithm. Indeed, Algorithm 1 is obtained by incorporating an efficient backtracking strategy into the accelerated multi-step scheme by [29, 30]. It provides a different way to look for a larger stepsize by employing a relaxed line search condition, instead of searching for the gradient Lipschitz constant of the data fidelity term. Steps 1 and 2 constitute the backtracking ones. Both $B_{i}^{md}$ and $B_{i}^{ag}$ are linear combinations of all past iterations of $B_{i}$, but based on different weights as can be seen from their updates, i.e. $B_{i}^{ag}=(1-\alpha_{i-1})B_{i-1}^{ag}+\alpha_{i-1}B_{i}$ and $B_{i}^{md}=(1-\alpha_{i})B_{i}^{ag}+\alpha_{i}B_{i}$. Here ‘ag’ simply denotes ‘aggregate’. The data fidelity term in the cost function is linearized at $B_{i}^{md}$. Since it is used after we have obtained $B_{i}^{ag}$ and before we obtain $B_{i+1}$, we use ‘md’ to denote a ‘middle’ update. Further, $\Gamma_{i}$ and $Q_{i}$ are the parameters most closely related to our line search conditions. Intuitively, if we set $\Gamma_{i}\geq 0$, we are certainly able to guarantee the accelerated rate of convergence. However, this will render the stepsize smaller. Fortunately, our convergence analysis enables us to relax this condition by utilizing the summing up procedure, with $Q_{i}$ corresponding to the part of the sum of $\Gamma_{i}$. Thus, we can impose a relaxed termination condition on $Q_{i}$ (see step 2 of Algorithm 1) without impacting the rate of convergence while being able to obtain a more aggressive stepsize. In fact, parameter $C$ in Step 2 plays an important role, i.e. the number of the trial steps can be reduced significantly when a relatively larger $C$ is selected. However, the value of $C$ can not be too large either, since it might impair the convergence rate in terms of the objective function value.

The convergence rate of the proposed algorithm 1 is established next.

Proposition 6**.**

Let $\{B_{k+1}^{ag}\}$ be generated by Algorithm 1. Then, for any $k\geq 1$

[TABLE]

where $\tilde{C}$ is a finite positive number independent of $k$.

See Appendix C for the detailed proof. It should be noted that the convergence rate of $\{B_{k+1}\}$ is still an open problem, which needs to be addressed further in future, considering its good performance for some cases.

Next, we enhance the algorithm for solving (18) in the general case. The accelerated convergence rate can be obtained by following the proof for Proposition 6. Similarly to the case in Algorithm 1, the initial trial step of $\eta_{i}$ is a safeguarded BB choice

[TABLE]

The update of the $L$ component is based on singular value thresholding, while that of the $R$ component on (group) soft-thresholding. Note that the most expensive computational operation corresponds to the singular value decomposition (SVD) when updating $L_{i+1}$. As mentioned earlier, the proposed algorithm is able to look for larger magnitude step sizes by conducting fewer number of line searches, due to employing more relaxed line search conditions. Actually, this is an important improvement considering the computational cost of SVD. Indeed, the efficiency of the proposed algorithm can be enhanced further if we employ the truncated SVD [32] instead of the full SVD.

The convergence rate of the proposed algorithm 2 (given in supplement due to limited space) is established next.

Proposition 7**.**

Let $(L_{k+1}^{ag},R_{k+1}^{ag})$ be a sequence of updates generated by Algorithm 2. Then, for any $k\geq 1$

[TABLE]

where $\tilde{C}$ is a finite positive number independent of $k$.

Proposition 7 is a direct extension of Proposition 6 and can be obtained by following the roadmap of the proof for the former result.

5 Performance Evaluation

Next, we present experimental results on both synthetic and real data. Specifically, the first two experiments focus on large-scale network learning with single penalty term to show the efficiency and effectiveness of the proposed algorithms, while the remaining ones assess the accuracy of recovering low rank plus structured sparse transition matrices $B$.

5.1 Performance metrics and experimental settings

We introduce the performance metrices used in the numerical work. For network estimation, we use the true positive rate (TPR) and false alarm rate (FAR) defined as:

• •

$\text{TPR}:=\frac{\sharp\{\hat{b}_{ij}\neq 0\ \text{and}\ b_{ij}\neq 0\}}{\sharp\{b_{ij}\neq 0\}}$

• •

$\text{FAR}:=\frac{\sharp\{b_{ij}=0\ \text{and}\ \hat{b}_{ij}\neq 0\}}{\sharp\{b_{ij}=0\}}$

where $b_{ij}$ and $\hat{b}_{ij}$ are the correspnding elemnets in $B$ and $\hat{B}$, respectively. The estimation error (EE) and out-of-sample prediction error (PE) are defined as

• •

$\text{EE}:=\frac{\|\hat{B}-B\|_{F}}{\|B\|_{F}}$

• •

$\text{PE}:=\|\hat{\mathcal{Y}}-\mathcal{Y}\|^{2}_{F}/\|\mathcal{Y}\|^{2}_{F}$

To select the optimal value of the tuning parameters, we combine the three (or two or one, respectively) -dimensional grid search method with the AIC/BIC/forward cross-validation criteria. We will specify the criterion on a case by case basis for the following experiments. In examples B and C, the tuning parameter $\lambda$ is selected by the AIC criterion. A grid of 100 values in the interval $[0,\|\mathcal{X}^{\prime}\mathcal{Y}\|_{\max}]$ is used for $\lambda$. In examples D, E and F, we utilize a two/three-dimensional grid search to select the optimal values of $\lambda_{N}$, $\mu_{N}$ and/or $\nu_{N}$ as that for $\lambda$. For the experiments employing synthetic data, the tuning parameters are selected by assuming the rank of the true low-rank transition matrix and/or the non-zero group-sparse components of the true group-sparse transition matrix are known. We will specify the forward cross-validation procedure for the real data case in example G.

For all the experiments, the parameters used in the proposed algorithms are depicted in table 1. Also we set $\Sigma_{\epsilon}=\epsilon^{2}I$ and $\epsilon^{2}=1$. We rescale the entries of $B$ to ensure stability of the process (the spectral radius $\rho\in(0.45,0.95)$. All the results are based on 50 replications. Finally, all algorithms are run in the MATLAB R2015a environment on a PC equipped with 12GB memory.

5.2 Large-scale sparse network learning

We start by comparing the performance of the proposed algorithm 1 with FISTA with line search [33] to solve problem (18) with a sparse transition matrix.

We consider three different VAR(1) models with $p=800,\,900$ and $1000$ variables. For each of these models, we generate $N=1000,1500$, and $2000$ observations from a Gaussian VAR(1) process (2). The $p\times p$ transition matrix $B$ with sparsity is generated in the following way. First, the topology is generated from a directed random graph $G(p,\xi)$, where the edge from one node to another node occurs independently with probability $\xi=10/p$. Then, the strength of the edges is generated independently from a Gaussian distribution. This process is repeated until we obtain a transition matrix $B$ with a desired spectral radius $\rho$. We compare TPR, FAR, EE, and computational time (denoted by $T$).

Table 2 shows the experimental results for sparse network structure with different network size $p$ and sample size $N$. It can be seen that the proposed algorithm performs similarly to FISTA in terms of TPR, FAR, and estimation error. To show the efficiency of the proposed algorithm, we also compare the computational time in seconds in terms of the convergence of the objective function value. Clearly, the proposed algorithm outperforms FISTA in efficiency, especially when the network and sample size become larger. This is mainly due to the relaxed line search scheme, as previously discussed. To further support our claim, we also show the graphs of the decreasing objective function value vs. CPU, see Figure 1 when $p=1000$ and $N=2000$.

Table 3 shows comparisons between a variant of FISTA and FNSL in terms of objective function value, CPU time in seconds, and the number of matrix products for a network of size $p=1000$ with sample size 1000, 1500 and 2000, respectively. For each data set, FISTA needs around 30 iterations to reach convergence and the total number of line searches is 3 for all iterations, while FNSL needs no more than 20 iterations and the total number of line search is no more than 4 for all iterations. This illustrates the computational savings of the proposed algorithm.

5.3 Network learning with a low-rank transition matrix

To further show the efficiency of the proposed algorithms, we compare the performance of a variant of FISTA [34] and FNSL on estimating low-rank transition matrices.

We consider three different VAR(1) models with $p=200,\,300$ and $400$ variables. For each of these models, we generate $N=400,1200$, and $2000$ observations from a Gaussian VAR(1) process (2). The $p\times p$ low-rank transition matrix $B$ is generated with rank $\lfloor p/25\rfloor+1$. Subsequently, we rescale the entries of $B$ to ensure the spectral radius $\rho$ lies in $(0.45,0.95)$. We compare the rank of the estimated transition matrix, denoted by $\hat{r}$, EE, and computational time $T$.

Table 4 shows the experimental results for low-rank network structure with different network and sample size. Both FISTA and the proposed algorithm achieve good recovery of the transition matrix $B$ with the correct rank and they have similar performance in terms of estimation error. Clearly, the proposed algorithm outperforms FISTA in efficiency for this case as well. The graphs of the decreasing objective function value vs. CPU are depicted in Figure 2.

The first two experiments focused on computational efficiency of the proposed algorithms, while retaining good network estimation properties. Next, we demonstrate their accuracy for learning structured sparse networks.

5.4 Sparse plus low-rank network learning

Next, we investigate estimation of sparse plus low-rank transition matrices and compare it to ordinary least square (OLS) and lasso estimates.

We consider three different VAR(1) models with $p=50,\,75$ and $100$ variables. For each of these models, we generate $N=100$ and $200$ observations from the model defined in (2) where $B$ can be decomposed into a low-rank matrix $L$ of rank $\lfloor p/25\rfloor+1$ and a sparse matrix $S$ with $2-4\%$ non-zero entries. We rescale the entries of $B$ to ensure stability of the process (the spectral radius is set to $\rho(B)=0.7$). We compare the estimation and out-of-sample prediction errors. The number of out of samples is set to 10.

First, we study the influence of $\alpha$ in (18) on this learning problem with $p=50$ and $N=200$. From Table 5, that a smaller $\alpha$ parameter leads to markedly improved identification of all the true nonzero entries in the sparse component, which consequently leads to better separating the sparse component $S$ from the low-rank component $L$.

The corresponding estimation errors are reported in Table 6. In all three settings, we find that the low-rank plus sparse VAR estimates outperform the estimates using ordinary least-squares (OLS) and lasso, as expected. We observe that as the ratio of $N/p$ increases, OLS may produce lower estimation error than lasso, even though the OLS model is not interpretable for this case. Also, we observe that the estimation errors of all three methods decrease with increasing sample sizes as expected and predicted by theory. Further, we illustrate how the squared Frobenius norm error in Proposition 4 of the VAR model with low rank plus sparse transition matrix scales with the sample size $N$ and dimension $p$, when the rank $r$ of $B$ is fixed. The network size $p$ is set to $50,100,150$ and $200$, respectively, while the rank $r$ is fixed to be 2 for all $p$. Sparsity $s$ and $\varepsilon$ are defined similarly as above, while the sample size is set to $N\in(150,5500)$. The squared Frobenius norm error of estimation given by $\|S-\hat{S}\|_{\text{F}}^{2}+\|L-\hat{L}\|_{\text{F}}^{2}$ is depicted in the left panel of Figure 3, which displays the errors for different values of $p$, plotted against the sample size $N$. As predicted by our theoretical result, the error is larger for larger $p$. The right panel of Figure 3 displays the error against the rescaled sample size $N/(s\,log(p)+rp)$. it can be seen that the corresponding error curves for different values of $p$ align well, which is consistent with the estimation error obtained in Proposition 4.

In addition to its improved estimation and prediction performance, the low-rank plus sparse modeling strategy aids in recovering the underlying Granger causal network after accounting for the latent structure. In Figures 4, we demonstrate this using a VAR(1) model with $p=50$ and $n=200$. The top panel of the Figures 4 displays the true transition matrix $B$, its low-rank component $L$ and the structure of its sparse component $S$. The bottom panel of the Figures 4 displays the structure of the Granger causal networks estimated by the method of Lasso and the low-rank plus sparse modeling strategy. As predicted by the theory, it can be that the lasso estimate of the Granger causal network selects many false positives due to its failure to account for the latent structure. On the other hand, the sparse component $S$ provides an estimate exhibiting significantly fewer false positives entries.

It is interesting to note that the estimation performance of the regularized estimates in low-rank plus sparse VAR models is worse than the performance of lasso in sparse VAR models of similar dimension [8], even for the same sample sizes. This is in line with the error bounds presented in Proposition 4. The estimation error in low-rank plus sparse models is of the order of $O(rp+s\log p)/N$, while the error of lasso in sparse VAR models scales at a faster rate of $O(s\log p/N)$. Further note that a $s$-sparse VAR requires estimating $s$ parameters in $S$, while the presence of $r$ factors introduces an additional $rp$ parameters in the loading matrix $\Lambda$.

5.5 Sparse plus group-sparse plus low-rank network learning

Finally, we conduct numerical experiments to assess the performance of low-rank plus sparse plus group-sparse modeling in VAR analysis and compare it to the performance of sparse plus group-sparse and low-rank plus sparse estimates.

We consider three different VAR(1) models with $p=50,\,100$ and $150$ variables. For each of these models, we generate $N=200$ and $300$ observations from the Gaussian VAR(1) process defined in (2), where $B$ can be decomposed into a low-rank matrix $L$ of rank $\lfloor p/25\rfloor+1$, a sparse matrix $S$ with $2-4\%$ non-zero entries, and a group-sparse matrix $G$ with each column corresponding to a different group for a total of $p$ groups. We rescale the entries of $B$ to ensure stability of the process (the spectral radius is set to $\rho(B)=0.7$) and compare the estimation and out-of-sample prediction errors, with the number of out-samples set to 10.

The corresponding estimation errors are reported in Table 7. In those three settings, we find that the low-rank plus sparse plus group-sparse VAR estimates performs only slightly better than low-rank plus sparse VAR estimates. One of the reasons lie in that the ability of the identification will degrade as more structures are involved. The other one is that multiple-times shrinkage for the multiple structures lead to severe bias estimation. Even though the group structures can be recovered completely, some non-zero elements in sparse component vanished. An ad-hoc way to improve the performance for this case is to combine these two methods together. However, both methods outperform the estimates using sparse plus group-sparse VAR. We also observe that, as the ratio of $N/p$ increases, the estimation errors of all three methods decrease with increasing sample sizes as expected and predicted by theory.

5.6 Structured network learning of asset pricing data

Finally, we employ the proposed framework to learn Granger causal networks of asset pricing data obtained from the University of Chicago’s Center for Research in Security Prices (CRSP) and retrieved from the Wharton Research Data Service (WRDS). Specifically, we examine the network structure of realized volatilities of financial institutions representing banks (BA), primary broker/dealers (PB) and insurance companies (INS). The analysis was performed across the following time periods: September 2002 - August 2005, September 2006 - August 2008 and September 2010 - August 2012 that correspond to instances before the financial crisis of 2008 (pre-crisis period), the build-up and apex of the crisis and the post-crisis period, respectively. For each period, we collected data on 75 firms with 25 companies in each of the three categories - BA, PB and INS based on the size of the average market capitalization of each firm, but dropped a few due to duplicate/missing observations. The final form of the variables used are based on the $\log$ transformation of the difference between the highest and lowest stock price during a day that acts as a proxy for realized volatility and subsequently detrending it.

To select the tuning parameters, we employ the following forward cross-validation procedure: (1) We use a time window of length $W$, the available number of time points in the data. Then, starting from time $t$, we use the most recent $W$ observations to estimate $B$ and denote the transition matrix estimate by $\hat{B}_{t}$. We use the next $W^{\prime}$ observations (right after $W$ observations) to validate $\hat{B}_{t}$. (2) We select the optimal tuning parameters $(\lambda_{N}^{{}^{\prime}},\mu_{N}^{{}^{\prime}})$ so that

[TABLE]

where $Err(\hat{B}_{t})=\|\mathcal{Y}_{W^{\prime}}-\mathcal{X}_{W^{\prime}}\hat{B}_{t}\|_{F}^{2}$ and $i=0,...$

In our analysis, we set $W=500$ and $W^{\prime}=50$. Further, to separate the sparse component from the low-rank component as much as possible, we set $\alpha=p/10$. The learned Granger causal network structures (the sparse component $\hat{S}$) estimated by a sparse plus low-rank model are depicted in Figure 5. It can be seen that even in the presence of a low-rank component, the sparse component exhibits a certain density (about 5% in the pre-crisis period, rising to 10% during the crisis and dropping down to about 6% in the post-crisis period). This increased connectivity during the crisis period has been observed in Granger causal networks for log-returns as well [4]. In the Supplement, we also provide the Granger causal network structures estimated by only assuming sparsity of the transition matrix $B$ (see Supplement Figure 5). A similar increased connectivity pattern is observed during the crisis period. Also note that after accounting for the low-rank component, the estimated sparse component is significantly more sparse than that estimated by a lasso approach, thus enabling us to better examine specific firms that are key drivers in the volatility network.

6 Discussion

Our modeling and technical developments were based on a VAR(1) model. However, it can be generalized to VAR(d) models in different ways, depending on the context of the application. One possible formulation where the low rank component stays the same across lags can be expressed as

[TABLE]

This model can be posed as a 1-order (L+S) model on the concatenated process

$\tilde{X}_{t}=\left[X_{t}^{\top},X_{t-1}^{\top},\ldots,X_{t-d+1}^{\top}\right]^{\top}$ using the standard transformation [7]:

[TABLE]

Since the matrices $\tilde{L}$ and $\tilde{S}$ are respectively low-rank and sparse, our proposed algorithm can be used and the error bounds will be applicable. Further, to maintain the special structure of these new matrices, additional constraints can be imposed to the posited objective function or the final estimates can be projected to this space.

Finally, motivated by [37, 38], it would be of interest to consider analogous developments in the frequency domain and address identifiability issues, as well as establish finite sample error bounds.

Acknowledgment

This work was supported by a UF Informatics Institute Fellowship to XL, NSF grant DMS-1812128 to SB, and NSF grants IIS 1632730, CCF 1540093, DMS 1545277 to GM.

Appendix A Estimation Error Bounds

The proof of Proposition 1 and part (a) of Proposition 2 can be easily obtained by the following proof for part (b) of Corollary 2.

Proof of part (b) of Corollary 2.

By the optimality of $(\hat{L},\hat{S},\hat{G})$ and the feasibility of $(L^{*},S^{*},G^{*})$, we have

[TABLE]

By setting $\hat{\Delta}^{L}=\hat{L}-L^{*}$, $\hat{\Delta}^{S}=\hat{S}-S^{*}$, and $\hat{\Delta}^{G}=\hat{G}-G^{*}$ and combining with $\mathcal{Y}=\mathcal{X}(L^{*}+S^{*}+G^{*})+E$, we have

[TABLE]

By Lemma 1 in [19] and lemma 2.3 in [35], we obtain

[TABLE]

where the matrices $(A,B)\in\left\{(\bar{A},\bar{B}):\bar{A}\bar{B}^{\prime}=0~{}~{}\&~{}~{}\bar{A}^{\prime}\bar{B}=0\right\}$, $(M,M^{\bot})$ and $(N,N^{\bot})$ denote an arbitrary subspace pair for which $\|S\|_{1}$ and $\|G\|_{2,1}$ are decomposable, respectively. Since

[TABLE]

Substituting (28) into (27) and recalling conditions for $\lambda_{N}$, $\mu_{N}$ and $\nu_{N}$, we have

[TABLE]

By the RSC condition, the constraints on $L$ and $G$ in (18), and the definition of $\mu_{N}$ and $\nu_{N}$, we have

[TABLE]

Inserting the above inequality into (29), we have

[TABLE]

By the compatibility constant in [19], we have

[TABLE]

By our assumptions, we have

[TABLE]

Combining with the inequality $\|\hat{\Delta}^{L}\|_{F}^{2}+\|\hat{\Delta}^{S}+\hat{\Delta}^{G}\|_{F}^{2}\leq 2(\|\hat{\Delta}^{L}\|_{F}^{2}+\|\hat{\Delta}^{S}\|_{F}^{2}+\|\hat{\Delta}^{G}\|_{F}^{2})$, we conclude part (b) of Corollary 2. ∎

Appendix B Deviation Bounds

Proof of Proposition 3.

We want to find upper bounds on $\|\mathcal{X}^{\prime}E/N\|_{\max}$, $\|\mathcal{X}^{\prime}E/N\|$ and $\|\mathcal{X}^{\prime}E/N\|_{2,\max}$ that hold with high probability. Note that such an upper bound for $\|\mathcal{X}^{\prime}E/N\|_{\max}$ has already been derived in [8]. Here we adopt a different technique that takes a unified approach to provide upper bounds on both quantities. To this end, note that the two norms have the following representations

[TABLE]

For any given $u,v\in\mathbb{S}^{p-1}$, we first provide a bound on $u^{\prime}(\mathcal{X}^{\prime}E/N)v$.

Using Proposition 2.3 of [8], we obtain

[TABLE]

for any $u,v\in\mathbb{S}^{p-1}$ and any $\eta>0$.

To derive the deviation bound on $\|\mathcal{X}^{\prime}E/N\|_{\max}$, we simply take a union bound over the $p^{2}$ possible choices of $u,v\in\{e_{1},e_{2},\ldots,e_{p}\}$. This leads to

[TABLE]

Since $N\succsim p$, we can set $\eta=\sqrt{(2+c_{1})\log p/cN}$ so that $\eta<1$ (i.e., $\eta^{2}<\eta$) will be satisfied for large enough $N$. This implies that

[TABLE]

for some universal constants $c_{i}>0$.

Similarly, for any group $G_{i}$ of size $m_{i}$, we have

[TABLE]

Taking a union bound over $K$ non-overlapping groups $G_{i}$ leads to

[TABLE]

where $m=\max_{i=1,\ldots,K}m_{i}$. As before, setting $\eta=\sqrt{\log p/N}$ implies

[TABLE]

for some $c_{i}>0$.

To derive the deviation bound on the spectral norm, we discretize the unit ball $S^{p-1}$ using an $\epsilon$-net $\mathcal{N}$ of cardinality at most $(1+2/\epsilon)^{p}$. An argument along the line of Supplementary Lemma F.2 in [8] then shows that for a small enough $\epsilon>0$,

[TABLE]

for some constant $\kappa>1$, possibly dependent on $\epsilon$. As before, taking a union bound over the $(1+2/\epsilon)^{2p}$ choices of $u$ and $v$, we get

[TABLE]

Since $N\succsim p$, choosing $\eta=\sqrt{(c_{1}+2\log(1+2/\epsilon))p/cN}$ ensures $\eta<1$ for large enough $N$. Setting $\eta$ as above concludes the proof. 2. 2.

We want to obtain a lower bound on the minimum eigenvalue of $\mathcal{X}^{\prime}\mathcal{X}/N$ that holds with high probability.

Since $\Lambda_{\min}\left(\mathcal{X}^{\prime}\mathcal{X}/N\right)=\inf_{v\in\mathbb{S}^{p-1}}v^{\prime}(\mathcal{X}^{\prime}\mathcal{X}/N)v$, we start with the single deviation bound of Proposition 2.3 in [8],

[TABLE]

for any $v\in\mathbb{S}^{p-1}$ and $\eta>0$.

The next step is to extend this single deviation bound uniformly on the set $\mathbb{S}^{p-1}$. As in the proof of part 1, we construct a $\epsilon$-net of cardinality at most $(1+2/\epsilon)^{p}$ and approximate the quadratic form using its values on the net. This yields the following deviation bound

[TABLE]

for some constant $\kappa>1$. Seting $\eta=\EuFrak{m}(f_{X})/4\kappa\pi\mathcal{M}(f_{X})<1$ and noting that $N\succsim\mathcal{M}^{2}(f_{X})/\EuFrak{m}^{2}(f_{X})p$, we conclude

[TABLE]

The result follows from the lower bound on $\EuFrak{m}(f_{X})$ presented in (6) and the fact that $v^{\prime}\Gamma_{X}(0)v\geq\EuFrak{m}(f_{X})$ for all $v\in\mathbb{S}^{p-1}$.

∎

Proof of Proposition 4 and 5.

Clearly, setting $\zeta$ to the lower bound on $\EuFrak{m}(f_{X})$ as in (6) satisfies the RSC. Combining the estimates of Proposition 1 and 3 leads to proposition 4 and 5 after simple algebraic computation. ∎

Appendix C Convergence Analysis

In the following proof, we denote $\frac{1}{2}\left\|\mathcal{Y}-\mathcal{X}B\right\|^{2}_{F}$ and the regularization term by $H(B)$ and $P_{B}(B,\lambda)$, respectively.

Proof of Proposition 6.

By the differentiability of $H$, we have

[TABLE]

Then, by the definition of $H(B)$ and $B_{i+1}^{ag}$, and the relationship $B_{i+1}^{ag}-B_{i}^{md}=\alpha_{i}(B_{i+1}-B_{i})$, we have

[TABLE]

By the convexity of $H(B)$ and (35), we have

[TABLE]

Subtracting $l(B)$ from both sides of (36) and rearranging some terms, we have

[TABLE]

where $\xi\in\partial P_{B}(B_{i+1},\lambda)$. On the other hand, by the first-order optimality conditions for the sequence $B_{i+1}$ in Algorithm 1, we have

[TABLE]

Combining (37) and (38), we obtain

[TABLE]

where we used the relationship $2\langle a-b,a-c\rangle=-\|b-c\|^{2}+\|a-c\|^{2}+\|a-b\|^{2}$ and the definition of $B_{i+1}^{e}$.

Dividing both sides of (39) by $\alpha_{i}\eta_{i}$, we have

[TABLE]

Adding $\frac{(\beta_{i}Q_{i}+\Gamma_{i})}{2}$ to both sides of (40), we have

[TABLE]

Since $Q_{i+1}=\beta_{i}Q_{i}+\Gamma_{i}$, $0\leq\beta_{i}\leq(1-\frac{1}{i})^{2}$, and $Q_{i}\geq-\frac{C}{(i-1)^{2}}$, we obtain

[TABLE]

Setting $B=\hat{B}$, by the relationship $\frac{1}{\alpha_{i}\eta_{i}}=\frac{1-\alpha_{i+1}}{\alpha_{i+1}\eta_{i+1}}$, and $\alpha_{1}=1$, we obtain

[TABLE]

after summing (42) from $i=1$ to $k$.

Next we show the upper bound of $\alpha_{k}\eta_{k}$. Since $\eta_{\min}\leq\eta_{0,1}$, we have $\eta_{\min}\leq\|\mathcal{X}^{T}\mathcal{X}\|_{2}$. Then, by definition of $\eta_{0,i}$, we get

[TABLE]

Denote $\sigma^{l}\eta_{0,i}$ by $\eta_{i}^{{}^{\prime}}$, where $l$ is the number of line search in Step 3 of Algorithm 1. By $\frac{1}{\alpha_{i}\eta_{i}}=\frac{1-\alpha_{i+1}}{\alpha_{i+1}\eta_{i+1}}$ and the definition of $\eta_{i}$, we have

[TABLE]

Then, by induction we can get, with $\alpha_{1}=1$,

[TABLE]

which implies

[TABLE]

where we used (45) and the definition of $\eta_{i}^{{}^{\prime}}$.

Combining (43) and (46), we obtain (22). ∎

Appendix D Minimizer of optimization problems

For convenience, we provide solutions for the minimization problems in Algorithm 1 when the penalty term has different norms.

When $P(B,\lambda)=\lambda\|B\|_{1}$,

[TABLE]

where $S(B)_{kl}=\max(0,|b_{kl}|-\frac{\lambda}{\eta_{i}})\text{sgn}(b_{kl})$ and $b_{kl}$ is the entry of matrix $B$ from $k$th row and $l$th column.

When $P(B,\lambda)=\lambda\|B\|_{2,1}$,

[TABLE]

where $GS((B)_{G_{k}})=\max(0,\frac{\lambda/\eta_{i}}{\|(B)_{G_{k}}\|_{2}})((B)_{G_{k}})$

When $P(B,\lambda)=\lambda\|B\|_{\ast}$,

[TABLE]

where $D_{\tau}(Z)=UD_{\tau}V^{\ast}$, which is singular value decomposition. $D_{\tau}=\text{diag}(\max(0,\sigma_{i}-\lambda/\eta_{i}))$ and $\{\sigma_{i}\}_{i=1}^{k}$ are the singular values of matrix $Z$.

Appendix E Algorithm 2

Due to space limitations, the detailed Algorithm 2 is given next.

Appendix F Additional Numerical Experiments

Sparse plus group-sparse network learning problem

We discuss the experimental setting for a sparse plus group-sparse transition matrix and compare the performance with methods that either assume pure sparsity (lasso) or pure group sparsity (group lasso, SGL).

We consider three different VAR(1) models with $p=50,\,100$ and $200$ variables. For each of these models, we generate $N=100$ and $200$ observations from a Gaussian VAR(1) process (see equation (9) in the main file) where $B$ can be decomposed into a sparse matrix $S$ with $5\%$ non-zero entries and a group-sparse matrix $G$ with each column corresponding to a different group (hence we have $p$ groups in $G$). We randomly select two columns (groups) consisting of two super hubs, in which the strength of the edges is generated independently from a Gaussian distribution. To be more consistent with our error bound analysis, we set $\gamma$ to $p/2$ in (18). The network topology of $S$ is generated the same way as that in subsection 5.2 except that the occurring probability of the edge from one node to another node $\xi$ is set to be $0.05$. Subsequently, the entries of the corresponding two columns in $G$ are set to be zero. Finally, we rescale the entries of $B$ so that a desired spectral radius is reached. We employ the TPR, FAR, and EE metrics in the comparisons.

Table 8 shows the experimental results for different network size and number of samples. It can be seen that utilizing an S+G model enables us to identify a larger portion of correct nonzero numbers in $B$, while achieving almost the same false alarm rate compared to lasso and SGL. Particularly, the S+G model can recover the group information perfectly while lasso and SGL miss a number of edges, as expected since they correspond to misspecified structures. Further, the S+G model exhibits the lowest estimation error amongst them. It should be noted that the advantage of the S+G model will be more evident if the strength of the edges within the groups is weaker.

In Supplement Figure 1 the true network structure of S+G, S and G with $p=50$ and $N=200$ is depicted. The recovered network structures by the S+G model are given in Supplement Figure 2 (top) and the group-sparse components are recovered perfectly. We also compare the recovered network structures by S+G, SGL and lasso models, see Supplement Figure 2 (bottom), from which we can see that S+G performs best.

Figure 6 shows the true network structure S+G, S and G with $p=50$ and $N=200$. The recovered network structures by S+G model are given in Figure 7 (top). Clearly, the group-sparse components are recovered perfectly. We also compare the recovered network structures by S+G, SGL and lasso models, see Figure 7 (bottom), from which we can see that S+G performs best.

Figures 9 shows the estimated Granger causal network using low-rank plus sparse plus group-sparse VAR estimates using a VAR(1) model with $p=50$ and $n=300$. The top panel of the Figures 9 displays the true the structure of sparse plus group-sparse components $S+G$, the structure of group-sparse component $G$, and the structure of low-rank component $L$. The bottom panel of the Figures 9 displays the structure of the Granger causal networks estimated by the method of $L+S$ and $S+G$ modeling strategy. It can be seen that the $S+G$ estimate selects many false positives due to its failure to account for the latent structure. On the other hand, the $L+S$ method provides an estimate exhibiting significantly fewer false positives entries as that by $L+S+G$.

The learned Granger causal network structures estimated by lasso are given in figure 10, which correspond to the case of pre-crisis, crisis and post-crisis, respectively. From the network structures by Lasso and S+L strategy, we can see that companies are highly connected when financial crisis is coming. After accounting for the low-rank component, the estimated sparse components by L+S strategy are usually more sparse than that by Lasso.