Linear Convergence of Gradient Methods for Estimating Structured Transition Matrices in High-dimensional Vector Autoregressive Models
Year: 2022

Abstract:
In this paper, we present non-asymptotic optimization guarantees of gradient descent methods for estimating structured transition matrices in high-dimensional vector autoregressive (VAR) models. We adopt the projected gradient descent (PGD) for single-structured transition matrices and the alternating projected gradient descent (AltPGD) for superposition-structured ones. Our analysis demonstrates that both gradient algorithms converge linearly to the statistical error even though the strong convexity of the objective function is absent under the high-dimensional settings. Moreover our result is sharp (up to a constant factor) in the sense of matching the phase transition theory of the corresponding model with independent samples. To the best of our knowledge, this analysis constitutes first nonasymptotic optimization guarantees of the linear rate for regularized estimation in high-dimensional VAR models. Numerical results are provided to support our theoretical analysis.

heading: 1 Introduction
text: Learning the network structure through high-dimensional time series data has been an important focus of research for the past decades. There are lots of application examples ranging from macroeconomic analysis [1–3] to connectivity measuring among financial firms [4], gene regularity network inference [5] and radar signals processing [6, 7]. To perform these tasks, vector autoregressive (VAR) models play a critical role in both theory and application. For example, VAR models are widely adopted to characterize the spatially and temporally colored disturbance for multichannel adaptive signal detection in [7–11].
Under the low-dimensional settings where the dimension of the transition matrix and the number of time series are relatively small, the theory of VAR models is well established, see e.g., [12]. However, lots of meaningful applications are under the high-dimensional settings where the problem dimension far exceeds the number of time series and additional structure information of parameters is required to guarantee successful recovery. For the cases where the samples are independent, both theoretical properties and practical algorithms of high-dimensional statistical problems have been studied by considerable literature during the past few years, including but not limited to [13–17]. For the correlated time series cases, the corresponding results are still developing. By assuming that the spectral norm of the transition matrix is less than 1, Loh and Wainwright consider VAR models regularized by the l1-norm in [18]. Under the double asymptotic framework, Han and Liu also consider VAR models under the similar assumption in [19]. In [20], Basu and Michailidis analyze sparse transition matrices estimation of VAR models with a milder stability assumption by
∗Corresponding author.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
introducing the spectral density. Recently, Melnyk and Banerjee extend the analysis to structured VAR models regularized by any suitable norm in [21].
Compared with the progress on the statistical analysis of VAR models, much less is known about their computational issues. For instance, Basu et al. [22] make exploration in this direction by proposing the fast network structure learning algorithm (FNSL) for the penalized recovery procedure. Their analysis does not take the structure information of the parameter into account and it only establishes a sub-linear convergence rate for the FNSL.
In this paper, we first provide the non-asymptotic optimization guarantee of PGD for VAR models with single-structured transition matrices. Our analysis illustrates that the distance between iteration points of PGD and the real transition matrix would converge linearly to the statistical error despite the objective function is not strongly convex under the high-dimensional settings. Our result is sharp in the sense that the minimal requirement of samples to guarantee the linear rate matches the phase transition of the model with independent samples up to a constant factor.
On the other hand, considering the parameter to be estimated only has one type of low-dimensional structure or is single-structured might be too oversimplified for messy real applications. So superposition-structured models have received more attention of researchers in last decade. Typical examples include robust PCA [23, 24], multi-task learning [25] and robust matrix sensing [26]. For these scenarios, we employ AltPGD to solve related optimization problems and establish the corresponding non-asymptotic optimization guarantee. Our results show that AltPGD also enjoys a linear convergence rate, which is much more efficient than the sub-linear rate in [22]. At the same time, our analysis avoids a drawback in [22] that the estimation error does not converge to zero when the number of measurements approaches infinity.
Last but not the least, we also illustrate that AltPGD is a practical algorithm to solve the general superposition-structured statistical model in [27]. Apart from the time series case, our analysis also adapts to multi-task learning and robust PCA.

