Matrix Estimation

Chapter Matrix estimation Over the past decade or so, matrices have entered the picture of high-dimensional statistics for several reasons. Perhaps the simplest explanation is that they are the most natural extension of vectors. While this is true, and we will see examples where the extension from vectors to matrices is straightforward, matrices have a much richer structure than vectors allowing “interaction” between their rows and columns. In particular, while we have been describing simple vectors in terms of their sparsity, here we can measure the complexity of a matrix by its rank. This feature was successfully employed in a variety of applications ranging from multi-task learning to collaborative filtering. This last application was made popular by the Netflix prize in particular. In this chapter, we study several statistical problems where the parameter of interest θ is a matrix rather than a vector. These problems include: multivariate regression, covariance matrix estimation and principal component analysis. Before getting to these topics, we begin by a quick reminder on matrices and linear algebra. 4.1 BASIC FACTS ABOUT MATRICES Matrices are much more complicated objects than vectors. In particular, while vectors can be identified with linear operators from IRd to IR, matrices can be identified to linear operators from IRd to IRn for n ≥ 1. This seemingly simple fact gives rise to a profusion of notions and properties as illustrated by Bernstein’s book [Ber09] that contains facts about matrices over more than a thousand pages. Fortunately, we will be needing only a small number of such properties, which can be found in the excellent book [GVL96], that has become a standard reference on matrices and numerical linear algebra. Singular value decomposition Let A = {aij , 1 ≤ i ≤ m, 1 ≤ j ≤ n} be a m × n real matrix of rank r ≤ min(m, n). The Singular Value Decomposition (SVD) of A is given by A = U DV ⊤ = r X λj uj vj⊤ , j=1 where D is a r × r diagonal matrix with positive diagonal entries {λ1 , . . . , λr }, U is a matrix with columns {u1 , . . . , ur } ∈ IRm that are orthonormal and V is a matrix with columns {v1 , . . . , vr } ∈ IRn that are also orthonormal. Moreover, it holds that AA⊤ uj = λ2j uj , and A⊤ Avj = λ2j vj for j = 1, . . . , r. The values λj > 0 are called singular values of A and are uniquely defined. If rank r < min(n, m) then the singular values of A are given by λ = (λ1 , . . . , λr , 0, . . . , 0)⊤ ∈ IRmin(n,m) where there are min(n, m) − r zeros. This way, the vector λ of singular values of a n × m matrix is a vector in IRmin(n,m) . In particular, if A is a n × n symmetric positive semidefinite (PSD), i.e. A⊤ = A and u⊤ Au ≥ 0 for all u ∈ IRn , then the singular values of A are equal to its eigenvalues. The largest singular value of A denoted by λmax (A) also satisfies the following variational formulation: λmax (A) = maxn x∈IR |Ax|2 y ⊤ Ax = maxn = max y ⊤ Ax . n−1 |x|2 x∈IR |y|2 |x|2 x∈S m y∈IR y∈S m−1 In the case of a n × n PSD matrix A, we have λmax (A) = max x⊤ Ax . n−1 x∈S Norms and inner product Let A = {aij } and B = {bij } be two real matrices. Their size will be implicit in the following notation. Vector norms The simplest way to treat a matrix is to deal with it as if it were a vector. In particular, we can extend ℓq norms to matrices: X 1/q |A|q = |aij |q , q > 0. ij The cases where q ∈ {0, ∞} can also be extended matrices: X |A|0 = 1I(aij = |A|∞ = max |aij | . 6 0) , ij ij The case q = 2 plays a particular role for matrices and |A|2 is called the Frobenius norm of A and is often denoted by kAkF . It is also the HilbertSchmidt norm associated to the inner product: hA, Bi = Tr(A⊤ B) = Tr(B ⊤ A) . Spectral norms Let λ = (λ1 , . . . , λr , 0, . . . , 0) be the singular values of a matrix A. We can define spectral norms on A as vector norms on the vector λ. In particular, for any q ∈ [1, ∞], kAkq = |λ|q , is called Schatten q-norm of A. Here again, special cases have special names: • q = 2: kAk2 = kAkF is the Frobenius norm defined above. • q = 1: kAk1 = kAk∗ is called the Nuclear norm (or trace norm) of A. • q = ∞: kAk∞ = λmax (A) = kAkop is called the operator norm (or spectral norm) of A. We are going to employ these norms to assess the proximity to our matrix of interest. While the interpretation of vector norms is clear by extension from the vector case, the meaning of “kA− Bkop is small” is not as transparent. The following subsection provides some inequalities (without proofs) that allow a better reading. Useful matrix inequalities Let A and B be two m × n matrices with singular values λ1 (A) ≥ λ2 (A) . . . ≥ λmin(m,n) (A) and λ1 (B) ≥ . . . ≥ λmin(m,n) (B) respectively. Then the following inequalities hold: max λk (A) − λk (B) ≤ kA − B kop , k X 2 λk (A) − λk (B) ≤ kA − Bk2F , Weyl (1912) Hoffman-Weilandt (1953) k hA, Bi ≤ kAkq kBkq , 4.2 1 1 + = 1, p, q ∈ [1, ∞] , Hölder p q MULTIVARIATE REGRESSION In the traditional regression setup, the response variable Y is a scalar. In several applications, the goal is not to predict a variable but rather a vector Y ∈ IRT , still from a covariate X ∈ IRd . A standard example arises in genomics data where Y contains T physical measurements of a patient and X contains the expression levels for d genes. As a result the regression function in this case f (x) = IE[Y |X = x] is a function from IRd to IRT . Clearly, f can be estimated independently for each coordinate, using the tools that we have developed in the previous chapter. However, we will see that in several interesting scenarios, some structure is shared across coordinates and this information can be leveraged to yield better prediction bounds. The model Throughout this section, we consider the following multivariate linear regression model: Y = XΘ∗ + E , (4.1) where Y ∈ IRn×T is the matrix of observed responses, X is the n × d observed design matrix (as before), Θ ∈ IRd×T is the matrix of unknown parameters and E ∼ subGn×T (σ 2 ) is the noise matrix. In this chapter, we will focus on the prediction task, which consists in estimating XΘ∗ . As mentioned in the foreword of this chapter, we can view this problem as T (univariate) linear regression problems Y (j) = Xθ∗,(j) +ε(j) , j = 1, . . . , T , where Y (j) , θ∗,(j) and ε(j) are the jth column of Y, Θ∗ and E respectively. In particular, an estimator for XΘ∗ can be obtained by concatenating the estimators for each of the T problems. This approach is the subject of Problem 4.1. The columns of Θ∗ correspond to T different regression tasks. Consider the following example as a motivation. Assume that the Subway headquarters want to evaluate the effect of d variables (promotions, day of the week, TV ads,. . . ) on their sales. To that end, they ask each of their T = 40, 000 restaurants to report their sales numbers for the past n = 200 days. As a result, franchise j returns to headquarters a vector Y(j) ∈ IRn . The d variables for each of the n days are already known to headquarters and are stored in a matrix X ∈ IRn×d . In this case, it may be reasonable to assume that the same subset of variables has an impact of the sales for each of the franchise, though the magnitude of this impact may differ from franchise to franchise. As a result, one may assume that the matrix Θ∗ has each of its T columns that is row sparse and that they share the same sparsity pattern, i.e., Θ∗ is of the form:   0 0 0 0  • • • •     • • • •      Θ∗ =  0 0 0 0  ,  .. .. .. ..   . . . .     0 0 0 0  • • • • where • indicates a potentially nonzero entry. It follows from the result of Problem 4.1 that if each task is performed ˆ such that individually, one may find an estimator Θ 1 kT log(ed) IEkXΘ̂ − XΘ∗ k2F . σ 2 , n n where k is the number of nonzero coordinates in each column of Θ∗ . We remember that the term log(ed) corresponds to the additional price to pay for not knowing where the nonzero components are. However, in this case, when the number of tasks grows, this should become easier. This fact was proved in [LPTVDG11]. We will see that we can recover a similar phenomenon when the number of tasks becomes large, though larger than in [LPTVDG11]. Indeed, rather than exploiting sparsity, observe that such a matrix Θ∗ has rank k. This is the kind of structure that we will be predominantly using in this chapter. Rather than assuming that the columns of Θ∗ share the same sparsity pattern, it may be more appropriate to assume that the matrix Θ∗ is low rank or approximately so. As a result, while the matrix may not be sparse at all, the fact that it is low rank still materializes the idea that some structure is shared across different tasks. In this more general setup, it is assumed that the columns of Θ∗ live in a lower dimensional space. Going back to the Subway example this amounts to assuming that while there are 40,000 franchises, there are only a few canonical profiles for these franchises and that all franchises are linear combinations of these profiles. Sub-Gaussian matrix model Recall that under the assumption ORT for the design matrix, i.e., X⊤ X = nId , then the univariate regression model can be reduced to the sub-Gaussian sequence model. Here we investigate the effect of this assumption on the multivariate regression model (4.1). Observe that under assumption ORT, 1 ⊤ 1 X Y = Θ ∗ + X⊤ E . n n Which can be written as an equation in IRd×T called the sub-Gaussian matrix model (sGMM): y = Θ∗ + F , (4.2) where y = n1 X⊤ Y and F = n1 X⊤ E ∼ subGd×T (σ 2 /n). Indeed, for any u ∈ S d−1 , v ∈ S T −1 , it holds u⊤ F v = 1 1 (Xu)⊤ Ev = √ w⊤ Ev ∼ subG(σ 2 /n) , n n √ ⊤ where w = Xu/ n has unit norm: |w|22 = u⊤ X n X u = |u|22 = 1. Akin to the sub-Gaussian sequence model, we have a direct observation model where we observe the parameter of interest with additive noise. This enables us to use thresholding methods for estimating Θ∗ when |Θ∗ |0 is small. However, this also follows from Problem 4.1. The reduction to the vector case in the sGMM is just as straightforward. The interesting analysis begins when Θ∗ is low-rank, which is equivalent to sparsity in its unknown eigenbasis. Consider the SVD of Θ∗ : X Θ∗ = λj uj vj⊤ . j and recall that kΘ∗ k0 = |λ|0 . Therefore, if we knew uj and vj , we could simply estimate the λj s by hard thresholding. It turns out that estimating these eigenvectors by the eigenvectors of y is sufficient. Consider the SVD of the observed matrix y: X ˆ j ûj v̂ ⊤ . y= λ j j Definition 4.1. The singular value thresholding estimator with threshold 2τ ≥ 0 is defined by X λ̂j 1I(|λ̂j | > 2τ )ûj v̂j⊤ . Θ̂svt = j Recall that the threshold for the hard thresholding estimator was chosen to be the level of the noise with high probability. The singular value thresholding estimator obeys the same rule, except that the norm in which the magnitude of the noise is measured is adapted to the matrix case. Specifically, the following lemma will allow us to control the operator norm of the matrix F . Lemma 4.2. Let A be a d × T random matrix such that A ∼ subGd×T (σ 2 ). Then p p kAkop ≤ 4σ log(12)(d ∨ T ) + 2σ 2 log(1/δ) with probability 1 − δ. Proof. This proof follows the same steps as Problem 1.4. Let N1 be a 1/4net for S d−1 and N2 be a 1/4-net for S T −1 . It follows from Lemma 1.18 that we can always choose |N1 | ≤ 12d and |N2 | ≤ 12T . Moreover, for any u ∈ S d−1 , v ∈ S T −1 , it holds 1 max u⊤ Av x∈N1 4 u∈S d−1 1 1 ≤ max max x⊤ Ay + max max x⊤ Av + max u⊤ Av x∈N1 y∈N2 4 x∈N1 v∈S T −1 4 u∈S d−1 1 ≤ max max x⊤ Ay + max max u⊤ Av x∈N1 y∈N2 2 u∈S d−1 v∈S T −1 u⊤ Av ≤ max x⊤ Av + It yields kAkop ≤ 2 max max x⊤ Ay x∈N1 y∈N2 So that for any t ≥ 0, by a union bound, X

