Missing Data in Sparse Transition Matrix Estimation for Sub-Gaussian Vector Autoregressive Processes
Amin Jalali Rebecca Willett

Identifying Young Kepler Planet Host Stars from Keck-HIRES Spectra of Lithium1
    Travis A. Berger2, Andrew W. Howard2,3, & Ann Merchant Boesgaard2
    
==============================================================================

empty




High-dimensional time series data exist in numerous areas such as finance, genomics, healthcare, and neuroscience. An unavoidable aspect of all such datasets is missing data, and dealing with this issue has been an important focus in statistics, control, and machine learning. In this work, we consider a high-dimensional estimation problem where a dynamical system, governed by a stable vector autoregressive model, is randomly and only partially observed at each time point. Our task amounts to estimating the transition matrix, which is assumed to be sparse. In such a scenario, where covariates are highly interdependent and partially missing, new theoretical challenges arise. While transition matrix estimation in vector autoregressive models has been studied previously, the missing data scenario requires separate efforts. Moreover, while transition matrix estimation can be studied from a high-dimensional sparse linear regression perspective, the covariates are highly dependent and existing results on regularized estimation with missing data from i.i.d. covariates are not applicable. At the heart of our analysis lies 1) a novel concentration result when the innovation noise satisfies the convex concentration property, as well as 2) a new quantity for characterizing the interactions of the time-varying observation process with the underlying dynamical system. 




§ INTRODUCTION


Consider a p-dimensional covariance-stationary vector autoregressive model of lag one, namely a (1), as

    w_t+1 =  w_t + ϵ_t  ,   t=0,…, -1,
	
where ∈ℝ^p× p is the corresponding transition matrix, and each ϵ_t is a p-dimensional vector of innovations, with zero mean and covariance Σ_ϵ, that is temporally uncorrelated with other noise vectors. The goal is to estimate  from partial observations of entries of w_0,…, w_, given prior knowledge on  being sparse. 

Concatenating all the vectors in (<ref>) as columns of matrices yields

    [ w_1   ⋯   w_ ]_𝒴
    	= [ w_0  ⋯  w_-1 ]_
    	+ [ ϵ_0  ⋯ ϵ_-1 ]_ℰ

where each of the brackets represent a (p×)-dimensional matrix. 
The available information from (<ref>) is in the form of entries in W = [ w_0  ⋯  w_ ] that are missing according to i.i.d. Bernoulli random variables with probability 0 ≤δ <1. Consider a new process {_t} where, for any i=1,…,p,

    (_t)_i = 
     	(w_t)_i    with probability  1-δ
    
     	0    with probability δ  ,

and observation is independent for different i=1,…,p and different t=1,…,. 
We use the bar notation for other objects constructed from _0,…,_. For example, W̅ = [ _0  ⋯ _ ], = [ _0  ⋯ _-1 ], and so forth. For simplicity, we consider a centered process, hence w_0=0. 


We handle the missing data by modifying the LASSO <cit.> using population information on the observation pattern, as described in the Appendix. 
More specifically, similar to <cit.>, we solve either of the two following constrained versions of this program: either 

    _B_1 ≤ b_0 √(k)1/B  - _F^2 
    	- δ/B_F^2
    	+ (1-δ)^2 λ_B_1