heading: 2 Problem formulation
text: In this paper, we consider a d-dimensional vector-valued stationary time series {x0, · · · ,xn} generated by a VAR model of lag 1 with serially uncorrelated Gaussian errors. The VAR(1) model is defined as
xt+1 = Γ T ? xt + et+1, t = 0, · · · , n− 1, (1)
where Γ? ∈ Rd×d is the transition matrix and et iid∼ N (0,Σe). This model can be reformulated in the matrix form
Y = XΓ? +E, (2)
where Y = [x1, · · · ,xn]T ∈ Rn×d, X = [x0, · · · ,xn−1]T ∈ Rn×d, and E = [e1, · · · , en]T ∈ Rn×d. The goal of the VAR(1) model is to recover the transition matrix Γ? from the observation matrix Y and the data matrixX . In the high-dimensional settings with n d2, tractable recovery is possible when the transition matrix Γ? is well structured. Thus we introduce a convex regularizer R(·) to promote the structure of Γ?. Then a popular way to estimate Γ? is to solve the following constrained least square problem
min Γ
1
2n ||Y −XΓ||2F,
s.t. R(Γ) ≤ R(Γ?), (3)
where || · ||F represents the Frobenius norm of a matrix. When the transition matrix Γ? is superposition-structured in which Γ? is the sum of two singlestructured components, i.e., Γ? = S? + L?, we adopt two convex functions RS(·) and RL(·) to characterize the structures of its two components and solve the following constrained problem to
estimate S? and L?
min S,L
1
2n ||Y −X(S +L)||2F,
s.t. RS(S) ≤ RS(S?), RL(L) ≤ RL(L?).
(4)

