Linear Regression Model

Chapter 2 Linear Regression Model In this chapter, we consider the following regression model: Yi = f (Xi ) + εi , i = 1, . . . , n , (2.1) where ε = (ε1 , . . . , εn )⊤ is sub-Gaussian with variance proxy σ 2 and such that IE[ε] = 0. Our goal is to estimate the function f under a linear assumption. Namely, we assume that x ∈ IRd and f (x) = x⊤ θ∗ for some unknown θ∗ ∈ IRd . 2.1 FIXED DESIGN LINEAR REGRESSION Depending on the nature of the design points X1 , . . . , Xn , we will favor a different measure of risk. In particular, we will focus either on fixed or random design. Random design The case of random design corresponds to the statistical learning setup. Let (X1 , Y1 ), . . . , (Xn+1 , Yn+1 ) be n+1 i.i.d. random couples. Given (X1 , Y1 ), . . . , (Xn , Yn ) the goal is construct a function fˆn such that fˆn (Xn+1 ) is a good predictor of Yn+1 . Note that when fˆn is constructed, Xn+1 is still unknown and we have to account for what value it is likely to take. Consider the following example from [HTF01, Section 3.2]. The response variable Y is the log-volume of a cancerous tumor, and the goal is to predict it based on X ∈ IR6 , a collection of variables that are easier to measure (age of patient, log-weight of prostate, . . . ). Here the goal is clearly to construct f for prediction purposes. Indeed, we want to find an automatic mechanism that outputs a good prediction of the log-weight of the tumor given certain inputs for a new (unseen) patient. A natural measure of performance here is the L2 -risk employed in the introduction: R(fˆn ) = IE[Yn+1 − fˆn (Xn+1 )]2 = IE[Yn+1 − f (Xn+1 )]2 + kfˆn − f k2L2 (PX ) , where PX denotes the marginal distribution of Xn+1 . It measures how good the prediction of Yn+1 is in average over realizations of Xn+1 . In particular, it does not put much emphasis on values of Xn+1 that are not very likely to occur. Note that if the εi are random variables with variance σ 2 then, one simply has R(fˆn ) = σ 2 + kfˆn − f k2L2 (PX ) . Therefore, for random design, we will focus on the squared L2 norm kfˆn − f k2 2 as a measure of accuracy. It measures L (PX ) how close fˆn is to the unknown f in average over realizations of Xn+1 . Fixed design In fixed design, the points (or vectors) X1 , . . . , Xn are deterministic. To emphasize this fact, we use lowercase letters x1 , . . . , xn to denote fixed design. Of course, we can always think of them as realizations of a random variable but the distinction between fixed and random design is deeper and significantly affects our measure of performance. Indeed, recall that for random design, we look at the performance in average over realizations of Xn+1 . Here, there is no such thing as a marginal distribution of Xn+1 . Rather, since the design points x1 , . . . , xn are considered deterministic, our goal is estimate f only at these points. This problem is sometimes called denoising since our goal is to recover f (x1 ), . . . , f (xn ) given noisy observations of these values. In many instances, fixed design can be recognized from their structured form. A typical example is the regular design on [0, 1], given by xi = i/n, i = 1, . . . , n. Interpolation between these points is possible under smoothness assumptions.
⊤ Note that in fixed design, we observe µ∗ +ε, where µ∗ = f (x1 ), . . . , f (xn ) ∈ IRn and ε = (ε1 , . . . , ε)⊤ ∈ IRn is sub-Gaussian with variance proxy σ 2 . Instead of a functional estimation problem, it is often simpler to view this problem as a vector problem in IRn . This point of view will allow us to leverage the Euclidean geometry of IRn . In the case of fixed design, we will focus on the Mean Squared Error (MSE) as a measure of performance. It is defined by n
2 1X ˆ fn (xi ) − f (xi ) . MSE(fˆn ) = n i=1 Equivalently, if we view our problem as a vector problem, it is defined by n
2 1X 1 MSE(µ̂) = µ̂i − µ∗i = |µ̂ − µ∗ |22 . n i=1 n Often, the design vectors x1 , . . . , xn ∈ IRd are stored in a n × d design matrix X, whose jth row is given by x⊤ j . With this notation, the linear regression model can be written Y = Xθ∗ + ε , (2.2) where Y = (Y1 , . . . , Yn )⊤ and ε = (ε1 , . . . , εn )⊤ . Moreover, MSE(Xθ̂) = 1 X⊤ X ˆ (θ − θ∗ ) . |X(θˆ − θ∗ )|22 = (θˆ − θ∗ )⊤ n n (2.3) A natural example of fixed design regression is image denoising. Assume that µ∗i , i ∈ 1, . . . , n is the grayscale value of pixel i of an image. We do not get to observe the image µ∗ but rather a noisy version of it Y = µ∗ + ε. Given a library of d images {x1 , . . . , xd }, xj ∈ IRn , our goal is to recover the original image µ∗ using linear combinations of the images x1 , . . . , xd . This can be done fairly accurately (see Figure 2.1). Figure 2.1. Reconstruction of the digit “6”: Original (left), Noisy (middle) and Reconstruction (right). Here n = 16 × 16 = 256 pixels. Source [RT11]. As we will see in Remark 2.3, choosing fixed design properly also ensures ˆ then |θˆ − θ∗ |2 is that if MSE(fˆ) is small for some linear estimator fˆ(x) = x⊤ θ, 2 also small. In this chapter we only consider the fixed design case. 2.2 LEAST SQUARES ESTIMATORS Throughout this section, we consider the regression model (2.2) with fixed design. Unconstrained least squares estimator Define the (unconstrained) least squares estimator θˆls to be any vector such that θ̂ls ∈ argmin |Y − Xθ|22 . θ∈IRd Note that we are interested in estimating Xθ∗ and not θ∗ itself, so by extension, we also call µ̂ls = Xθˆls least squares estimator. Observe that µ̂ls is the projection of Y onto the column span of X. It is not hard to see that least squares estimators of θ∗ and µ∗ = Xθ∗ are maximum likelihood estimators when ε ∼ N (0, σ 2 In ). Proposition 2.1. The least squares estimator µ̂ls = Xθˆls ∈ IRn satisfies X⊤ µ̂ls = X⊤ Y . Moreover, θˆls can be chosen to be θ̂ls = (X⊤ X)† X⊤ Y , where (X⊤ X)† denotes the Moore-Penrose pseudoinverse of X⊤ X. Proof. The function θ 7→ |Y − Xθ|22 is convex so any of its minima satisfies ∇θ |Y − Xθ|22 = 0 Where ∇θ is the gradient operator. Using matrix calculus, we find  ∇θ |Y − Xθ|22 = ∇θ |Y |22 + −2Y ⊤ Xθ + θ⊤ X⊤ Xθ = −2(Y ⊤ X − θ⊤ X⊤ X)⊤ . Therefore, solving ∇θ |Y − Xθ|22 = 0 yields X⊤ Xθ = X⊤ Y . It concludes the proof of the first statement. The second statement follows from the definition of the Moore-Penrose pseudoinverse. We are now going to prove our first result on the finite sample performance of the least squares estimator for fixed design. Theorem 2.2. Assume that the linear model (2.2) holds where ε ∼ subGn (σ 2 ). Then the least squares estimator θˆls satisfies   1 r IE MSE(Xθˆls ) = IE|Xθˆls − Xθ∗ |22 . σ 2 , n n where r = rank(X⊤ X). Moreover, for any δ > 0, with probability 1 − δ, it holds MSE(Xθ̂ls ) . σ 2 r + log(1/δ) . n Proof. Note that by definition |Y − Xθ̂ls |22 ≤ |Y − Xθ∗ |22 = |ε|22 . (2.4) Moreover, |Y − Xθ̂ls |22 = |Xθ∗ + ε − Xθ̂ls |22 = |Xθ̂ls − Xθ∗ |22 − 2ε⊤ X(θˆls − θ∗ ) + |ε|22 . Therefore, we get ε⊤ X(θˆls − θ∗ ) |Xθ̂ls − Xθ∗ |22 ≤ 2ε⊤ X(θˆls − θ∗ ) = 2|Xθˆls − Xθ∗ |2 |X(θˆls − θ∗ )|2 (2.5) Note that it is difficult to control ε⊤ X(θˆls − θ∗ ) |X(θˆls − θ∗ )|2 as θˆls depends on ε and the dependence structure of this term may be complicated. To remove this dependency, a traditional technique is “sup-out” θˆls . This is typically where maximal inequalities are needed. Here we have to be a bit careful. Let Φ = [φ1 , . . . , φr ] ∈ IRn×r be an orthonormal basis of the column span of X. In particular, there exists ν ∈ IRr such that X(θˆls − θ∗ ) = Φν. It yields ε⊤ X(θˆls − θ∗ ) ε⊤ Φν ε⊤ Φν ν = = = ε̃⊤ ≤ sup ε̃⊤ u , |Φν|2 |ν|2 |ν|2 u∈B2 |X(θˆls − θ∗ )|2 where B2 is the unit ball of IRr and ε̃ = Φ⊤ ε. Thus |Xθ̂ls − Xθ∗ |22 ≤ 4 sup (ε̃⊤ u)2 , u∈B2 Next, note that for any u ∈ S r−1 , it holds |Φu|22 = u⊤ Φ⊤ Φu = u⊤ u = 1 so that for any s ∈ IR, we have ⊤ IE[esε̃ u ⊤ ] = IE[esε Φu ]≤e s2 σ 2 2 . Therefore, ε̃ ∼ subGr (σ 2 ). To conclude the bound in expectation, observe that Lemma 1.4 yields r X   4IE sup (ε̃⊤ u)2 = 4 IE[ε̃2i ] ≤ 16σ 2 r . u∈B2 i=1 Moreover, with probability 1 − δ, it follows from the last step in the proof1 of Theorem 1.19 that sup (ε̃⊤ u)2 ≤ 8 log(6)σ 2 r + 8σ 2 log(1/δ) . u∈B2 Remark 2.3. If d ≤ n and B := X⊤ X n |θˆls − θ∗ |22 ≤ has rank d, then we have MSE(Xθ̂ls ) , λmin (B) and we can use Theorem 2.2 to bound |θˆls − θ∗ |22 directly. 1 we could use Theorem 1.19 directly here but at the cost of a factor 2 in the constant. Constrained least squares estimator Let K ⊂ IRd be a symmetric convex set. If we know a priori that θ∗ ∈ K, we ls may prefer a constrained least squares estimator θˆK defined by ls θ̂K ∈ argmin |Y − Xθ|22 . θ ∈K Indeed, the fundamental inequality (2.4) would still hold and the bounds on the MSE may be smaller. Indeed, (2.5) can be replaced by ls ls |XθˆK − Xθ∗ |22 ≤ 2ε⊤ X(θˆK − θ∗ ) ≤ 2 sup (ε⊤ Xθ) , θ∈K−K where K − K = {x − y : x, y ∈ K}. It is easy to see that if K is symmetric (x ∈ K ⇔ −x ∈ K) and convex, then K − K = 2K so that 2 sup (ε⊤ Xθ) = 4 sup (ε⊤ v) v∈XK θ∈K−K where XK = {Xθ : θ ∈ K} ⊂ IRn . This is a measure of the size (width) of XK. If ε ∼ N (0, Id ), the expected value of the above supremum is actually called Gaussian width of XK. Here, ε is not Gaussian but sub-Gaussian and similar properties will hold. ℓ1 constrained least squares Assume here that K = B1 is the unit ℓ1 ball of IRd . Recall that it is defined by d n X B1 = x ∈ IRd : |xi | ≤ 1 , i=1 and it has exactly 2d vertices V = {e1 , −e1 , . . . , ed , −ed }, where ej is the j-th vector of the canonical basis of IRd and is defined by ej = (0, . . . , 0, |{1z} , 0, . . . , 0)⊤ . jth position It implies that the set XK = {Xθ, θ ∈ K } ⊂ IRn is also a polytope with at most 2d vertices that are in the set XV = {X1 , −X1 , . . . , Xd , −Xd } where Xj is the j-th column of X. Indeed, XK is a obtained by rescaling and embedding (resp. projecting) the polytope K when d ≤ n (resp., d ≥ n). Note that some columns of X might not be vertices of XK so that XV might be a strict superset of the set of vertices of XK. Theorem 2.4. Let K = B1 be the unit ℓ1 ball of IRd , d ≥ 2 and assume that θ∗ ∈ B1 . Moreover, assume the conditions of Theorem√2.2 and that the columns of X are normalized in such a way that maxj |Xj |2 ≤ n. Then the constrained ls least squares estimator θˆB satisfies 1 r   1 log d ls ls ∗ 2 ˆ ˆ IE MSE(XθB1 ) = IE|XθB1 − Xθ |2 . σ , n n Moreover, for any δ > 0, with probability 1 − δ, it holds r log(d/δ) ls MSE(Xθ̂B1 ) . σ . n Proof. From the considerations preceding the theorem, we got that ls |Xθ̂B − Xθ∗ |22 ≤ 4 sup (ε⊤ v) 1 v∈XK Observe √ now that since ε ∼ subGn (σ 2 ),then for any column Xj such that |Xj |2 ≤ n, the random variable ε⊤ Xj ∼ subG(nσ 2 ). Therefore, applying   ls Theorem 1.16, we get the bound on IE MSE(XθˆK ) and for any t ≥ 0,     nt2 ls IP MSE(XθˆK ) > t ≤ IP sup (ε⊤ v) > nt/4 ≤ 2de− 32σ2 v ∈XK To conclude the proof, we find t such that nt2 2de− 32σ2 ≤ δ ⇔ t2 ≥ 32σ 2 log(2d) log(1/δ) + 32σ 2 . n n ls (exercise!) so that Note that the proof of Theorem 2.2 also applies to θˆB 1 ls θ̂B1 benefits from the best of both rates.  r r log d    ls ˆ IE MSE(XθB1 ) . min , . n n √ This is called an elbow effect. The elbow takes place around r ≃ n (up to logarithmic terms). ℓ0 constrained least squares We abusively call ℓ0 norm of a vector θ ∈ IRd it number of non-zero coefficient. It is denoted by d X |θ|0 = 1I(θj = 6 0) . j=1 We call a vector θ with “small” ℓ0 norm a sparse vector. More precisely, if |θ|0 ≤ k, we say that θ is a k-sparse vector. We also call support of θ the set  supp(θ) = j ∈ {1, . . . , d} : θj = 6 0 so that |θ|0 = card(supp(θ)) =: | supp(θ)| . Remark 2.5. The ℓ0 terminology and notation comes from the fact that lim+ q→0 d X j=1 |θj |q = |θ|0 Therefore it is really limq→0+ |θ|qq but the notation |θ|00 suggests too much that it is always equal to 1. By extension, denote by B0 (k) the ℓ0 ball of IRd , i.e., the set of k-sparse vectors, defined by B0 (k) = {θ ∈ IRd : |θ|0 ≤ k} . ls In this section, our goal is to control the MSE of θˆK when K = B0 (k). Note that
d ls ˆ computing θ essentially requires computing least squares estimators, B0 (k) k which is an exponential number in k. In practice this will be hard (or even impossible) but it is interesting to understand the statistical properties of this estimator and to use them as a benchmark. Theorem 2.6. Fix a positive integer k ≤ d/2. Let K = B0 (k) be set of k-sparse vectors of IRd and assume that θ∗ ∈ B0 (k). Moreover, assume the conditions of Theorem 2.2. Then, for any δ > 0, with probability 1 − δ, it holds   d σ2 k σ2 σ2 ls ˆ MSE(XθB0 (k) ) . log + + log(1/δ) . n 2k n n Proof. We begin as in the proof of Theorem 2.2 to get (2.5): ls ε⊤ X(θˆK − θ∗ ) ls ls ls |Xθ̂K − Xθ∗ |22 ≤ 2ε⊤ X(θˆK − θ∗ ) = 2|Xθ̂K − X θ ∗ |2 . ls − θ ∗ )| |X(θˆK 2 ls ls We know that both θˆK and θ∗ are in B0 (k) so that θˆK − θ∗ ∈ B0 (2k). For any S ⊂ {1, . . . , d}, let XS denote the n × |S| submatrix of X that is obtained from the columns of Xj , j ∈ S of X. Denote by rS ≤ |S| the rank of XS and let ΦS = [φ1 , . . . , φrS ] ∈ IRn×rS be an orthonormal basis of the column span of XS . Moreover, for any θ ∈ IRd , define θ(S) ∈ IR|S| to be the vector with ls ˆ ≤ 2k coordinates θj , j ∈ S. If we denote by Sˆ = supp(θˆK − θ∗ ), we have |S| rŜ and there exists ν ∈ IR such that ls ls ˆ ˆ = Φ ˆν . X(θˆK − θ∗ ) = XSˆ (θˆK (S) − θ∗ (S)) S Therefore, ls ε⊤ ΦSˆ ν ε⊤ X(θˆK − θ∗ ) = ≤ max sup [ε⊤ ΦS ]u ls − θ ∗ )| |ν | 2 |S|=2k u∈BrS |X(θˆK 2 2 where B2rS is the unit ball of IRrS . It yields ls 2 |Xθ̂K − Xθ∗ |22 ≤ 4 max sup (ε̃⊤ S u) , |S|=2k u∈BrS 2 2 ε̃S = Φ⊤ S ε ∼ subGrS (σ ). Using a union bound, we get for any t > 0, X
IP max sup (ε̃⊤ u)2 > t ≤ IP |S|=2k u∈BrS 2 |S |=2k sup (ε̃⊤ u)2 > t r u∈B2S It follows from the proof of Theorem 1.19 that for any |S| ≤ 2k,
t t IP sup (ε̃⊤ u)2 > t ≤ 6|S| e− 8σ2 ≤ 62k e− 8σ2 . r u∈B2S
Together, the above three displays yield ls IP(|XθˆK − Xθ∗ |22 > 4t) ≤   t d 62k e− 8σ2 . 2k (2.6) To ensure that the right-hand side of the above inequality is bounded by δ, we need   n o d + k log(6) + log(1/δ) . t ≥ Cσ 2 log 2k How large is log d 2k
? It turns out that it is not much larger than k. Lemma 2.7. For any integers 1 ≤ k ≤ n, it holds     n en k ≤ k k Proof. Observe first that if k = 1, since n ≥ 1, it holds,    en 1 n = n ≤ en = 1 1 Next, we proceed by induction and assume that it holds for some k ≤ n − 1.     n en k ≤ k k Observe that     n n n − k  en k n − k ek nk+1  1 k = ≤ = 1 + , k+1 k k+1 k k+1 (k + 1)k+1 k where we used the induction hypothesis in the first inequality. To conclude, it suffices to observe that  1 k 1+ ≤e k It immediately leads to the following corollary: Corollary 2.8. Under the assumptions of Theorem 2.6, for any δ > 0, with probability 1 − δ, it holds ls MSE(Xθ̂B ). 0 (k)  ed  σ 2 k σ2 k σ2 log + log(6) + log(1/δ) . n 2k n n Note that for any fixed δ, there exits a constant Cδ > 0 such that for any n ≥ 2k,  ed  σ2 k ls MSE(XθˆB ) ≤ Cδ log . 0 (k) n 2k Comparing this result with Theorem 2.2 with r = k, we see that the price to pay for not knowing the support of θ∗ but only its size, is a logarithmic factor in the dimension d. This result immediately leads the following bound in expectation. Corollary 2.9. Under the assumptions of Theorem 2.6,  ed    σ2 k ls IE MSE(Xθ̂B ) . log . 0 (k) n k Proof. It follows from (2.6) that for any H ≥ 0, Z ∞   ls ls IE MSE(Xθ̂B0 (k) ) = IP(|Xθ̂K − Xθ∗ |22 > nu)du 0 Z ∞ ls ≤H+ IP(|Xθ̂K − Xθ∗ |22 > n(u + H))du 0 Z 2k   X d 2k ∞ − n(u+H) ≤H+ 6 e 32σ2 , j 0 j=1 =H+ 2k   X d 2k − nH2 32σ 2 6 e 32σ du . j n j=1 Next, take H to be such that 2k   X d 2k − nH2 6 e 32σ = 1 . j j=1 In particular, it yields H. which completes the proof.  ed  σ2 k log , n k 2.3 THE GAUSSIAN SEQUENCE MODEL The Gaussian Sequence Model is a toy model that has received a lot of attention, mostly in the eighties. The main reason for its popularity is that it carries already most of the insight of nonparametric estimation. While the model looks very simple it allows to carry deep ideas that extend beyond its framework and in particular to the linear regression model that we are interested in. Unfortunately we will only cover a small part of these ideas and the interested reader should definitely look at the excellent books by A. Tsybakov [Tsy09, Chapter 3] and I. Johnstone [Joh11]. The model is as follows: Yi = θi∗ + εi , i = 1, . . . , d (2.7) 2 where ε1 , . . . , εd are i.i.d N (0, σ ) random variables. Note that often, d is taken equal to ∞ in this sequence model and we will also discuss this case. Its links to nonparametric estimation will become clearer in Chapter 3. The goal here is to estimate the unknown vector θ∗ . The sub-Gaussian Sequence Model Note first that the model (2.7) is a special case of the linear model with fixed design (2.1) with n = d, f (x) = x⊤ θ∗ , x1 , . . . , xn form the canonical basis of IRn and ε has a Gaussian distribution. Therefore, n = d is both the dimension of the parameter θ and the number of observation and it looks like we have chosen to index this problem by d rather than n somewhat arbitrarily. We can bring n back into the picture, by observing that this model encompasses slightly more general choices for the design matrix X as long as it satisfies the following assumption. Assumption ORT The design matrix satisfies X⊤ X = Id , n where Id denotes the identity matrix of IRd . Assumption ORT allows for cases where d ≤ n but not d > n (high dimensional case) because of obvious rank constraints. In particular, √ it means that the d columns of X are orthogonal in IRn and all have norm n. Under this assumption, it follows from the linear regression model (2.2) that y := 1 ⊤ X⊤ X ∗ 1 ⊤ X Y = θ + X ε n n n = θ∗ + ξ , where ξ = (ξ1 , . . . , ξd ) ∼ subGd (σ 2 /n). As a result, under the assumption ORT, the linear regression model (2.2) is equivalent to the sub-Gaussian Sequence Model (2.7) up to a transformation of the data Y and a change of variable for the variance. Moreover, for any estimator θˆ ∈ IRd , under ORT, it follows from (2.3) that ˆ = (θˆ − θ∗ )⊤ MSE(Xθ) X⊤ X ˆ (θ − θ∗ ) = |θˆ − θ∗ |22 . n Furthermore, for any θ ∈ IRd , the assumption ORT yields, 1 |y − θ|22 = | X⊤ Y − θ|22 n 2 1 = |θ|22 − θ⊤ X⊤ Y + 2 Y ⊤ XX⊤ Y n n 1 2 1 = |Xθ|22 − (Xθ)⊤ Y + |Y |22 + Q n n n 1 = |Y − Xθ|22 + Q , n (2.8) where Q is a constant that does not depend on θ and is defined by Q= 1 ⊤ ⊤ 1 Y XX Y − |Y |22 n2 n This implies in particular that the least squares estimator θˆls is equal to y. We introduce a sightly more general model called sub-Gaussian sequence model : y = θ∗ + ξ 2 ∈ IRd (2.9) where ξ ∼ subGd (σ /n). In this section, we can actually completely “forget” about our original model (2.2). In particular we can define this model independently of Assumption ORT and thus for any values of n and d. The sub-Gaussian sequence model, like the Gaussian sequence model are called direct (observation) problems as opposed to inverse problems where the goal is to estimate the parameter θ∗ only from noisy observations of its image through an operator. The linear regression model one such inverse problem where the matrix X plays the role of a linear operator. However, in these notes, we never try to invert the operator. See [Cav11] for an interesting survey on the statistical theory of inverse problems. Sparsity adaptive thresholding estimators If we knew a priori that θ was k sparse, we could employ directly Corollary 2.8 to obtain that with probability 1 − δ, we have ls MSE(Xθ̂B ) ≤ Cδ 0 (k)  ed  σ2 k log . n 2k As we will see, the assumption ORT gives us the luxury to not know k and yet adapt to its value. Adaptation means that we can construct an estimator that does not require the knowledge of k (the smallest such that |θ∗ |0 ≤ k) and yet, ls , up to a multiplicative constant. perform as well as θˆB 0 (k) Let us begin with some heuristic considerations to gain some intuition. Assume the sub-Gaussian sequence model (2.9). If nothing is known about θ∗ it is natural to estimate it using the least squares estimator θˆls = y. In this case, σ2 d MSE(Xθ̂ls ) = |y − θ∗ |22 = |ξ|22 ≤ Cδ , n where the last inequality holds with probability at least 1 − δ. This is actually what we are looking for if k = Cd for some positive constant C ≤ 1. The problem with this approach is that it does not use the fact that k may be much smaller than d, which happens when θ∗ has many zero coordinate. If θj∗ = 0, then, yj = ξj , which is a sub-Gaussian random variable with variance proxy σ 2 /n. In particular, we know from Lemma 1.3 that with probability at least 1 − δ, r 2 log(2/δ) |ξj | ≤ σ =τ. (2.10) n The consequences of this inequality are interesting. One the one hand, if we observe |yj | ≫ τ , then it must correspond to θj∗ = 6 0. On the other hand, if |yj | ≤ τ is smaller, then, θj∗ cannot be very large. In particular, by the triangle inequality, |θj∗ | ≤ |yj | + |ξj | ≤ 2τ . Therefore, we loose at most 2τ by choosing θ̂j = 0. It leads us to consider the following estimator. Definition 2.10. The hard thresholding estimator with threshold 2τ > 0 is denoted by θˆhrd and has coordinates  yj if |yj | > 2τ , θ̂jhrd = 0 if |yj | ≤ 2τ , for j = 1, . . . , d. In short, we can write θˆjhrd = yj 1I(|yj | > 2τ ). From our above consideration, we are tempted to choose τ as in (2.10). Yet, this threshold is not large enough. Indeed, we need to choose τ such that |ξj | ≤ τ simultaneously for all j. This can be done using a maximal inequality. Namely, Theorem 1.14 ensures that with probability at least 1 − δ, r 2 log(2d/δ) max |ξj | ≤ σ 1≤j≤d n It yields the following theorem. Theorem 2.11. Consider the linear regression model (2.2) under the assumption ORT or, equivalenty, the sub-Gaussian sequence model (2.9). Then the hard thresholding estimator θˆhrd with threshold r 2 log(2d/δ) 2τ = 2σ , (2.11) n enjoys the following two properties on the same event A such that IP(A) ≥ 1−δ: (i) If |θ∗ |0 = k, k log(2d/δ) MSE(Xθˆhrd ) = |θˆhrd − θ∗ |22 . σ 2 . n (ii) if minj∈supp(θ∗ ) |θj∗ | > 3τ , then supp(θˆhrd ) = supp(θ∗ ) . Proof. Define the event A= n max |ξj | ≤ τ , j and recall that Theorem 1.14 yields IP(A) ≥ 1 − δ. On the event A, the following holds for any j = 1, . . . , d. First, observe that |yj | > 2τ ⇒ |θj∗ | ≥ |yj | − |ξj | > τ (2.12) |yj | ≤ 2τ ⇒ |θj∗ | ≤ |yj | + |ξj | ≤ 3τ (2.13) and It yields |θ̂jhrd − θj∗ | = |yj − θj∗ |1I(|yj | > 2τ ) + |θj∗ |1I(|yj | ≤ 2τ ) ≤ τ 1I(|yj | > 2τ ) + |θj∗ |1I(|yj | ≤ 2τ ) ≤ τ 1I(|θj∗ | > τ ) + |θj∗ |1I(|θj∗ | ≤ 3τ ) ≤ by (2.12) and (2.13) 4 min(|θj∗ |, τ ) It yields |θ̂hrd − θ∗ |22 = d X j=1 |θˆjhrd − θj∗ |2 ≤ 16 d X j=1 min(|θj∗ |2 , τ 2 ) ≤ 16|θ∗ |0 τ 2 . This completes the proof of (i). To prove (ii), note that if θj∗ 6= 0, then |θj∗ | > 3τ so that |yj | = |θj∗ + ξj | > 3τ − τ = 2τ . Therefore, θˆjhrd 6= 0 so that supp(θ∗ ) ⊂ supp(θˆhrd ). Next, if θˆjhrd 6= 0, then |θˆjhrd | = |yj | > 2τ . It yields |θj∗ | ≥ |yj | − τ > τ Therefore, |θj∗ | 6= 0 and supp(θˆhrd ) ⊂ supp(θ∗ ). 1 0 −2 −1 y −2 −1 y 0 1 2 Soft 2 Hard −2 −1 0 1 2 −2 −1 x 0 1 2 x Figure 2.2. Transformation applied to yj with 2τ = 1 to obtain the hard (left) and soft (right) thresholding estimators Similar results can be obtained for the soft thresholding estimator θˆsft defined by   yj − 2τ if yj > 2τ , θ̂jsft = yj + 2τ if yj < −2τ ,  0 if |yj | ≤ 2τ , In short, we can write 2.4  2τ  θ̂jsft = 1 − yj |yj | + HIGH-DIMENSIONAL LINEAR REGRESSION The BIC and Lasso estimators It can be shown (see Problem 2.5) that the hard and soft thresholding estimators are solutions of the following penalized empirical risk minimization problems: n o θ̂hrd = argmin |y − θ|22 + 4τ 2 |θ|0 θ∈IRd n o θ̂sft = argmin |y − θ|22 + 4τ |θ|1 θ∈IRd In view of (2.8), under the assumption ORT, the above variational definitions can be written as n1 o θ̂hrd = argmin |Y − Xθ|22 + 4τ 2 |θ|0 n θ∈IRd n1 o θˆsft = argmin |Y − Xθ|22 + 4τ |θ|1 n θ∈IRd When the assumption ORT is not satisfied, they no longer correspond to thresholding estimators but can still be defined as above. We change the constant in the threshold parameters for future convenience. Definition 2.12. Fix τ > 0 and assume the linear regression model (2.2). The BIC2 estimator of θ∗ in is defined by any θˆbic such that n1 o θ̂bic ∈ argmin |Y − Xθ|22 + τ 2 |θ|0 n θ∈IRd Moreover the Lasso estimator of θ∗ in is defined by any θˆL such that n1 o θ̂L ∈ argmin |Y − Xθ|22 + 2τ |θ|1 n θ∈IRd Remark 2.13. Numerical considerations. Computing the BIC estimator can be proved to be NP-hard in the worst case. In particular, no computational method is known to be significantly faster than the brute force search among all 2d sparsity patterns. Indeed, we can rewrite: n1 o n o 1 min |Y − Xθ|22 + τ 2 |θ|0 = min min |Y − Xθ|22 + τ 2 k 0≤k≤d θ : |θ|0 =k n θ ∈IRd n
To compute minθ : |θ|0 =k n1 |Y − Xθ|22 , we need to compute kd least squares estimators on a space of size k. Each costs O(k 3 ) (matrix inversion). Therefore the total cost of the brute force search is d   X d 3 C k = Cd3 2d . k k=0 Instead the the Lasso estimator is convex problem and there exists many efficient algorithms to compute it. We will not describe this optimization problem in details but only highlight a few of the best known algorithms: 1. Probably the most popular method among statisticians relies on coordinate gradient descent. It is implemented in the glmnet package in R [FHT10], 2 Note that it minimizes the Bayes Information Criterion (BIC) employed in the tradip tional literature of asymptotic statistics if τ = log(d)/n. We will use the same value below, up to multiplicative constants (it’s the price to pay to get non asymptotic results). 2. An interesting method called LARS [EHJT04] computes the entire regularization path, i.e., the solution of the convex problem for all values of τ . It relies on the fact that, as a function of τ , the solution θˆL is a piecewise linear function (with values in IRd ). Yet this method proved to be too slow for very large problems and has been replaced by glmnet which computes solutions for values of τ on a grid much faster. 3. The optimization community has made interesting contribution to this field by using proximal methods to solve this problem. It exploits the structure of the form: smooth (sum of squares) + simple (ℓ1 norm). A good entry point to this literature is perhaps the FISTA algorithm [BT09]. 4. There has been recently a lot of interest around this objective for very large d and very large n. In this case, even computing |Y − Xθ|22 may be computationally expensive and solutions based on stochastic gradient descent are flourishing. Note that by Lagrange duality computing θˆL is equivalent to solving an ℓ1 constrained least squares. Nevertheless, the radius of the ℓ1 constraint is unknown. In general it is hard to relate Lagrange multipliers to the size constraints. The name “Lasso” was given to the constrained version this estimator in the original paper of Robert Tibshirani [Tib96]. Analysis of the BIC estimator While computationally hard to implement, the BIC estimator gives us a good benchmark for sparse estimation. Its performance is similar to that of θˆhrd but without assumption ORT. Theorem 2.14. Assume that the linear model (2.2) holds where ε ∼ subGn (σ 2 ). Then, the BIC estimator θˆbic with regularization parameter τ 2 = 16 log(6) σ2 σ 2 log(ed) + 32 . n n satisfies 1 log(ed/δ) |Xθ̂bic − Xθ∗ |22 . |θ∗ |0 σ 2 n n with probability at least 1 − δ. MSE(Xθ̂bic ) = Proof. We begin as usual by noting that 1 1 |Y − Xθ̂bic |22 + τ 2 |θˆbic |0 ≤ |Y − Xθ∗ |22 + τ 2 |θ∗ |0 . n n It implies |Xθ̂bic − Xθ∗ |22 ≤ nτ 2 |θ∗ |0 + 2ε⊤ X(θˆbic − θ∗ ) − nτ 2 |θˆbic |0 . (2.14) First, note that 2ε⊤ X(θˆbic − θ∗ ) = 2ε⊤  Xθ̂bic − Xθ∗  |Xθ̂bic − Xθ∗ |2 |Xθ̂bic − Xθ∗ |2 h  Xθ̂bic − Xθ∗ i2 1 ≤ 2 ε⊤ + |Xθ̂bic − Xθ∗ |22 , 2 |Xθˆbic − Xθ∗ |2 where we use the inequality 2ab ≤ 2a2 + 12 b2 . Together with the previous display, it yields  2 |Xθ̂bic − Xθ∗ |22 ≤ 2nτ 2 |θ∗ |0 + 4 ε⊤ U (θˆbic − θ∗ ) − 2nτ 2 |θˆbic |0 (2.15) where Xθ̂bic − Xθ∗ U (θˆbic − θ∗ ) = |Xθ̂bic − Xθ∗ |2 Next, we need to “sup out” θˆbic . To that end, we decompose the sup into a max over cardinalities as follows: sup = max max sup 1≤k≤d |S |=k supp(θ)=S θ∈IRd . Applied to the above inequality, it yields  2 4 ε⊤ U (θˆbic − θ∗ ) − 2nτ 2 |θ̂bic |0   2 ≤ max max sup 4 ε⊤ U(θ − θ∗ ) − 2nτ 2 k 1≤k≤d ≤ max 1≤k≤d  |S|=k supp(θ)=S max |S|=k  2 sup 4 ε⊤ ΦS,∗ u − 2nτ 2 k , r u∈B2S,∗ where ΦS,∗ = [φ1 , . . . , φrS,∗ ] is an orthonormal basis of the set {Xj , j ∈ S ∪ supp(θ∗ )} of columns of X and rS,∗ ≤ |S| + |θ∗ |0 is the dimension of this column span. Using union bounds, we get for any t > 0,    2 IP max max sup 4 ε⊤ ΦS,∗ u − 2nτ 2 k ≥ t) 1≤k≤d ≤ |S|=k d X X k=1 |S|=k IP r u∈B2S,∗    ⊤ 2 t 1 ε ΦS,∗ u ≥ + nτ 2 k r 4 2 u∈B S,∗ sup 2 Moreover, using the ε-net argument from Theorem 1.19, we get for |S| = k, t     2 + 1 nτ 2 k  t 1 IP sup ε⊤ ΦS,∗ u ≥ + nτ 2 k ≤ 2 · 6rS,∗ exp − 4 2 2 r 4 2 8σ u∈B2S,∗   t nτ 2 k ∗ ≤ 2 exp − − + (k + |θ | ) log(6) 0 32σ 2 16σ 2   t ∗ ≤ exp − − 2k log(ed) + |θ | log(12) 0 32σ 2 where, in the last inequality, we used the definition (2.14) of τ . Putting everything together, we get   IP |Xθ̂bic − Xθ∗ |22 ≥ 2nτ 2 |θ∗ |0 + t ≤ d X X k=1 |S|=k = d   X d k=1 ≤ =  exp − d X k=1 d X k=1 k  t ∗ − 2k log(ed) + |θ | log(12) 0 32σ 2  exp −  exp −  t ∗ | | log(12) − 2k log(ed) + θ 0 32σ 2  t − k log(ed) + |θ∗ |0 log(12) 2 32σ by Lemma 2.7  (ed)−k exp −  ≤ exp − t 32σ 2  t ∗ + |θ | lo g (12) 0 32σ 2  + |θ∗ |0 log(12) . To conclude the proof, choose t = 32σ 2 |θ∗ |0 log(12)+32σ 2 log(1/δ) and observe that combined with (2.15), it yields with probability 1 − δ, |Xθ̂bic − Xθ∗ |22 ≤ 2nτ 2 |θ∗ |0 + t = 64σ 2 log(ed)|θ∗ |0 + 64 log(12)σ 2 |θ∗ |0 + 32σ 2 log(1/δ) ≤ 224|θ∗ |0 σ 2 log(ed) + 32σ 2 log(1/δ) . It follows from Theorem 2.14 that θˆbic adapts to the unknown sparsity of θ , just like θˆhrd . Moreover, this holds under no assumption on the design matrix X. ∗ Analysis of the Lasso estimator Slow rate for the Lasso estimator The properties of the BIC estimator are quite impressive. It shows that under no assumption on X, one can mimic two oracles: (i) the oracle that knows the support of θ∗ (and computes least squares on this support), up to a log(ed) term and (ii) the oracle that knows the sparsity |θ∗ |0 of θ∗ , up to a smaller logarithmic term log(ed/|θ∗ |0 ) is replaced by log(ed). Actually the latter can even be removed by using a modified BIC estimator (see Problem 2.6). The Lasso estimator is a bit more difficult because, by construction, it should more naturally adapt to the unknown ℓ1 -norm of θ∗ . This can be easily shown as in the next theorem, analogous to Theorem 2.4. Theorem 2.15. Assume that the linear model (2.2) holds where ε ∼ subGn (σ 2 ). Moreover, assume that the columns of X are normalized in such a way that √ maxj |Xj |2 ≤ n. Then, the Lasso estimator θˆL with regularization parameter r r 2 log(2d) 2 log(1/δ) 2τ = 2σ + 2σ . (2.16) n n satisfies 1 MSE(Xθ̂ ) = |XθˆL − Xθ∗ |22 ≤ 4|θ∗ |1 σ n L r 2 log(2d) + 4|θ∗ |1 σ n r 2 log(1/δ) n with probability at least 1 − δ. Moreover, there exists a numerical constant C > 0 such that r   log(2d) L ∗ IE MSE(Xθ̂ ) ≤ C|θ |1 σ . n Proof. From the definition of θˆL , it holds 1 1 |Y − XθˆL |22 + 2τ |θˆL |1 ≤ |Y − Xθ∗ |22 + 2τ |θ∗ |1 . n n Using Hölder’s inequality, it implies
|Xθ̂L − Xθ∗ |22 ≤ 2ε⊤ X(θˆL − θ∗ ) + 2nτ |θ∗ |1 − |θˆL |1 ≤ 2|X⊤ ε|∞ |θˆL |1 − 2nτ |θˆL |1 + 2|X⊤ ε|∞ |θ∗ |1 + 2nτ |θ∗ |1 = 2(|X⊤ ε|∞ − nτ )|θˆL |1 + 2(|X⊤ ε|∞ + nτ )|θ∗ |1 Observe now that for any t > 0, t2 − 2nσ2 IP(|X⊤ ε|∞ ≥ t) = IP( max |X⊤ j ε| > t) ≤ 2de 1≤j≤d Therefore, taking t = σ probability 1 − δ, p 2n log(2d) + σ p 2n log(1/δ) = nτ , we get that with |Xθ̂L − Xθ∗ |22 ≤ 4nτ |θ∗ |1 . The bound in expectation follows using the same argument as in the proof of Corollary 2.9. Notice that the regularization parameter (2.16) depends on the confidence level δ. This not the case for the BIC estimator (see (2.14)). p The rate in Theorem 2.15 if of order (log d)/n (slow rate), which is much slower than the rate of order (log d)/n (fast rate) for the BIC estimator. Hereafter, we show that fast rates can be achieved by the computationally efficient Lasso estimator but at the cost of a much stronger condition on the design matrix X. Incoherence Assumption INC(k) We say that the design matrix X has incoherence k for some integer k > 0 if X⊤ X 1 − Id ∞ ≤ n 14k where the |A|∞ denotes the largest element of A in absolute value. Equivalently, 1. For all j = 1, . . . , d, |Xj |22 1 −1 ≤ . n 14k 2. For all 1 ≤ i, j ≤ d, i 6= j, we have X⊤ i Xj ≤ 1 . 14k Note that Assumption ORT arises as the limiting case of INC(k) as k → ∞. However, while Assumption ORT requires d ≤ n, here we may have d ≫ n as illustrated in Proposition 2.16 below. To that end, we simply have to show that there exists a matrix that satisfies INC(k) even for d > n. We resort to the probabilistic method [AS08]. The idea of this method is that if we can find a probability measure that puts a positive probability of objects that satistify a certain property, then there must exist objects that satisfy said property. In our case, we consider the following probability distribution on random matrices with entries in {±1}. Let the design matrix X have entries that are i.i.d Rademacher (±1) random variables. We are going to show that most realizations of this random matrix satisfy Assumption INC(k) for large enough n. Proposition 2.16. Let X ∈ IRn×d be a random matrix with entries Xij , i = 1, . . . , n, j = 1, . . . , d that are i.i.d Rademacher (±1) random variables. Then, X has incoherence k with probability 1 − δ as soon as n ≥ 392k 2 log(1/δ) + 784k 2 log(d) . It implies that there exists matrices that satisfy Assumption INC(k) for n & k 2 log(d) , for some numerical constant C. Proof. Let εij ∈ {−1, 1} denote the Rademacher random variable that is on the ith row and jth column of X. Note first that the jth diagonal entries of X⊤ X/n is given by n 1X 2 ε =1 n i=1 i,j 6 k, the (j, k)th entry of the d × d matrix Moreover, for j = n X⊤ X n is given by n 1X 1 X (j,k) εi,j εi,k = ξ , n i=1 n i=1 i (j,k) (j,k) where for each pair, (j, k), ξi = εi,j εi,k so that the random variables ξ1 are iid Rademacher random variables. Therefore, we get that for any t > 0, IP X⊤ X − Id n ∞ n  
1 X (j,k) > t = IP max ξi >t 6 j=k n i=1 ≤ ≤ X IP 6 j=k X n   1X (j,k) >t ξi n i=1 2e− 6 j=k ≤ d2 e− nt2 2 (Union bound) (Hoeffding: Theorem 1.9) nt2 2 Taking now t = 1/(14k) yields IP X⊤ X − Id n ∞ for > n 1
≤ d2 e− 392k2 ≤ δ 14k n ≥ 392k 2 log(1/δ) + 784k 2 log(d) . For any θ ∈ IRd , S ⊂ {1, . . . , d} define θS to be the vector with coordinates  θj if j ∈ S , θS,j = 0 otherwise . In particular |θ|1 = |θS |1 + |θS c |1 . The following lemma holds Lemma 2.17. Fix a positive integer k ≤ d and assume that X satisfies assumption INC(k). Then, for any S ∈ {1, . . . , d} such that |S| ≤ k and any θ ∈ IRd that satisfies the cone condition |θS c |1 ≤ 3|θS |1 , it holds |θS |22 ≤ 2 |Xθ|22 n (j,k) , . . . , ξn (2.17) Proof. We have |XθS |22 X⊤ X |Xθ|22 1 = |XθS + XθS c |22 ≥ + 2θS⊤ θS c n n n n If follows now from the incoherence condition that |θS |21 |XθS |22 X⊤ X X⊤ X θS = |θS |22 + θS⊤ ( − Id )θS ≥ |θS |22 − = θS⊤ 14k n n n and X⊤ X 1 3 θS c ≤ |θS |1 |θS c |1 ≤ |θS |21 n 14k 14k Observe now that it follows from the Cauchy-Schwarz inequality that θS⊤ |θS |21 ≤ |S||θS |22 Thus for |S| ≤ k,  |Xθ|22 7|S|  1 ≥ 1− |θS |22 ≥ |θS |22 n 14k 2 Fast rate for the Lasso Theorem 2.18. Fix n ≥ 2. Assume that the linear model (2.2) holds where ε ∼ subGn (σ 2 ). Moreover, assume that |θ∗ |0 ≤ k and that X satisfies assumption INC(k). Then the Lasso estimator θˆL with regularization parameter defined by r r log(2d) log(1/δ) 2τ = 8σ + 8σ n n satisfies 1 ˆL log(2d/δ ) |Xθ − Xθ∗ |22 . kσ 2 n n and r log(2d/δ) |θ̂L − θ∗ |1 . kσ . n with probability at least 1 − δ. Moreover, r  L    log(2d/δ) L 2 log(2d) ∗ ˆ ˆ IE MSE(Xθ ) . kσ , and IE |θ − θ |1 . kσ . n n MSE(Xθ̂L ) = Proof. From the definition of θˆL , it holds 1 1 |Y − Xθ̂L |22 ≤ |Y − Xθ∗ |22 + 2τ |θ∗ |1 − 2τ |θˆL |1 . n n Adding τ |θˆL − θ∗ |1 on each side and multiplying by n, we get |XθˆL −Xθ∗ |22 +nτ |θˆL −θ∗ |1 ≤ 2ε⊤ X(θˆL −θ∗ )+nτ |θˆL −θ∗ |1 +2nτ |θ∗ |1 −2nτ |θˆL |1 . Applying Hölder’s inequality and using the same steps as in the proof of Theorem 2.15, we get that with probability 1 − δ, we get ε⊤ X(θˆL − θ∗ ) ≤ |ε⊤ X|∞ |θˆL − θ∗ | nτ L ≤ |θ̂ − θ∗ |1 , 2 where we used the fact that |Xj |22 ≤ n + 1/(14k) ≤ 2n. Therefore, taking S = supp(θ∗ ) to be the support of θ∗ , we get |Xθ̂L − Xθ∗ |22 + nτ |θ̂L − θ∗ |1 ≤ 2nτ |θˆL − θ∗ |1 + 2nτ |θ∗ |1 − 2nτ |θˆL |1 = 2nτ |θ̂L − θ∗ |1 + 2nτ |θ∗ |1 − 2nτ |θˆL |1 S ≤ 4nτ |θ̂SL S ∗ − θ |1 (2.18) In particular, it implies that |θ̂SLc − θS∗ c |1 ≤ 3|θ̂SL − θS∗ |1 . so that θ = θˆL − θ∗ satisfies the cone condition (2.17). Using now the CauchySchwarz inequality and Lemma 2.17 respectively, we get since |S| ≤ k, r p 2k ˆL L ∗ L ∗ ˆ |Xθ − Xθ∗ |2 . |θ̂S − θ |1 ≤ |S||θS − θ |2 ≤ n Combining this result with (2.18), we find |Xθ̂L − Xθ∗ |22 ≤ 32nkτ 2 . Moreover, it yields L ∗ r 2k ˆL |Xθ − Xθ∗ |2 n r 2k √ ≤4 32nkτ 2 ≤ 32kτ n |θ̂ − θ |1 ≤ 4 The bound in expectation follows using the same argument as in the proof of Corollary 2.9. Note that all we required for the proof was not really incoherence but the conclusion of Lemma 2.17: inf inf |S|≤k θ∈CS |Xθ|22 ≥κ n|θS |22 where κ = 1/2 and CS is the cone defined by  CS = |θS c |1 ≤ 3|θS |1 . (2.19) Condition (2.19) is sometimes called restricted eigenvalue (RE) condition. Its name comes from the following observation. Note that all k-sparse vectors θ are in a cone CS with |S| ≤ k so that the RE condition implies that the smallest eigenvalue of XS satisfies λmin (XS ) ≥ nκ for all S such that |S| ≤ k. Clearly, the RE condition is weaker than incoherence and it can actually be shown that a design matrix X of i.i.d Rademacher random variables satisfies the RE conditions as soon as n ≥ Ck log(d) with positive probability.