where 
= (diag( '))^1/2∈ℝ^p× p is a diagonal matrix of sample autocovariances for each of the p covariates, 
k is the number of nonzero entries of , i.e., 
k = _0, 
b_0 is any value at least equal to _F, 
and λ_ is the regularization parameter that will be chosen according to the parameters of the problem, 
or

    _B_1 ≤_1  1/B  - _F^2 
    	- δ/B_F^2  .

Note that, with a possibly non-convex quadratic optimization program, we need a suitably constrained feasible set to avoid an unbounded optimization problem and hope for recovering the target model, hence the constraints in (<ref>) and (<ref>).  


Through employing a well-known machinery for the analysis of LASSO (summarized as Theorem <ref> in the Appendix), and by developing new concentration results for the random processes of interest in this work (Proposition <ref>), we provide guarantees on the ℓ_1- and ℓ_2-norm estimation errors for (<ref>) and (<ref>) in Theorem <ref>. Before stating the main result, we discuss the prior art and how our setup leads to new challenges. We then discuss certain characteristics of the processes in (<ref>) and (<ref>) that are used in our guarantees, in Section <ref>. More specifically, we introduce a new quantity, namely () in (<ref>), which is used in characterizing the interplay between the dependence among covariates and the difficulty of recovery from partial information. We elaborate on these characteristics in Section <ref>. Section <ref> provides a sketch of the proof for providing error bounds on LASSO and its variants, and is similar to many other works on LASSO in the literature. Section <ref> contains our main contribution on establishing the required concentration inequalities for providing error bounds for LASSO through concentration of sub-Gaussian quadratic forms. 




 §.§ Prior Art and New Challenges

Estimators (<ref>) and (<ref>) can be seen as modifications of the LASSO <cit.>, with constraints that help remedy the possible non-convexity of the estimator. Similar estimators have been considered in the literature for several sparse regression tasks and a similar framework has been used to analyze such estimators; e.g., in <cit.>, and many more. However, the distinguishing aspect of different works in this area is the difference in the data generation processes and the required concentration analysis. In a simple data generation scheme as y = X + ϵ where ϵ∼𝒩(0,I) and X has i.i.d. random entries drawn from 𝒩(0,1) independently from ϵ, the analysis of LASSO, 

    _β 1/Xβ -y _2^2 
    	+ λβ_1

boils down to understanding the spectrum of the random matrix X. A more complicated case of correlated Gaussian designs is considered in <cit.>. More involved data generation scenarios require more involved probabilistic arguments to establish the conditions that guarantee (near-)optimality of  for (<ref>). For example, <cit.> extends the above to transition matrix estimation in Gaussian vector autoregressive models. Authors in <cit.> extend (<ref>) to the case of non-convex quadratic optimization programs when X has one of several interesting dependency patterns. 

We note that our focus is different from <cit.> which considers sparse regression with missing data when the design matrix is generated by an autoregressive process with known transition matrix A satisfying A_2<1. In our case, the interactions among covariates depend on the unknown transition matrix, making the problem even more challenging. 

Assuming Gaussian innovations, ϵ_0,…,ϵ_-1, makes {w_t} of (<ref>) a Gaussian process. However, the partially observed process {_t} will no longer be a Gaussian process. Nonetheless, {_t} belongs to the family of sub-Gaussian processes for which many properties are known. In this work, we consider a subset of sub-Gaussian processes for the innovations: those with the convex concentration property in Definition <ref>. 





 §.§ Matrix Quantities



We now review some important quantities associated to the transition matrix of interest, . Consider  as the adjacency matrix of a weighted directed graph on p nodes as in the left panel of Figure <ref>. It is intuitive that not only the number of edges in such graph but also the configuration of edges, e.g., the degree distribution, plays an important role in any inverse problem for identifying such graph. In the following, we present important quantities associated with  that help in reflecting nuances of the corresponding graph. From a dynamical systems point of view, we review the notions of the spectral radius and the spectral norm, as well as new quantities presented in Section <ref>, which allow us to capture the stability of the autoregressive process as well as the information content of our observations from this system.  



   
   
   
Denote by ℂ^p× p the space of all p by p complex matrices. The spectral radius of B∈ℂ^p× p is the non-negative real number 

    ρ(B) =max{λ: λ∈σ(B)}

where |·| denotes the magnitude and σ(B) is the set of all eigenvalues of B. A norm on ℂ^p× p is called a matrix norm if it satisfies the submultiplicative property as AB≤AB for all A,B∈ℂ^p× p. If · is any matrix norm, then ρ(B) ≤B. 
For any ϱ,ι≥ 1, consider the corresponding ℓ_ϱ and ℓ_ι vector norms. The corresponding induced operator norm is then defined as 

    B_ι→ϱsup_x≠ 0Bx_ϱ/x_ι.

For example, B_1→ 2 is the maximum ℓ_2 norm among the columns of B, and B_2→∞ is the maximum ℓ_2 norm among the rows of B. When ι=ϱ, we simply denote B_ιB_ι→ι. For example, B_1 = max_1≤ j ≤ p∑_i=1^pB_ij and B_∞ = max_1≤ i ≤ p∑_j=1^pB_ij. When ι=2, B_2 = √(Λ_max(B'B)) is the spectral norm which is also simply referred to as the operator norm. Notice that B_2 ≠ρ(B) in general. The extension of vector ℓ_ι norms to matrices is denoted by B_ιvec(B)_ι. Finally, we abuse notation to denote by B_0 the number of nonzero elements in the vector or matrix input. The interested reader is referred to Section 5.6 of <cit.> for a more comprehensive account of these matrix quantities.


Going back to the directed graph interpretation of B, we can view B_1→ 2 as the maximum energy that any one node can exert on other nodes which is related to the maximum out-degree of nodes, while B_2→∞ can be viewed as the maximum energy that is being exerted upon each node and is related to the maximum in-degree among nodes. 




 §.§ Salient Characteristics of the Model


The autoregressive process in (<ref>) is called stable if and only if ρ()<1: all eigenvalues of  have modulus less than one. This is equivalent to

    (I- z)≠ 0   for all   |z|≤1  ,

where z is a complex-valued scalar variable. In such case, we define three main quantities 

    (B)    max_z=1 I-Bz_2
    
    (B)    max_z=1 (I - Bz)^-1_2
    
     (B)    max_z=1 (I - Bz)^-1_1→ 2 .

The first two quantities are related to the least and the largest singular values of the transfer function on the unit circle. In other words, for all z with z=1, they quantify the least and the largest values of (I-Bz)^-1u_2 when u_2 = 1. The third quantity, on the other hand, characterizes the largest value of (I-Bz)^-1u_2 when u_1 = 1, for all z with z=1. 

Moreover, for k = B_0, we have (see Section <ref>)

    (B) ≤ (B) ≤√(2k)  (B) .

While (<ref>) and (<ref>) have been considered in other works on autoregressive models (e.g., see μ_max and μ_min in Equation (2.6) of <cit.>), the definition of (<ref>) in the context of autoregressive model estimation is, to the best of our knowledge, new and motivated by the missing data setup. 




§ MAIN RESULTS

In this section, we state our main result followed by a discussion on the main quantities, a sketch of the proof, and a list of ingredients for this proof that we establish in the subsequent sections. In essence, we would show that the error scales with

    θ_0 ()^2/()^2∈ [1/2k, 1] .

We define a few more quantities to simplify the presentation of our main result. First, consider an innovation condition number defined as 

    κ_ϵ    36 √(c_ac_ϵ^2Σ_ϵ_2)Σ_ϵ^-1_2

where c_a is a global constant and c_ϵ is a function of the innovation process ϵ and will be defined later. We also consider a a transition condition number defined as

    κ_0    ()^2 ()^2 = 
    max_z=1σ_max^2(I- z)/min_z=1σ_min^2(I- z) .

Finally, the nonzero pattern of  and the quality of our choice for b_0 can be measured through the followings: 

    h b_0/7 (_2→∞^2 + 1) √(k)  ,  ζ1 + δθ_0  k/hk .

Note that 
_F ≤√(k)_2→∞≤√(k)_2→∞^2/ρ(), which helps in understanding h through 

    _F/(_2→∞^2 + 1) √(k)   ≤min{1/ρ() , _2→∞} ,

as b_0 is chosen to be at least _F.



Consider the p-dimensional autoregressive process in (<ref>) satisfying ρ()<1. Suppose _0 = k and the innovations are temporally uncorrelated with zero mean and a positive definite covariance matrix Σ_ϵ, and satisfy the convex concentration property (Definition <ref>) with constant c_ϵ. Suppose we have partially observed the process, with missing probability δ, for time length  satisfying 

    √(/log p)≥ κ_ϵκ_0 ζ/(1-δ)^21/(1/27ζ  - δθ_0  )^2

while ζ > 27 δθ_0. 
Define 

    Φc_0 b_0√(k)/7κ_ϵκ_0 ζ/(1-δ)^2√(log p/) 	.

Consider any b_0 ≥_F and any λ_ satisfying λ_≥2Φ/φ_0 where 

    φ_0     c_0()^2
    	Σ_ϵ^-1_2

and c_0 and c_1 are universal constants. Then, with probability at least 1-10 p^-1, for any optimal  in (<ref>) we have 

    - _F  ≤  2 √(k)φ_0
    	λ_ ,    - _1  ≤  16 k φ_0
    	λ_

and, for B̃[_ij1_ |_ij|>λ_]_i,j=1,…,p we have

    supp(B̃) ∖supp()   ≤   112 k φ_0 .

The same bounds, where λ_ is replaced by 2Φ/φ_0, apply to (<ref>).




As mentioned before, there is a well-developed theory for providing guarantees on the performance of the LASSO and its variants which we summarize as Theorem <ref> in the Appendix. This framework requires establishing certain concentration properties for the underlying data generation process. We provide the required concentration results in Section <ref> and combine them with the framework of Theorem <ref> to prove Theorem <ref>. In the following, we first provide further intuition into the results of Theorem <ref> in Section <ref>. We then elaborate on the required conditions for proving Theorem <ref> and motivate our concentration result in Section <ref>. 




§ REMARKS ON THE MAIN QUANTITIES

The quantities appearing in Theorem <ref> worth further discussion. In this section, we provide further details on different quantities we defined in relation to the transition matrix of interest, . 




 §.§ The Support

First, as we will see in the concentration result, () appears because of the missing data setup; specifically, due to the term - B_F^2 in (<ref>) and (<ref>). Intuitively, we expect that the support of , and how each covariate affects the value of other covariates in the next time step (see Figure <ref>), should play an important role in our ability in recovery from missing data. For example, if  is diagonal, then covariates are temporally uncorrelated (do not directly affect each other over time) and the entries of  have to be estimated independently. On the other hand, for more distributed supports of , we experience two competing phenomena:


  * when each covariate is influenced by many covariates from the previous time point, residuals between  and  are generally smaller because the value of missing covariates play less of a role, so recovery is more robust to missing data. This is captured by θ_0 (defined in (<ref>)) in our results. 

  * higher dependence among covariates makes observations more highly correlated and the resulting inverse problem becomes more ill-posed even when we have complete data. This is captured by κ_0 (defined in (<ref>)) in our results. 



In short, not all k-sparse  are equally easy or difficult to infer from incomplete data. For example, if only one column of  is nonzero (in-star graph in Figure <ref>), then one element of w_t is influenced by the previous realizations of the process, while the other covariates are not. If only one row of  is nonzero (out-star graph in Figure <ref>), then all covariates are being influenced by the same single covariate and there are no other influences. Finally, if  is nonzero on a single off-diagonal, then the i-th covariate is only influencing covariate i+1, for i=1,…,p-1, corresponding to a chain graph representation of influence structure. 








   





 §.§ Dimension-independence

Denote the set of nonzero rows of B by J_r⊆{1,…,p} and the set of its nonzero columns by J_c⊆{1,…,p}. For k = B_0, it is easy to see that J_r≤ k and J_c≤ k, with · denoting the size of the set. Moreover, for J J_r∪ J_c, all of the nonzero entries of B are in a principal submatrix indexed by J. Therefore, for any integer value t≥ 1, all of the nonzero entries of B^t are in the same principal sub-matrix indexed by J. Considering the Neumann series (I-A)^-1=∑_t=0^∞ A^t when ρ(A)<1, the above implies that (B), (B), and (B), are only concerned with the smallest principal submatrix of B containing all of its nonzero entries and are independent of the dimension of B: embedding B into a larger zero matrix does not change these values, as desired. 

It is worth mentioning that the same conclusion, of independence from the ambient dimension, cannot be made about the quantity ℳ(f_w,) used in <cit.> as the innovations could make the time series fully supported over all entries. 



 §.§ Bounds

While B and the quantities in (<ref>), (<ref>), and (<ref>), do not directly scale with each other, they have a close relationship that can be used in better understanding the main theorem. 

	Suppose (I-Bz)≠ 0 for all |z|≤1. Then, 
	
    (B) 
    		≤ 1 + B_2
    		≤ 1+B_1 +  B_∞/2.

Moreover, if B is diagonalizable, then

    (B) 
    	≤1/1-ρ(B)R_2R^-1_2

where the columns of R are the eigenvectors of B. 

Moreover, submultiplicativity of induced operator norms provides

    (B) 
    ≥max_z=1 (I-Bz)_1→ 2^-1≥ (1+ B_1→ 2)^-1 .

Therefore, when B is diagonalizable, 

    √(θ_0)=	()/()≥1-ρ(B)/ (1+ B_1→ 2) R_2R^-1_2

and 

    √(κ_0) 
    =	(B) (B) 
     ≤R_2R^-1_2/1-ρ(B)  (1+B_1 +  B_∞/2).






 §.§ A Restrictive Assumption We Avoid
 
In this paper, we only assume stability, i.e., ρ() < 1. This assumption is milder than the more stringent condition _2 < 1 prevalent in the literature. Only requiring the milder assumption used in this paper has important practical consequences. While _2<1 implies ρ()<1 (hence the stability of the corresponding autoregressive process), _2<1 is necessary only when  is symmetric. In other words, an asymmetric matrix  with _2 ≥ 1 can correspond to a stable autoregressive process; e.g., see Lemma E.1 in <cit.>. For example, for some 0<a<1, the matrix

    = [   a 1/a;   0   a ]

has eigenvalues equal to a, hence a spectral radius of a<1, but an operator norm that is slightly larger than 1/a > 1. Not assuming a spectral norm bound on the transition matrix becomes important in the study of vector autoregressive processes with a lag larger than one, defined as 
	
    w_t = B_1 w_t-1 + B_2 w_t-2 + … + B_d w_t-d + ϵ_t ,

	where d≥1 is the lag. It is easy to see that the above can be reformulated as a vector autoregressive process with lag one, as
	
    [     w_t;   w_t-1;       ⋮; w_t-d+1 ]= [   B_1   B_2     ⋯ B_d-1   B_d;   I_p     0     ⋯     0     0;     ⋮     ⋮     ⋱     ⋮     ⋮;     0     0     ⋯   I_p     0 ]_B[ w_t-1; w_t-2;     ⋮; w_t-d ]+[ ϵ_t;   0;   ⋮;   0 ].

Lemma E.2 in <cit.> establishes the fact that d>1 implies B_2≥ 1, even when ρ(B) < 1, illustrating the restrictiveness of operator norm bound assumptions. 





§ ESTIMATION ERROR FOR NON-CONVEX LASSO
 
Both (<ref>) or (<ref>) can be viewed as constrained quadratic optimization programs, 

    ∈_ B∈ℬ tr(B  B') - 2⟨ B ,  ⟩ + λB_1 ,

for corresponding choices of the constraint set ℬ⊂ℝ^p× p and regularization parameter λ, where 

    = 1/ ( '  - δdiag( ') )  
     ,    =1/'	.

In this work, we are not concerned with the possible non-convexity of these estimators from a computational point of view and focus on the statistical performance. Nonetheless, simple algorithms such as variants of projected gradient descent can be used for convergence to a small neighborhood of the set of all global minimizers, similar to <cit.>. We postpone such convergence guarantees to future work. 

LASSO <cit.> and its variants have been used and studied extensively in the literature. We specifically use a popular approach for providing guarantees on the estimation performance of LASSO and its variants presented in <cit.>. We present a version of this result, tailored to norm-constrained ℓ_1-regularized non-convex quadratic optimization, as Theorem <ref> in the Appendix. 

Theorem <ref> is essentially the same in any work on guarantees for LASSO and its variants, but lumps all the mechanical, and now well-known, parts of the process in one theorem and is discussed for clarity of our exposition. For example, the theorem can be seen as an extension of <cit.> and <cit.> for non-convex LASSO and an extension of <cit.> for transition matrix estimation in multivariate time series. 


To provide ℓ_1 and ℓ_2 norm error bounds for such estimates, following the framework developed in <cit.>, we need to establish the so-called lower restricted eigenvalue condition, stated equivalently <cit.> as 

    v'  v ≥α_v_2^2 - τ_v_1^2 
       for all  v∈ℝ^p	,

as well as a deviation bound, 

    - _∞≤ 	c √(log p/) ,

where c depends on the parameters of the problem instance. 

Since Q and L come from samples generated by the partial observation of a vector autoregressive model, they are random objects and reasonable values of α_, τ_, and c, in (<ref>) and (<ref>), may be used only with high probability. Therefore, we use relevant concentration results to establish these bounds with high probability. 

In the following, we expand the conditions in (<ref>) and (<ref>) and represent them as simple functions of the autoregressive process, which will then be bounded in Section <ref> using results on concentration of sub-Gaussian quadratic forms. Let us fix some notation first. For a p-dimensional discrete-time, centered, covariance-stationary (wide-sense stationary) process {w_t}, denote the autocovariance function by 

    Γ_w(h) = (w_t, w_t+h) .

For a matrix A, the transpose is denoted by A'. Denote by ⊙ and ⊘ the Hadamard (element-wise) product and division respectively, and by ⊗ the Kronecker product. The covariance matrix for the Bernoulli mask characterized in (<ref>) is given by 

    P = (1-δ)^2 1 + δ(1-δ) I,

so that = 1/' ⊘ P. 




 §.§ Restricted Eigenvalue Condition

Observe that = Γ_w(0) = Γ_(0) ⊘, which gives 

    - 
    = (1/' - Γ_(0)) 
    - δ (1/' - Γ_(0)) ⊙ I .

Then, bounding v' ( - )v, for all v∈ℝ^p, allows for establishing (<ref>) through the application of the triangle inequality. Suppose we established the following condition for a fixed value of  which will be determined later: 




  * 
For any fixed v∈ℝ^p with v_0 ≤ 2 and v_2 = 1, there exists η(s) such that v' ( - )v ≤η() with probability at least 1-p_1(). 


Then, such concentration can be stated over the set of 2-sparse vectors using a discretization argument, as in Lemma F.2 of <cit.>, followed by a simple argument that relates the set of sparse vectors to those with a bounded ℓ_1 norm, as in Lemma 12 of <cit.>. As the above calculations depend on the free parameter , it should be chosen in a way that makes η() as small as possible while maintaining the probability for (<ref>), which depends on p_1() and the union bound in the discretization step, at a desired level. We specify our choice of  for the proof of Theorem <ref> right after the statement of Theorem <ref> in the Appendix. 




 §.§ Deviation Bound

The matrix of interest in (<ref>) is given by

    - 
    = (1/' ⊘ - Γ_w(0)) 
     - 1/(1-δ)^2 (1/' - Γ_(1))	'

where we used the fact that Γ_w(0) = Γ_w(1)' and Γ_w(1) = 1/(1-δ)^2Γ_(1). 
The first assertion considers full information and is related to the interaction of {w_t} and {ϵ_t} processes. In fact, using the original process in (<ref>) we get

    Γ_w(1) -Γ_w(0) '
       = (w_t, w_t+1)  - (w_t, w_t) '	
       = (w_t, w_t+1)  - (w_t,  w_t) 	
       = (w_t, ϵ_t)
 
which is zero in our setup. 
For clarity, we state (<ref>) as another condition:



  *  
There exists φ>0 such that - _∞≤φ with probability at least 1-p_2.


Therefore, to derive the desired bounds in (<ref>) and (<ref>) and provide ℓ_1 and ℓ_2 norm error bounds for (<ref>) and (<ref>), we can establish <ref> and <ref>; a complete description of this procedure is stated as Theorem <ref> in the Appendix. This amounts to computing concentration bounds on the four terms in (<ref>) and (<ref>). To that end, we rewrite our process in matrix form and leverage classical linear time invariant dynamical systems to establish several quantities that characterize the process in (<ref>) and will appear in those concentration bounds. These relationships are summarized in Lemma <ref> and lead to the main concentration result given in Proposition <ref>. 




§ CONCENTRATION OF SUB-GAUSSIAN QUADRATIC FORMS

As mentioned before, establishing either of the conditions (<ref>) and (<ref>) relies on certain concentration properties for the underlying data generation process that defines  and . In the following, we make this relationship concrete and provide the main concentration result in Proposition <ref>. 

To establish <ref> for (<ref>) (to get (<ref>)), we are interested in the concentration of v'  ' v and v'  (' ⊙ I) v around their mean, for any fixed v. In the following, we manipulate these quantities into convex quadratic forms in terms of the noise vector 

    _' [ w_0'  ϵ_0'  ϵ_1'  ⋯ ϵ_-2' ] .

Define I_Ω∈{0,1}^p× p as the diagonal matrix whose (pt+ j)-th diagonal entry is one if (w_t)_j is observed and zero otherwise, for t=0,…,-1 and j=1,…,p. Moreover, define

    Ψ_(B) = [    I    0    0    ⋯    0;    B    I    0    ⋯    0;    ⋮    ⋮    ⋮    ⋱    ⋮; B^-1 B^-2 B^-3    ⋯    I ]

which is a block-Toeplitz matrix. Then, v' ' v  = _' Ψ_(1) ' Ψ_(1)_ and v' (' ⊙ I) v = _' Ψ_(2) ' Ψ_(2)_ where 

    Ψ_(1)    (I_n ⊗ v)' I_ΩΨ_(B) 
    Ψ_(2)    (I_n ⊗diag(v)) I_ΩΨ_(B) .

The latter is because of the following, 

    v' ( ' ⊙ I) v 
       = ⟨' ⊙ I, vv' ⟩ 
    = ⟨' , vv' ⊙ I ⟩
       = ⟨' , diag(v)^2 ⟩ 
    = diag(v) _F^2 
       = vec(diag(v) )_2^2
       = (I_⊗diag(v) )vec()_2^2 
        = (I_⊗diag(v) ) I_ΩΨ_(B)  __2^2  .

The concentration of the above two quadratic forms, in Ψ_(1) and Ψ_(2), can be studied when we assume the so-called convex concentration property on noise vectors ϵ_t, for t=0,…, -1, or equivalently on the noise vector _. 
 
Let x be a random vector in ℝ^n. We will say that x has the convex concentration property with constant c_x if for every  convex function g:ℝ^n→ℝ, we have g(x)<∞ and for every t>0, 

    ℙ[   g(x) -  g(x) ≥ t  ] ≤ 2 exp(-t^2/c_x^2)  .


If the above tail bound holds for all functions g(x) = ⟨ x, u ⟩ where u∈ℝ^ is any vector with u_2=1, then x is called a sub-Gaussian random vector <cit.>. However, the convex concentration property requires such tail bound to hold for every 1-Lipschitz convex function, and characterizes a subclass for the sub-Gaussian random vectors. See <cit.> for examples of such random vectors. As pointed out by <cit.>, 2c_x^2 ≥Σ_x_2 always holds. 

Improving upon a bound in <cit.>, 
Theorem 2.5 in <cit.> allows for bounding the deviations of our quadratic forms from their mean, as a function of 
Ψ_(i)'Ψ_(i)_2 =Ψ_(i)_2^2 and Ψ_(i)'Ψ_(i)_F^2, which is at most v_0 Ψ_(i)_2^4, for i=1,2. 

These operator norms can be related to certain norms of the block-Toeplitz matrix Ψ_(B). This matrix can in turn be related to a a transfer function that is used in the definitions of (B) and (B). The result is summarized in the next lemma whose proof is given in the Appendix. 



With the above notation, the followings hold

Ψ_(1)_2 
     ≤  Ψ_(B)_2
        ≤  (B)

 
Ψ_(2)_2 
     ≤  Ψ_(B)_1→2
        ≤  (B).


All in all, we get the following concentration result. 

Consider the autoregressive time series in (<ref>) where all ϵ_t, for t=0,1,…, -1, are temporally uncorrelated, have zero mean and variance Σ_ϵ, and satisfy the convex concentration property with constant c_ϵ. Moreover, consider {_t} as the partially observed time series corresponding to {w_t}, as characterized in (<ref>). Then, for any fixed vector v∈ℝ^p with v_0 ≥ 2_0, any t>0, and any r>0,
	
    ℙ[  v' (1/' - Γ_(0) ) v ≥ t (B)^2 Σ_ϵ_2  ]
    	≤ 2exp( -Σ_ϵ_2 /c_a c_ϵ^2 min{ t^2,  t  })

and 
	
    ℙ[  v' ((1/' - Γ_(0))⊙ I ) v ≥ r v_0 (B)^2 Σ_ϵ_2   ] 
    	≤ 2exp( -v_0 Σ_ϵ_2 /c_a c_ϵ^2 min{ r^2,  r  })
	
	where c_a is a universal constant. 



The proof is by plugging the bounds of Lemma <ref> in Theorem 2.5 of <cit.> followed by some algebraic manipulations. For the first bound, we bound the operator norm by a scaled Frobenius norm via Ψ_(1)'Ψ_(1)_F^2 ≤Ψ_(1)_2^4. For the second bound, we use

    Ψ_(2)'Ψ_(2)_F^2 
       ≤rank(Ψ_(2)'Ψ_(2)) Ψ_(2)'Ψ_(2)_2^2
       = rank(Ψ_(2)) Ψ_(2)_2^4 
       ≤rank(I_n ⊗diag(v)) Ψ_(2)_2^4
       = v_0 ·Ψ_(2)_2^4 
     .



As it is evident from (<ref>) and (<ref>), they can be directly used to bound the quadratic forms in both terms in (<ref>) and in the first term in (<ref>). For the last term in (<ref>), we can derive another concentration result from (<ref>). Observe that 

    2u' (1/' - Γ_(1) ) v    = 2/ u''v - 2u'Γ_(1)v 
       = 1/ ('u + 'v)' ('u + 'v) 
       - 	[ u' v' ][  Γ_(0)  Γ_(1); Γ_(1)'  Γ_(0) ][ u; v ]
       - ( 1/ u''u - u'Γ_(0)u )
    - ( 1/ v''v - v'Γ_(0)v )  .

Remember W̅ = [ _0  ⋯ _ ] and observe that  and  are simply subsets of this matrix. Hence, u' and v' can be expressed similarly through Ψ_+1(B), and choosing certain rows (corresponding to  and  being subsets of W̅) does not increase the operator norm. 



§ CONCLUSIONS

This paper presented a new methodology and associated performance guarantees for estimating the parameters of linear vector autoregressive processes by leveraging 1) ideas from sparse regression and the LASSO, 2) estimators designed for robustness to missing data, and 3) concentration results from empirical process theory. Note that optimization problems in (<ref>) and (<ref>) are possibly non-convex because of the -B_F^2 term. Without this term we would have a convex formulation, but would not have a consistent estimator. Our approach generalizes to other measurement schemes beyond multiplication by i.i.d. Bernoulli masks as in (<ref>). In fact, we can adapt our analysis to any covariance-stationary observation process independent of the underlying process whose autocovariance matrices of lag 0 and 1 have no zero entries. 

alpha







§ DERIVATION OF THE ESTIMATORS IN (<REF>) AND (<REF>)

In this section, we motivate the design of the proposed estimators in (<ref>) and (<ref>). First, we review the relevant notation. For a p-dimensional discrete-time, centered, covariance-stationary (wide-sense stationary) process {w_t}, denote the autocovariance function by Γ_w(h) = (w_t, w_t+h). For a matrix A, the transpose is denoted by A' and the conjugate transpose is denoted by A^†. Denote by ⊙ and ⊘ the Hadamard (element-wise) product and division respectively, and by ⊗ the Kronecker product. Assuming the process is stationary and ignoring for the moment the fact that Σ_ϵ might not be the identity matrix, for any t, the best linear estimator for  is given by

    B^⋆   = _B  w_t+1 - B w_t_2^2 
       = _B  w_t+1_2^2 
     	+  tr(B w_t w_t' B' ) 
     	- 2  tr(Bw_t w_t+1') 
       = _B  ⟨Γ_w(0) , B'B ⟩ - 2 ⟨Γ_w(1) , B' ⟩

and replacing the autocovariance with its sample approximation yields

    B^⋆≃_B  ⟨1/' , B'B ⟩ - 2 ⟨1/' , B' ⟩ 
     	=_B  1/B  - _F^2  .

Given the prior information that  is sparse, and provided that we have complete information on  and , we can solve either of the following convex optimization problems to estimate :

    _ =_B  1/B  - _F^2  	+ λ_B_1  
     	    or    _ =_B_1≤_1 1/B  - _F^2  	 .

Guarantees on support recovery as well as different error measures for _ with respect to  can be derived through establishing the now well-known lower restricted eigenvalue condition and deviation bound for sample statistics ' and ' <cit.>. 



  
Multiplicative Corruption. When  and  are not fully observed, the above estimators cannot be used anymore. However, going back to (<ref>), we can design a new estimator from scratch if we can estimate the autocovariance matrices Γ_w(0) and Γ_w(1) from the given partial data. Suppose that the underlying process {w_t} is observed through the lens of another covariance-stationary process {m_t}, independent of {w_t}: 

    _t = w_t ⊙ m_t  .

In this case, for any integer value h, we have: 

    Γ_(h) 
    = cov(_t,_t+h)
    = ((w_t ⊙ m_t)(w_t+h⊙ m_t+h)')
    = Γ_w(h) ⊙Γ_m(h)

where we used the fact that w_t = 0 and the independence of w_t and m_t implies _t = 0 regardless of m_t being centered or not. Suppose that for the observation process {m_t}, the autocovariance matrices Γ_m(0) and Γ_m(1) have no zero entries. In this case, we have

    Γ_w(0) = Γ_(0) ⊘Γ_m(0) 
       and   Γ_w(1) = Γ_(0) ⊘Γ_m(1)

which can be plugged in (<ref>) to yield 

    B^⋆   = _B  w_t+1 - B w_t_2^2 
       = _B  ⟨Γ_(0) ⊘Γ_m(0)  , B'B ⟩ - 2 ⟨Γ_(0) ⊘Γ_m(1)  , B' ⟩

whose approximation via the sample autocovariance matrices gives

    B^⋆   ≃_B  ⟨1/' ⊘Γ_m(0)  , B'B ⟩ - 2 ⟨1/' ⊘Γ_m(1) , B' ⟩ .

Observe that the Hadamard division by Γ_m(0) can make the quadratic term non-convex. 



  
Missing Data. 
A simple scenario for partial observations is when each m_t in (<ref>) has entries drawn i.i.d. from a Bernoulli distribution of parameter 1-δ, for some δ∈[0,1). In this case,

    Γ_m(0) = (1-δ)^2 1 + δ(1-δ) I
       and   Γ_m(1) = (1-δ)^2 1

have no zero entries and (<ref>) can be simply expressed as 

    B^⋆   ≃_B  1/tr(B  ( '  - δdiag( ') )  B') 
     	- 2/tr(' B  )	
       = _B  1/B  - _F^2 
    - δBD̅_F^2

where D̅ = (1/diag( '))^1/2∈ℝ^p× p is a diagonal matrix of sample autocovariances for each of the p covariates. In this case, with a possibly non-convex quadratic optimization program, we need a constrained optimization program to hope for recovering the target model. Again, given the prior information that  is sparse, we can use regularization or an ℓ_1-norm constraint. With this consideration, we arrive at the problems in (<ref>) and (<ref>). 




§ ESTIMATION ERROR FOR NON-CONVEX LASSO


Consider two random matrices, a symmetric matrix ∈ℝ^p× p and a matrix ∈ℝ^p× p, as well as a reference matrix ∈ℝ^p× p with _0 = k, and an integer ≥ 1. Suppose the following conditions hold:


  * For any v∈ℝ^p with v_0 ≤ 2 and v_2 = 1, there exists η(s) such that v' ( - )v ≤η() with probability at least 1-p_1(), 


  * There exists φ>0 such that - _∞≤φ with probability at least 1-p_2.

Consider either of the following estimators, 

    ∈_B_1 ≤ b_0 √(k) tr(B  B') - 2⟨ B ,  ⟩ + λB_1 
       ∈_B_1 ≤_1  tr(B  B') - 2⟨ B ,  ⟩ .

Consider the largest value of  that satisfies 

    η() ≤1/27min{Λ_min() ·/ 128 k+ , φ·/b_0√(k)}

while 

    p_3()  p_1() ·exp (2min{log p , log21ep/2}) ≪ 1

and define α_Λ_min() -27 η(). 
Then, for any  with _0 ≤ k, there is a universal positive constant c_0 such that any global optimum  of (<ref>) with any b_0 ≥_F and λ≥ 2φ satisfies the bounds

    - _F ≤c_0 √(k)/α_λ   ,    - _1 ≤8c_0 k/α_λ

with probability at least 1 - p_3(s) -p_2. The same bounds,  where λ is replaced by φ, apply to (<ref>). Further, a threshold variant of (<ref>), defined as B̃ = {_ij1_ |_ij|>λ}_i,j=1,…,p, satisfies 

    supp(B̃) ∖supp()  ≤56c_0 k/α_ .


We omit the proof of Theorem <ref> for brevity. 

Choosing an appropriate  in establishing <ref> might require a lot of algebraic manipulations. Hence, we mention our choice of  in the proof of Theorem <ref> using Theorem <ref>:

    = (1-δ)^2/κ_ϵκ_04hk / 1+4kθ_0√(/log p) ,

where all notations have been defined in Section <ref>. 




§ PROOF OF LEMMA <REF>
 
Using the submultiplicativity of operator norms, we have 

    Ψ_(1)_2 
       ≤(I_n ⊗ v)' I_ΩΨ_(B) _2 
       ≤I_n ⊗ vI_Ω_2 Ψ_(B) _2 
       ≤v_2 Ψ_(B) _2

as well as

    Ψ_(2)_2
       = (I_⊗diag(v) ) I_ΩΨ_(B)   _2 
       = I_Ω(I_⊗diag(v) ) Ψ_(B)   _2 
       ≤(I_⊗diag(v) ) Ψ_(B)   _2.

in which we used the commutativity of diagonal matrices. Notice that we do not bound the last term with Ψ_(B) _2 as a possibly tighter bound is possible. 
It is easy to see that Lemma <ref> implies 

    Ψ_(2)_2 ≤Ψ_(B)_1→2

Therefore, it remains to upper bound Ψ_(B) _2 and Ψ_(B)_1→2. 
However, a closer look reveals that these two quantities are input-output gains, in specific norms, of the following discrete-time linear time-invariant system, 

	x_t+1    = B x_t + u_t 
	 ,       t=0,1,… .

Since we have assumed ρ(B)<1, this system is stable. Moreover, the transfer matrix from u to x is given by 

    G(z) = (zI-B)^-1

where z is a complex number. Therefore, we get the right-most set of inequalities in Lemma <ref>.



Given a matrix A, for any v we have

    diag(v) A_2 
     ≤ v_2  ·I_supp(v)A_1→2
	
where ·_1→ 2 denotes the largest ℓ_2 norm of columns. 


For v∈ℝ^p with v_2=1, observe that 	

    diag(v) A)_2
       = sup_u_2 = 1u' diag(v) A_2
       = sup_u_2 = 1 (u ⊙ v)' I_supp(v)A_2 
       ≤I_supp(v)A_1→2sup_u_2 = 1 u ⊙ v_1   .

Then, u ⊙ v_1 = ∑_i=1^p u_iv_i = u' v≤u_2 ·v_2= 1 establishes the claim.