heading: 3 Single-structured transition matrices estimation via PGD
text: In this part, we consider the case where the transition matrix to be estimated in (1) only has one type of low-dimensional structure which is characterized by a convex functionR(·). To estimate the transition matrix Γ?, we solve the problem (3) via the PGD update (summarized in Algorithm 1).
By setting fn(Γ) = ||Y −XΓ||2F/(2n), we could write the iteration of PGD as Γk+1 = PK(Γk − µ∇fn(Γk)), (5)
where µ is the step size, K = {Γ | R(Γ) ≤ R(Γ?)} is the descent set and PK represents the orthogonal projection onto the set K. We also introduce the concept of the descent cone C = cone(K − Γ?), where cone(·) is the conic hull of a set.
Algorithm 1 PGD for single-structured transition matrices estimation Input: Initial point Γ0, step size µ, iteration number K. for k = 0 to K − 1 do
Γk+1 = PK(Γk − µ∇fn(Γk)) end for Output: ΓK
To guarantee consistent estimation, we propose the following stability assumption for the VAR model (1), which is also imposed in [20, 22]. Assumption 1 (Stability). The characteristic polynomial of the VAR model (1) satisfies det(A(z)) 6= 0 on the unit circle of the complex plane {z ∈ C : |z| = 1}, where A(z) = Id×d − ΓT? z.
In the non-asymptotic analysis of the VAR model (1), we use the following two quantities M(fx) = ess sup
θ∈[−π,π] λmax(fx(θ)), (6)
m(fx) = ess inf θ∈[−π,π] λmin(fx(θ)), (7)
where fx(θ) is the spectral density function defined as
fx(θ) := 1
2π ∞∑ l=−∞ Σx(l)e −ilθ, θ ∈ [−π, π]. (8)
Here we useΣx(l) to represent
Σx(l) = E[xtxTt+l], t, l ∈ Z. (9) Specially, we writeΣx(0) = Σx for simplicity.
Compared with the model with independent samples, there is dependency among the rows of data matrix X in the VAR model (2), which is the main challenge when deriving the deviation bounds required by the estimation problem. For Gaussian processes, this dependency could be characterized by the covariance matrix Υx = E[vec(XT )vec(XT )T ], where vec(·) represents the column-wise vectorization of a matrix. The following lemma indicates that the concept of the spectral density function is a convenient tool to bound the extreme eigenvalues of Υx. Lemma 1 (Proposition 2.3 in [20]). Donate Υx = E[vec(XT )vec(XT )T ], where X is the data matrix in the VAR model (2). We could bound the extreme eigenvalues of Υx as
2πm(fx) ≤ λmin(Υx) ≤ λmax(Υx) ≤ 2πM(fx). (10) In particular, we also have
2πm(fx) ≤ λmin(Σx) ≤ λmax(Σx) ≤ 2πM(fx). (11)
For stable and invertible ARMA processes which include the model (1), the spectral density (8) has a closed form expression based on the matrix valued polynomials [20, Equation (2.4)]. Furthermore, the concrete calculation of the upper bound ofM(fx) and the lower bound of m(fx) for the model (1) is provided in [20, Proposition 2.2] which indicates m(fx) andM(fx) could be bounded away from zero and infinity. In this way, we could introduce the quantities κmin and κmax to simplify the expression of our analysis. Assumption 2 (Boundness). Suppose there are positive constants κmin and κmax satisfying
0 < κmin 2π ≤ m(fx) ≤M(fx) ≤ κmax 2π . (12)
With the above assumption, we could represent the extreme eigenvalues of Υx andΣx in a concise way
κmin ≤ λmin(Υx) ≤ λmax(Υx) ≤ κmax, (13) κmin ≤ λmin(Σx) ≤ λmax(Σx) ≤ κmax. (14)
In our analysis, we use the Gaussian width to quantify the size of a set T
ω(T ) := Esup x∈T 〈g,x〉, where g ∼ N (0, I).
We are now ready to present the non-asymptotic optimization guarantee of PGD for the problem (3). Theorem 1. Consider the VAR model (1) satisfying Assumptions 1 and 2. Suppose Γ? is singlestructured and R(·) is a convex function. Starting from a point Γ0 satisfying R(Γ0) ≤ R(Γ?), we solve the optimization problem (3) via PGD with the step size µ = 1/κmax. If the number of measurements satisfies
√ n > 2C κmax κmin (ω(C ∩ SF ) + u), (15)
then the PGD update (5) would obey
||Γk+1 − Γ?||F ≤ ρ k+1||Γ0 − Γ?||F +
ξ
1− ρ (16)
with probability at least 1− c exp(−u2). Here
ρ = 1− κmin κmax + C ω(C ∩ SF ) + u√ n < 1− κmin 2κmax , (17)
ξ = C ′ 1
√ κmax
||Σe|| 1 2 ω(C ∩ SF ) + u√
n , (18)
ξ
1− ρ =
1 1− ρ · C ′ 1√ κmax ||Σe|| 1 2 ω(C ∩ SF ) + u√ n < 2C ′ √ κmax κmin ||Σe|| 1 2 ω(C ∩ SF ) + u√ n , (19)
|| · || represents the spectral norm of a matrix, SF represents the sphere with unit Frobenius norm and c, C, C ′ are absolute constants. Remark 1 (Sharpness). Our result demonstrates that PGD can converge linearly to the statistical error despite the objective function fn(Γ) is not strongly convex under the high-dimensional settings. And the linear convergence is achieved when the number of measurements is of order ω(C ∩ SF )2, which matches the phase transition of the model with independent samples [15, 16]. Remark 2 (Impact of correlated samples). Our result also provide some insights for the impact of correlated samples. It is not hard to find that the temporal dependency is characterized by κmax and κmin which appear in the convergence rate, the estimation error, and the required number of samples. Let κ = κmax/κmin. Clearly, a smaller κ will lead to a faster convergence rate with smaller required samples and estimation error. Remark 3 (Comparison with related works). The existing works on structured VAR models such as [19–21] are mainly concerned with establishing the statistical error bounds for different recovery procedures. Our work focuses on algorithmic analysis and reveals how many samples would ensure that the algorithm achieves a fast convergence rate. On the other hand, our work could be regarded
as a generalization of the result in [17] from independent samples to time series. This generalization is nontrivial, since new mathematical tools (e.g., two deviation inequalities: Lemmas 4 and 5) are required to address time-series settings. To the best of our knowledge, Lemma 5 seems to be the first non-asymptotic result for the VAR models which yields the unified estimation error bounds for both independent and correlated samples. Furthermore, our analysis demonstrates PGD enjoys a linear convergence for structured VAR models, which is much more efficient than the sub-linear convergence rate of FNSL in [22]. Additionally, the linear convergence rate of PGD is also illustrated in [18] for the VAR(1) model with a sparse transition matrix, when the spectral norm of the transition matrix is strictly less than 1. Compared with the result in [18], our analysis adapts to general structured signals with the milder Assumption 1. Remark 4 (Extension). In Theorem 1, we set the step size µ = 1/κmax for a concise expression of the result. In fact, any step size satisfying µ ≤ 1/κmax could achieve a linear convergence rate by providing the corresponding number of measurements. In [20, 21], VAR(d) models are reformulated as VAR(1) models. With the same reformulation, our analysis also adapts to VAR(d) models.
Our analysis for the VAR model (1) also adapts to the multi-task learning problem with independent samples.1 Different from the time series setting, the measurements Y = XΓ? + E in this case are generated from xt iid∼ N (0,Σx) and et iid∼ N (0,Σe), for t = 1, · · · , n, where X = [x1, · · · ,xn]T ∈ Rn×d and E = [e1, · · · , en]T ∈ Rn×d are independent. Corollary 1. Consider the multi-task learning problem with the above conditions. Under Assumption 2, we solve the optimization problem (3) via PGD with the step size µ = 1/κmax and a starting point Γ0 satisfyingR(Γ0) ≤ R(Γ?). If the number of measurements satisfies
√ n > 2C κmax κmin (ω(C ∩ SF ) + u), (20)
then the update (5) would obey
||Γk+1 − Γ?||F < (1− κmin
2κmax )k+1||Γ0 − Γ?||F + 2C ′ √ κmax κmin ||Σe|| 1 2 ω(C ∩ SF ) + u√ n (21)
with probability at least 1− c exp(−u2), where c, C, C ′ are absolute constants.