IP kAkop > t ≤ IP x⊤ Ay > t/2 x∈N1 y∈N2 Next, since A ∼ subGd×T (σ 2 ), it holds that x⊤ Ay ∼ subG(σ 2 ) for any x ∈ N1 , y ∈ N2 .Together with the above display, it yields
t2
IP kAkop > t ≤ 12d+T exp − 2 ≤ δ 8σ for t ≥ 4σ p p log(12)(d ∨ T ) + 2σ 2 log(1/δ) . The following theorem holds. Theorem 4.3. Consider the multivariate linear regression model (4.1) under the assumption ORT or, equivalently, the sub-Gaussian matrix model (4.2). ˆ svt with threshold Then, the singular value thresholding estimator Θ r r log(12)(d ∨ T ) 2 log(1/δ) 2τ = 8σ + 4σ , (4.3) n n satisfies 1 kXΘ̂svt − XΘ∗ k2F = kΘ̂svt − Θ∗ k2F ≤ 144 rank(Θ∗ )τ 2 n  σ 2 rank(Θ∗ )  . d ∨ T + log(1/δ) . n with probability 1 − δ. Proof. Assume without loss of generality that the singular values of Θ∗ and y ˆ1 ≥ λ ˆ2 ≥ . . . . are arranged in a non increasing order: λ1 ≥ λ2 ≥ . . . and λ ˆ j | > 2τ }. Define the set S = {j : |λ Observe first that it follows from Lemma 4.2 that kF kop ≤ τ for τ chosen as in (4.3) on an event A such that IP(A) ≥ 1 − δ. The rest of the proof is on A. ˆ j − λj | ≤ kF kop ≤ τ . It Note that it follows from Weyl’s inequality that |λ implies that S ⊂ {j : |λj | > τ }P and S c ⊂ {j : |λj | ≤ 3τ }. ⊤ ¯ = Next define the oracle Θ j∈S λj uj vj and note that ˆ svt − Θk ¯ 2F + 2kΘ ¯ − Θ∗ k2F kΘ̂svt − Θ∗ k2F ≤ 2kΘ Using Cauchy-Schwarz, we control the first term as follows ˆ svt − Θ̄)kΘ̂svt − Θ̄k2 ≤ 2|S|kΘ ˆ svt − Θ̄k2 kΘ̂svt − Θ̄k2F ≤ rank(Θ op op (4.4) Moreover, ˆ svt − Θ ¯ kop ≤ kΘ ˆ svt − ykop + ky − Θ∗ kop + kΘ∗ − Θk ¯ op kΘ ≤ maxc |λ̂j | + τ + maxc |λj | ≤ 6τ . j∈S j∈S Therefore, kΘ̂svt − Θ̄k2F ≤ 72|S|τ 2 = 72 X τ2 . j ∈S The second term in (4.4) can be written as X kΘ̄ − Θ∗ k2F = |λj |2 . j ∈S c Plugging the above two displays in (4.4), we get X X kΘ̂svt − Θ∗ k2F ≤ 144 τ2 + |λj |2 j ∈S j ∈S c Since on S, τ 2 = min(τ 2 , |λj |2 ) and on S c , |λj |2 ≤ 3 min(τ 2 , |λj |2 ), it yields, X kΘ̂svt − Θ∗ k2F ≤ 432 min(τ 2 , |λj |2 ) j rank(Θ∗ ) ≤ 432 X τ2 j=1 = 432 rank(Θ∗ )τ 2 . In the next subsection, we extend our analysis to the case where X does not necessarily satisfy the assumption ORT. Penalization by rank The estimator from this section is the counterpart of the BIC estimator in the spectral domain. However, we will see that unlike BIC, it can be computed efficiently. ˆ rk be any solution to the following minimization problem: Let Θ min Θ∈IRd×T n1 n o kY − XΘk2F + 2τ 2 rank(Θ) . This estimator is called estimator by rank penalization with regularization parameter τ 2 . It enjoys the following property. Theorem 4.4. Consider the multivariate linear regression model (4.1). Then, ˆ rk with regularization parameter τ 2 , where the estimator by rank penalization Θ τ is defined in (4.3) satisfies  1 σ 2 rank(Θ∗ )  2 kXΘ̂rk − XΘ∗ kF ≤ 8 rank(Θ∗ )τ 2 . d ∨ T + log(1/δ) . n n with probability 1 − δ. Proof. We begin as usual by noting that ˆ rk k2F + 2nτ 2 rank(Θ ˆ rk ) ≤ kY − XΘ∗ k2F + 2nτ 2 rank(Θ∗ ) , kY − XΘ which is equivalent to ˆ rk ) + 2nτ 2 rank(Θ∗ ) . ˆ rk − XΘ∗ k2F ≤ 2hE, XΘ ˆ rk − XΘ∗ i − 2nτ 2 rank(Θ k XΘ Next, by Young’s inequality, we have 1 ˆ rk 2hE, XΘ̂rk − XΘ∗ i = 2hE, U i2 + kXΘ − XΘ∗ k2F , 2 where U= Write XΘ̂rk − XΘ∗ kXΘ̂rk − XΘ∗ kF . XΘ̂rk − XΘ∗ = ΦN , where Φ is a n × r, r ≤ d matrix whose columns form orthonormal basis of the column span of X. The matrix Φ can come from the SVD of X for example: X = ΦΛΨ⊤ . It yields ΦN U= kN kF and ˆ rk ) + 4nτ 2 rank(Θ∗ ) . kXΘ̂rk − XΘ∗ k2F ≤ 4hΦ⊤ E, N/kN kF i2 − 4nτ 2 rank(Θ (4.5) ˆ rk ) + rank(Θ∗ ). Therefore, by Hölder’s inNote that rank(N ) ≤ rank(Θ equality, we get hE, U i2 = hΦ⊤ E, N/kN kF i2 ≤ kΦ⊤ Ek2op kN k21 kN k2F ≤ rank(N )kΦ⊤ Ek2op   ˆ rk ) + rank(Θ∗ ) . ≤ kΦ⊤ Ek2 rank(Θ op Next, note that Lemma 4.2 yields kΦ⊤ Ek2op ≤ nτ 2 so that   ˆ rk ) + rank(Θ∗ ) . hE, U i2 ≤ nτ 2 rank(Θ Together with (4.5), this completes the proof. It follows from Theorem 4.4 that the estimator by rank penalization enjoys the same properties as the singular value thresholding estimator even when X does not satisfies the ORT condition. This is reminiscent of the BIC estimator which enjoys the same properties as the hard thresholding estimator. However this analogy does not extend to computational questions. Indeed, while the rank penalty, just like the sparsity penalty, is not convex, it turns out that XΘ̂rk can be computed efficiently. Note first that n1 o 1 min kY − XΘk2F + 2τ 2 rank(Θ) = min min kY − XΘk2F + 2τ 2 k . k n Θ∈IRd×T Θ∈IRd×T n rank(Θ)≤k Therefore, it remains to show that min Θ∈IRd×T rank(Θ)≤k kY − XΘk2F ¯ = X(X⊤ X)† X⊤ Y denote the orthogcan be solved efficiently. To that end, let Y onal projection of Y onto the image space of X: this is a linear operator from IRd×T into IRn×T . By the Pythagorean theorem, we get for any Θ ∈ IRd×T , 2 2 kY − XΘkF = kY − ȲkF + kȲ − XΘk2F . ¯ Next consider the SVD of Y: Ȳ = X λj uj vj⊤ j ˜ by where λ1 ≥ λ2 ≥ . . .. The claim is that if we define Y Ỹ = k X λj uj vj⊤ j=1 which is clearly of rank at most k, then it satisfies 2 kȲ − ỸkF = min Z:rank(Z)≤k Indeed, kȲ − Ỹk2F = X ¯ − Z k2 . kY F λ2j , j>k and for any matrix Z such that rank(Z) ≤ k with SVD Z= k X µj xj yj⊤ , j=1 where µ1 ≥ µ2 ≥ . . ., we have by Hoffman-Weilandt X X kZ − Ȳ k2F ≥ |λj − µj |2 ≥ λ2j . j≥1 j>k Therefore, any minimizer of XΘ 7→ kY − XΘk2F over matrices of rank at most ¯ at order k. k can be obtained by truncating the SVD of Y rk ˆ ˆ rk by least Once XΘ has been found, one may obtain a corresponding Θ squares but this is not necessary for our results. Remark 4.5. While the rank penalized estimator can be computed efficiently, it is worth pointing out that a convex relaxation for the rank penalty can also ˆ is defined to be any be used. The estimator by nuclear norm penalization Θ solution to the minimization problem n1 o min kY − XΘk2F + τ kΘk1 n Θ∈IRd×T Clearly this criterion is convex and it can actually be implemented efficiently using semi-definite programming. It has been popularized by matrix completion problems. Let X have the following SVD: X= r X λj uj vj⊤ , j=1 with λ1 ≥ λ2 ≥ . . . ≥ λr > 0. It can be shown that for some appropriate choice of τ , it holds 1 λ1 σ 2 rank(Θ∗ ) 2 ˆ − XΘ ∗ k F kXΘ . d∨T n λr n with probability .99. However, the proof of this result is far more involved than a simple adaption of the proof for the Lasso estimator to the matrix case (the readers are invited to see that for themselves). For one thing, there is no assumption on the design matrix (such as INC for example). This result can be found in [KLT11]. 4.3 COVARIANCE MATRIX ESTIMATION Empirical covariance matrix Let X1 , .. . , Xn be n i.i.d. copies of a random vector X ∈ IRd such that IE XX ⊤ = Σ for some unknown matrix Σ ≻ 0 called covariance matrix. This matrix contains information about the moments of order 2 of the random vector X. A natural candidate to estimate Σ is the empirical covariance matrix Σ̂ defined by n 1X Σ̂ = Xi Xi⊤ . n i=1 Using the tools of Chapter 1, we can prove the following result. Theorem 4.6. Let X1 , . . . , Xn be n i.i.d. sub-Gaussian random vectors such that IE[XX ⊤] = Σ and X ∼ subGd (kΣkop ). Then r  d + log(1/δ) d + log(1/δ)  kΣ̂ − Σkop . kΣkop ∨ , n n with probability 1 − δ. Proof. Observe first that without  loss of generality we can assume that Σ = Id . Indeed, note that since IE XX ⊤ = Σ ≻ 0, then X ∼ subGd (kΣkop ). Moreover, Y = Σ−1/2 X ∼ subGd (1) and IE[Y Y ⊤ ] = Σ−1/2 ΣΣ−1/2 = Id . Therefore, Pn k1 Xi Xi⊤ − Σkop kΣ̂ − Σkop = n i=1 kΣkop kΣkop P 1/2 kΣ kop k n1 ni=1 Yi Yi⊤ − Id kop kΣ1/2 kop ≤ kΣkop n 1X =k Yi Yi⊤ − Id kop . n i=1 Let N be a 1/4-net for S d−1 such that |N | ≤ 12d. It follows from the proof of Lemma 4.2 that ˆ − Id )y kΣ̂ − Id kop ≤ 2 max x⊤ (Σ x,y∈N So that for any t ≥ 0, by a union bound, X

ˆ − Id kop > t ≤ ˆ − Id )y > t/2 . IP kΣ IP x⊤ (Σ (4.6) x,y ∈N It holds, n X    ⊤ ˆ − Id )y = 1 x⊤ (Σ (Xi x)(Xi⊤ y) − IE (Xi⊤ x)(Xi⊤ y) . n i=1 Using polarization, we also have (Xi⊤ x)(Xi⊤ y) = 2 2 − Z− Z+ , 4 here Z+ = Xi⊤ (x + y) and Z− = Xi⊤ (x − y). It yields  

i  IE exp s (Xi⊤ x)(Xi⊤ y) − IE (Xi⊤ x)(Xi⊤ y)
s 2

i  s 2 2 2 − IE[Z− = IE exp Z+ − IE[Z+ ] − Z− ]) 4 4  

 