heading: 4 Superposition-structured transition matrices estimation with AltPGD
text: In this part, we consider the case where the transition matrix Γ? to be estimated in (1) is superpositionstructured, that is, Γ? = S? +L?. To estimate S? and L?, we solve the optimization problem (4) via AltPGD (summarized in Algorithm 2).
We set fn(S,L) = ||Y −X(S +L)||2F/(2n) and write the update of AltPGD as Sk+1 = PKS (Sk − µ∇Sfn(Sk,Lk)), Lk+1 = PKL(Lk − µ∇Lfn(Sk,Lk)),
(22)
where KS = {S | RS(S) ≤ RS(S?)} and KL = {L | RL(L) ≤ RL(L?)}. We also introduce two descent cones CS = cone(KS − S?) and CL = cone(KL − L?), which would be used in our analysis.
Algorithm 2 AltPGD for superposition-structured transition matrices estimation Input: Initial points S0 and L0, step size µ, iteration number K. for k = 0 to K − 1 do Sk+1 = PKS (Sk − µ∇Sfn(Sk,Lk)) Lk+1 = PKL(Lk − µ∇Lfn(Sk,Lk))
end for Output: SK and LK
In this part, we considerRS(·) andRL(·) both belong to decomposable norms defined in [14]. 1In addition, our analysis framework could also be used in the problem to estimate Γ? and the precision matrixΣ−1e simultaneously [28].
Definition 1 (Decomposable norm). A regularization functionR(·) is decomposable with respect to a subspace pair (M,M⊥), if
R(α+ β) = R(α) +R(β), ∀α ∈M, β ∈M⊥. (23) Here,M is referred to as the model subspace which captures the constraints determined by the model andM ⊆ M. M⊥ is called the perturbation subspace indicating the deviation from the model subspaceM.
For common structure priors such as sparsity, group-sparsity and low-rank property, the corresponding regularization functions l1-norm, l1,2-norm and nuclear norm all belong to decomposable norms with low-dimensional model subspaces.
For the superposition-structured transition matrix Γ? = S? +L?, we assumeRS(·) is decomposable with respect to a subspace pair (MS ,M⊥S ), which is suit for the single-structured parameter S?. Similarly,RL(·) is decomposable with respect to a subspace pair (ML,M⊥L ) suit for L?. Due to the superposition-structured property of the transition matrix, we need to impose an additional assumption about the interaction between the two different structured components to guarantee the separate estimation. Assumption 3 (Structural incoherence). Given the subspace pairs (MS ,M⊥S ) and (ML,M⊥L ) for the two parameters S? and L?. Suppose the covariance matrixΣx defined in (9) satisfies
max { σ̄max(PMSΣxPML), σ̄max(PM⊥SΣxPML),
σ̄max(PMSΣxPM⊥L ), σ̄max(PM⊥SΣxPM⊥L ) } ≤ κmin 8 , (24)
where σ̄max(·) for a matrix Σ is defined as σ̄max(Σ) = sup V ,U∈SF 〈V ,ΣU〉 and κmin is defined in (12). Here PMS , PML , PM⊥S and PM⊥L donate the orthogonal projection operators onto the corresponding subspaces. Remark 5 (Related works). Several similar assumptions have been imposed in [27, 29, 30]. This type of assumptions is first proposed by Yang and Ravikumar in [27], where they use the structural incoherence assumption to restrict the interaction between different components of superpositionstructured statistical models and the C-Linear condition guarantees the structural incoherence under the linear regression setting and the Gaussian design. Meng et al. generalize the C-Linear assumption to the Structural Fisher Incoherence assumption in [29] for the estimation of sparse plus low-rank matrices in Gaussian Graphical Models. In [30], Greenewald and Hero introduce the structural incoherence assumption to robust Kronecker product PCA models. Our Assumption 3 is also motivated by [27] and would reduce to the C-Linear condition in [27] when we consider Σ-Gaussian ensemble where the rows ofX in (2) are generated independently from N (0,Σx). Remark 6 (Nonidentifiability). There are also other conditions used in literature to deal with the nonidentifiability concern. In [22], Basu et al. refer to the spikiness condition, which is first introduced in [31] for matrix completion and then is extended to matrix decomposition in [32]. In [27], Yang and Ravikumar compare the structure incoherence used here with the spikiness condition and illustrate the structure incoherence could address a drawback of the spikiness condition that the estimation error does not approach zero when the number of samples approaches infinity and requires weaker conditions at the same time. Another common condition used in [23–26] for sparse plus low-rank matrices recovery is the incoherence condition which is first proposed in [33, 34]. In [31, 32], the authors illustrate that the spikiness condition is a milder condition than the incoherence condition and is more suitable for the noisy models because of the consideration of singular values.
We now present the non-asymptotic optimization guarantee of AltPGD for the problem (4). Theorem 2. Consider the VAR model (1) satisfying Assumptions 1,2 and 3. Suppose Γ? is superposition-structured and Γ? = S? +L?, where S? and L? are two single-structured parameters whose structures are characterized by two decomposable normsRS(·) andRL(·) respectively. Starting from points S0 and L0 satisfying RS(S0) ≤ RS(S?) and RL(L0) ≤ RL(L?), we solve the optimization problem (4) via AltPGD with the step size µ = 1/κmax. If the number of measurements satisfies
√ n > 4C κmax κmin (ω(CS ∩ SF ) + ω(CL ∩ SF ) + u), (25)
then the update (22) would obey
||Sk+1 − S?||F + ||Lk+1 −L?||F ≤ ρ k+1(||S0 − S?||F + ||L0 −L?||F) +
ξ
1− ρ (26)
with probability at least 1− c exp(−u2). Here
ρ = 1− 3 4 κmin κmax + C ω(CL ∩ SF ) + ω(CS ∩ SF ) + u√ n < 1− κmin 2κmax , (27)
ξ
1− ρ =
1 1− ρ · C ′ 1√ κmax ||Σe|| 1 2 ω(CL ∩ SF ) + ω(CS ∩ SF ) + u√ n
< 2C ′ √ κmax κmin ||Σe|| 1 2 ω(CL ∩ SF ) + ω(CS ∩ SF ) + u√ n , (28)
and c, C, C ′ are absolute constants. Remark 7 (Related works). Our analysis makes progress on three aspects compared with the result in [22]. First, we illustrate the linear convergence rate of AltPGD compared with the sub-linear rate of FNSL. Second, our analysis indicates the requirement of samples for the linear rate and the statistical error, which are absent in the analysis of optimization in [22]. Third, Theorem 2 addresses a drawback of the result in [22] that the estimation error does not converge to zero when the number of samples approaches infinity.
Our analysis framework is also valid for robust PCA considered in [32, 27]. Suppose we have n i.i.d. sample zi ∈ Rd, where zi = ui + vi, ui ∼ N (0,L?), vi ∼ N (0,S?), ui and vi are independent. Here we set L? is a low-rank matrix and S? is a sparse matrix. We could write the sample matrix as Y = 1n ∑n i=1 ziz T i = L? + S? +E, where E = 1 n ∑n i=1 ziz T i − (L? + S?) is a Wishart noise matrix. In this setting, we have the data matrixX = Id×d and solve the problem
min S,L
1 2 ||Y − S −L||2F,
s.t. ||vec(ST )||1 ≤ ||vec(S T ? )||1, ||L||? ≤ ||L?||?,
(29)
where || · ||1 represents the l1-norm of a vector and || · ||? represents the nuclear norm of a matrix. Corollary 2. Consider the robust PCA model where S? is a sparse matrix with s? non-zero entries and L? is a r?-rank matrix. Under Assumption 3 where Σx = Id×d and κmin = κmax = 1 in this setting, we solve the optimization problem (29) via AltPGD with the step size µ = 1 and starting points S0 and L0 satisfying ||vec(ST0 )||1 ≤ ||vec(ST? )||1 and ||L0||? ≤ ||L?||?. If the number of measurements satisfies
√ n > C ′( √ s? log d+ √ r?d+ u), (30)
then the update of AltPGD would obey
||Sk+1 − S?||F + ||Lk+1 −L?||F
≤ (1 4 )k+1(||S0 − S?||F + ||L0 −L?||F) + 4 3 C||S? +L?||
√ s? log d+ √ r?d+ u√
n , (31)
with probability at least 1− c exp(−u2). Here c, C and C ′ are absolute constants.

heading: 5 Numerical results

heading: 5.1 Synthetic data
text: In this section, we apply PGD and AltPGD1 to network learning problems and compare the performance with FNSL proposed in [22]. We regularize sparse matrices by the l1-norm and low-rank matrices by the nuclear norm.2 All simulations are run on a PC with Intel i5-6500 and 16GB memory.
1The step size is selected as the inverse of the maximum eigenvalue of the estimated covariance matrix of the samples and the estimated covariance matrix is derived fromXTX/n.
2The projection onto the l1-norm ball follows the procedure in [35] and the projection onto the nuclear norm ball is a union of the singular value decomposition and the projection onto the l1-norm ball.
We first introduce several performance metrics for network estimation. For the estimated transition matrix Γ̂ and the real transition matrix Γ? whose entries are denoted by γ̂ij and γ?ij respectively, we define the true positive rate (TPR) and false alarm rate (FAR) as
TPR := ]{γ̂ij 6= 0 and γ?ij 6= 0}
]{γ?ij 6= 0} , FAR := ]{γ̂ij 6= 0 and γ?ij = 0} ]{γ?ij = 0} .
We also introduce the estimation error (EE), where EE := ||Γ̂− Γ?||F/||Γ?||F.

heading: 5.1.1 Network learning with a sparse transition matrix
text: First we consider Γ? ∈ Rd×d is a sparse matrix with s? non-zero entries. We suppose each row of Γ? has s?/d non-zero entries whose values follow a standard normal distribution. Then we rescale Γ? to guarantee the stability of the process. In this simulation, we set d = 100 and s? = 3500. To illustrate the effect of the numbers of samples, we perform the simulation under three scenarios n = 1000, 1500, 2000 and each scenario is repeated for 100 trials. For FNSL, we choose the regularization parameter λS as O( √ n log d) according to Propositions 1 and 3 in [22]. Both algorithms start from Γ0 = 0.
In Table 1, we record the experimental results for the two algorithms under different numbers of samples. The results illustrate that PGD enjoys better performance with much less computation time than FNSL and support our analysis in Theorem 1.