1/2 s 2 s 2 2 2 ≤ IE exp Z+ − IE[Z+ ] IE exp − Z− − IE[Z− ] , 2 2 where in the last inequality, we used Cauchy-Schwarz. Next, since X ∼ subGd (1), we have Z+ , Z− ∼ subG(2), and it follows from Lemma 1.12 that 2 2 Z+ − IE[Z+ ] ∼ subE(32) , and 2 2 − IE[Z− ] ∼ subE(32) Z− Therefore for any s ≤ 1/16, we have for any Z ∈ {Z+ , Z− }, we have 

 2 s 2 IE exp Z − IE[Z 2 ] ≤ e128s , 2 It yields that   (Xi⊤ x)(Xi⊤ y) − IE (Xi⊤ x)(Xi⊤ y) ∼ subE(16) . Applying now Bernstein’s inequality (Theorem 1.13), we get h  2
t i ˆ − Id )y > t/2 ≤ exp − n ( t IP x⊤ (Σ ∧ ) . 2 32 32 Together with (4.6), this yields h  2
t i ˆ − Id kop > t ≤ 144d exp − n ( t IP kΣ ∧ ) . 2 32 32 (4.7) In particular, the right hand side of the above inequality is at most δ ∈ (0, 1) if  2d   2d 1/2 t 2 2 ≥ log(144) + log(1/δ) ∨ log(144) + log(1/δ) 32 n n n n This concludes our proof. Theorem 4.6 indicates that for fixed d, the empirical covariance matrix is a consistent estimator of Σ (in any norm as they are all equivalent in finite dimension). However, the bound that we got is not satisfactory in high-dimensions when d ≫ n. To overcome this limitation, we can introduce sparsity as we have done in the case of regression. The most obvious way to do so is to assume that few of the entries of Σ are non zero and it turns out that in this case thresholding is optimal. There is a long line of work on this subject (see for example [CZZ10] and [CZ12]). Once we have a good estimator of Σ, what can we do with it? The key insight is that Σ contains information about the projection of the vector X onto any direction u ∈ S d−1 . Indeed, we have that var(X ⊤ u) = u⊤ Σu, which d ˆ Observe that it follows from can be readily estimated by V ar(X ⊤ u) = u⊤ Σu. Theorem 4.6 d ˆ − Σ)u V ar(X ⊤ u) − Var(X ⊤ u) = u⊤ (Σ ≤ kΣ̂ − Σkop r  d + log(1/δ) d + log(1/δ)  . kΣkop ∨ n n with probability 1 − δ. The above fact is useful in the Markowitz theory of portfolio section for example [Mar52], where a portfolio of assets is a vector u ∈ IRd such that |u|1 = 1 and the risk of a portfolio is given by the variance Var(X ⊤ u). The goal is then to maximize reward subject to risk constraints. In most instances, the empirical covariance matrix is plugged into the formula in place of Σ. 4.4 PRINCIPAL COMPONENT ANALYSIS Spiked covariance model Estimating the variance in all directions is also useful for dimension reduction. In Principal Component Analysis (PCA), the goal is to find one (or more) directions onto which the data X1 , . . . , Xn can be projected without loosing much of its properties. There are several goals for doing this but perhaps the most prominent ones are data visualization (in few dimensions, one can plot and visualize the cloud of n points) and clustering (clustering is a hard computational problem and it is therefore preferable to carry it out in lower dimensions). An example of the output of a principal component analysis is given in Figure 4.1. In this figure, the data has been projected onto two orthogonal directions PC1 and PC2, that were estimated to have the most variance (among all such orthogonal pairs). The idea is that when projected onto such directions, points will remain far apart and a clustering pattern will still emerge. This is the case in Figure 4.1 where the original data is given by d = 500, 000 gene expression levels measured on n ≃ 1, 387 people. Depicted are the projections of these 1, 387 points in two dimension. This image has become quite popular as it shows that gene expression levels can recover the structure induced by geographic clustering. How is it possible to “compress” half a million dimensions into only two? The answer is that the data is intrinsically low dimensional. In this case, a plausible assumption is that all the 1, 387 points live close to a two-dimensional linear subspace. To see how this assumption (in one dimension instead of two for simplicity) translates into the structure of the covariance matrix Σ, assume that X1 , . . . , Xn are Gaussian random variables generated as follows. Fix a direction v ∈ S d−1 and let Y1 , . . . , Yn ∼ Nd (0, Id ) so that v ⊤ Yi are i.i.d. N (0, 1). In particular, the vectors (v ⊤ Y1 )v, . . . , (v ⊤ Yn )v live in the one-dimensional space spanned by v. If one would observe such data the problem would be easy as only two observations would suffice to recover v. Instead, we observe X1 , . . . , Xn ∈ IRd where Xi = (v ⊤ Yi )v + Zi , and Zi ∼ Nd (0, σ 2 Id ) are i.i.d. and independent of the Yi s, that is we add a isotropic noise to every point. If the σ is small enough, we can hope to recover the direction v (See Figure 4.2). The covariance matrix of Xi generated as such is given by     Σ = IE XX ⊤ = IE ((v ⊤ Y )v + Z)((v ⊤ Y )v + Z)⊤ = vv ⊤ + σ 2 Id . This model is often called the spiked covariance model. By a simple rescaling, it is equivalent to the following definition. Definition 4.7. A covariance matrix Σ ∈ IRd×d is said to satisfy the spiked covariance model if it is of the form Σ = θvv ⊤ + Id , where θ > 0 and v ∈ S d−1 . The vector v is called the spike. Courtesy of Macmillan Publishers Ltd. Used with permission. Figure 4.1. Projection onto two dimensions of 1, 387 points from gene expression data. Source: Gene expression blog. v Figure 4.2. Points are close to a one dimensional space space by v. This model can be extended to more than one spike but this extension is beyond the scope of these notes. Clearly, under the spiked covariance model, v is the eigenvector of the matrix Σ that is associated to its largest eigenvalue 1 + θ. We will refer to this vector simply as largest eigenvector. To estimate it, a natural candidate ˜ where Σ ˜ is any estimator of Σ. There is a is the largest eigenvector v̂ of Σ, caveat: by symmetry, if u is an eigenvector, of a symmetric matrix, then −u is also an eigenvector associated to the same eigenvalue. Therefore, we may only estimate v up to a sign flip. To overcome this limitation, it is often useful to describe proximity between two vectors u and v in terms of the principal angle between their linear span. Let us recall that for two unit vectors the principal angle between their linear spans is denoted by ∠(u, v) and defined as ∠(u, v) = arccos(|u⊤ v|) . The following result form perturbation theory is known as the Davis-Kahan sin(θ) theorem as it bounds the sin of the principal angle between eigenspaces. This theorem exists in much more general versions that extend beyond onedimensional eigenspaces. Theorem 4.8 (Davis-Kahan sin(θ) theorem). Let Σ satisfy the spiked covari˜ be any PSD estimator of Σ. Let ṽ denote the largest ance model and let Σ ˜ eigenvector of Σ. Then we have
8 ˜ 2 − Σkop . min |εṽ − v|22 ≤ 2 sin2 ∠(ṽ, v) ≤ 2 kΣ θ ε∈{±1} Proof. Note that for any u ∈ S d−1 , it holds under the spiked covariance model that u⊤ Σu = 1 + θ(v ⊤ u)2 = 1 + θ cos2 (∠(u, v)) . Therefore, v ⊤ Σv − ṽ ⊤ Σṽ = θ[1 − cos2 (∠(ṽ, v))] = θ sin2 (∠(ṽ, v)) . Next, observe that
˜ − ṽ ⊤ Σṽ − v ⊤ Σ ˜ −Σ v v ⊤ Σv − ṽ ⊤ Σṽ = v ⊤ Σv
˜ − ṽ ⊤ Σṽ − v ⊤ Σ ˜ −Σ v ≤ ṽ ⊤ Σṽ = hΣ̂ − Σ, ṽṽ ⊤ − vv ⊤ i ⊤ (4.8) ⊤ ≤ kΣ̃ − Σkop kṽṽ − vv k1 √ ≤ 2kΣ̃ − Σkop kṽṽ ⊤ − vv ⊤ kF (Hölder) (Cauchy-Schwarz) . where in the first inequality, we used the fact that ṽ is the largest eigenvector ˜ and in the last one, we used the fact that the matrix ṽṽ ⊤ − vv ⊤ has rank of Σ at most 2. Next, we have that kṽṽ ⊤ − vv ⊤ k2F = 2(1 − (v ⊤ ṽ)2 ) = 2 sin2 (∠(ṽ, v)) . Therefore, we have proved that ˜ − Σkop sin(∠(ṽ, v)) , θ sin2 (∠(ṽ, v)) ≤ 2kΣ so that 2 ˜ kΣ − Σkop . θ To conclude the proof, it remains to check that sin(∠(ṽ, v)) ≤ min |εṽ − v|22 = 2 − 2|ṽ ⊤ v| ≤ 2 − 2(ṽ ⊤ v)2 = 2 sin2 (∠(ṽ, v)) . ε∈{±1} Combined with Theorem 4.6, we immediately get the following corollary. Corollary 4.9. Let X1 , . . . , X  n be n i.i.d. copies of a sub-Gaussian random vector X ∈ IRd such that IE XX ⊤ = Σ and X ∼ subGd (kΣkop ). Assume further that Σ = θvv ⊤ + Id satisfies the spiked covariance model. Then, the ˆ satisfies, largest eigenvector v̂ of the empirical covariance matrix Σ r 1 + θ  d + log(1/δ) d + log(1/δ)  min |εv̂ − v |2 . ∨ θ n n ε∈{±1} with probability 1 − δ. ˆ as a reThis result justifies the use of the empirical covariance matrix Σ placement for the true covariance matrix Σ when performing PCA in low dimensions, that is when d ≪ n. In the high-dimensional case, where d ≫ n, the above result is uninformative. As before, we resort to sparsity to overcome this limitation. Sparse PCA In the example of Figure 4.1, it may be desirable to interpret the meaning of the two directions denoted by PC1 and PC2. We know that they are linear combinations of the original 500,000 gene expression levels. A natural question to ask is whether only a subset of these genes could suffice to obtain similar results. Such a discovery could have potential interesting scientific applications as it would point to a few genes responsible for disparities between European populations. In the case of the spiked covariance model this amounts to have v to be sparse. Beyond interpretability as we just discussed, sparsity should also lead to statistical stability as in the case of sparse linear regression for example. To enforce sparsity, we will assume that v in the spiked covariance model is k-sparse: |v|0 = k. Therefore, a natural candidate to estimate v is given by v̂ defined by v̂ ⊤ Σ̂v̂ = max u⊤ Σ̂u . u∈S d−1 |u|0 =k ˆ = v̂ ⊤ Σv̂ ˆ is the largest of all leading eigenvalues It is easy to check that λkmax (Σ) ˆ among all k × k sub-matrices of Σ so that the maximum is indeed attained, ˆ the k-sparse leading though there my be several maximizers. We call λkmax (Σ) ˆ eigenvalue of Σ and v̂ a k-sparse leading eigenvector. Theorem 4.10. Let X1 , . . . ,Xn be n i.i.d. copies of a sub-Gaussian random vector X ∈ IRd such that IE XX ⊤ = Σ and X ∼ subGd (kΣkop ). Assume further that Σ = θvv ⊤ + Id satisfies the spiked covariance model for v such that |v|0 = k ≤ d/2. Then, the k-sparse largest eigenvector v̂ of the empirical covariance matrix satisfies, r 1 + θ  k log(ed/k) + log(1/δ) k log(ed/k) + log(1/δ)  ∨ . min |εv̂ − v |2 . θ n n ε∈{±1} with probability 1 − δ. Proof. We begin by obtaining an intermediate result of the Davis-Kahan sin(θ) theorem (Theorem 4.8). Note that we get from (4.8) that ˆ − Σ, v̂v̂ ⊤ − vv ⊤ i v ⊤ Σv − v̂ ⊤ Σv̂ ≤ hΣ Since both v̂ and v are k sparse, there exists a (random) set S ⊂ {1, . . . , d} such that |S| ≤ 2k and {v̂v̂ ⊤ − vv ⊤ }ij = 0 if (i, j) ∈ / S 2 . It yields ˆ hΣ̂ − Σ, v̂v̂ ⊤ − vv ⊤ i = hΣ(S) − Σ(S), v̂(S)v̂(S)⊤ − v(S)v(S)⊤ i Where for any d × d matrix M , we defined the matrix M (S) to be the |S| × |S| sub-matrix of M with rows and columns indexed by S and for any vector x ∈ IRd , x(S) ∈ IR|S| denotes the sub-vector of x with coordinates indexed by S. It yields, by Hölder’s inequality that ˆ v ⊤ Σv − v̂ ⊤ Σv̂ ≤ kΣ(S) − Σ(S)kop kv̂(S)v̂(S)⊤ − v(S)v(S)⊤ k1 . Following the same steps as in the proof of Theorem 4.8, we get now that
8 min |εv̂ − v|22 ≤ 2 sin2 ∠(v̂, v) ≤ 2 θ ε∈{±1} sup S : |S|=2k ˆ kΣ(S) − Σ(S)kop . ˆ To conclude the proof, it remains to control supS : |S |=2k kΣ(S) − Σ(S)kop . To that end, observe that h i ˆ − Σ(S)kop > tkΣkop IP sup kΣ(S) S : |S|=2k ≤ ≤ X S : |S |=2k IP h sup S : |S |=2k kΣ̂(S) − Σ(S)kop > tkΣ(S)kop  h n  t 2 d t i ∧ ) . 1442k exp − ( 2k 2 32 32  i where we used (4.7) in the second inequality. Using Lemma 2.7, we get that the right-hand side above is further bounded by h n  t 2 t ed
i exp − ( ∧ ) + 2k log(144) + k log 2 32 32 2k Choosing now t such that r k log(ed/k) + log(1/δ) k log(ed/k) + log(1/δ) t≥C ∨ , n n for large enough C ensures that the desired bound holds with probability at least 1 − δ.