heading: 5.1.2 Network learning with a low-rank transition matrix
text: Then we consider Γ? ∈ Rd×d is a low-rank matrix whose rank is r?. Suppose Γ? is constructed by Γ? = UV
T , where U ,V ∈ Rd×r? are matrices with independent standard Gaussian entries. We also rescale Γ? to guarantee the stability of the process. We set d = 100, n = 8000, r? = 2 and repeat the scenario for 100 times. For FNSL, we choose the regularization parameter λL as O( √ nd) according to Propositions 1 and 3 in [22]. The result in Figure 1(a) illustrates that the updates of PGD enjoy a faster convergence rate than those of FNSL, as predicted in Theorem 1.
5 10 15 20 25 30 35 40 45 50
10 0
10 1 FNSL
PGD
1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
d = 50 d = 75 d = 100
0.5 1 1.5 2 2.5 3
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
d = 50 d = 75 d = 100
We also perform simulations under different dimensions d and different numbers of samples n to verify the order of estimation error, where we set r? = 4 and each scenario is repeated for 100 times. The results in Figure 1(b) and 1(c) indicate that all the squared empirical error curves behave as f(t) ∝ t−1 and support our theoretical results in Theorem 1.

heading: 5.1.3 Network learning with a superposition-structured transition matrix
text: In this part, we suppose the transition matrix in the VAR model (1) is superposition-structured. Specially, we consider Γ? = S? + L?, where S? is a sparse matrix with s? non-zero entries and L? is a rank-r? matrix. The construction of S? and L? follows the same procedure as the above simulations. First, we set d = 100, s? = 3500 and r? = 2. To illustrate the effect of the numbers of samples, we perform the simulation under three scenarios n = 1500, 2000, 2500 and each scenario is repeated for 100 trials. For FNSL, we choose the regularization parameters λS as O( √ n log d) and
λL as O( √ nd) according to Propositions 1 and 3 in [22]. Both algorithms start from S0 = 0 and L0 = 0.
In Table 2, we record the experimental results for the two algorithms under different numbers of samples. The results illustrate that AltPGD enjoys better performance with much less computation time than FNSL.
Then we compare the convergence rates of AltPGD and FNSL. We set d = 100, n = 8000, s? = 3500, r? = 3 and repeat the scenario for 100 times. The result in Figure 2(a) illustrates the efficiency of AltPGD.
We also perform simulations under different dimensions d and different numbers of samples n to verify the order of estimation error, where we set s? = 300, r? = 3 and each scenario is repeated for 100 times. The results in Figure 2(b) and 2(c) indicate that all the squared empirical error curves behave as f(t) ∝ t−1 and support our theoretical results in Theorem 2.
0 50 100 150
10 0
10 1
FNSL AltPGD
1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
0.5
1
1.5 d = 50
d = 75 d = 100
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.5
1
1.5 d = 50
d = 75 d = 100

heading: 5.2 Real data
text: Next, we analyze the temporal dynamics of the log-returns of stocks in the S&P 500 index. The stock data consists of 1259 daily closing prices for 434 companies in the S&P 500 index between February 8, 2013 and February 7, 2018 [36]. In this way, we get 1259 data vectors, each of which contains the closing prices of all stocks on a trading day. To ensure the data stationary, we calculate the log-returns {rt}T−1t=1 of stocks by
rt,i = log( pt+1,i pt,i ), t = 1, · · · , T − 1, (32)
where pt,i is the closing price of stock i at day t. In this way, we construct the data matrix X = [r1, · · · , rT−2]T and the observation matrix Y = [r2, · · · , rT−1]T . We adopt the VAR model (1) with the regularizer R(·) = || · ||1 to study the evolution of stock log-returns over the 2013-2018 period and then solve the model with PGD. The constrained parameter is selected through 5-fold cross validation. In Figure 3, we present the sparsity patterns of two parts of the transition matrix Γ̂ estimated by PGD, which indicate meaningful Granger causal effects [37, 38] among the log-returns of stocks. The 434 companies belong to 10 different sectors, such as materials (22 stocks), energy (29 stocks), consumer staples (31 stocks) and financials (83 stocks). In Figure 3(a) and 3(b), the log-returns of stocks in the energy sector have stronger Granger causal effects, and the Granger causal effects of the consumer staples sector and the financials sector are weaker. For comparison, we present the sparsity patterns of the transition matrix estimated by FNSL in [22], which illustrate a similar Granger causal network.
Acknowledgments and Disclosure of Funding
The authors would like to thank the anonymous reviewers for their constructive feedback and suggestions. This work was supported by the National Natural Science Foundation of China under Grants 61971044 and 62025103.

Review summaries:
This paper presents a recovery guarantee for gradient-based optimization method applied to high-dimensional vector autoregressive (VAR) models. It considers both single-structured and superposition-structured transition matrix, and shows that the distance between the true parameter and the iterates generated by projected gradient descent (PGD) algorithm will decrease linearly to certain statistical error controlled by sample size and Gaussian width of the assumed structure. In the experiments, PGD is shown to outperform the baseline algorithm FNSL empirically.

The paper proposes projected gradient descent algorithms for estimating high-dimensional vector auto-regressive models, and conducts convergence analysis following the restricted strong convexity lens.  Specifically, sparse, low-rank and sparse plus low-rank models are of interest.  While many of these topics have been studied "to death", the novelty here seems to be (a) there's limited work in applying sparse plus low-rank approaches to time-series settings, where there is an additional constraint that the transition matrix has to lead to a stable dynamical system (b) Typical convergence analysis is done under iid setting, but for VAR models, samples are naturally dependent with strong auto-correlation.

This paper considers the problem of estimating the transition matrix of a lag-1 vector autoregressive model. The authors propose projected gradient descent as well as alternating projected gradient descent algorithms to learn both single- and supervision-structured transition matrices, respectively. Non-asymptotic analysis is also provided which shows the algorithms convergence linearly tp the statistical error in the high-dimensional settings.

This manuscript proposes two algorithms for estimation of transition matrices in single-structured as well as superposition-structured VAR models. Linear convergence to true parameters are verified under certain stability conditions, and additional assumptions on the data distribution, step size and minimum sample size. Some empirical results are reported to illustrate good performance of proposed algorithms in finite sample cases.