Lecture notes on ordinary least squares1

François Portier

November 24, 2020

1
Please let me know if you find typos or mistakes ([Link]@[Link]). The author is grateful to Joseph
Salmon ([Link] for some help on the writing of this course and for sharing some materials.
2
Contents

1 Definition of ordinary least-squares and first properties 7

1.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Existence and uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 To centre the data or not to centre the data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 The determination coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Statistical model 13
2.1 The fixed-design model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.1 Bias, variance and risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.2 Best linear unbiased estimator (BLUE) . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1.3 Noise estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 The Gaussian model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.1 The Cochran lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.2 A concentration inequality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3 The random design model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Confidence intervals and hypothesis testing 21

3.1 Confidence intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.1 Gaussian model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.2 Nongaussian case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Hypothesis testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.2 Test of no effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 Forward variable selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4 Ridge regularization 29
4.1 PCA before OLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Definition of the Ridge estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Bias and variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.4 Choice of the regularization parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5 The LASSO 33
5.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2 Theoretical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.3 Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.4 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

A Elementary results from linear algebra 39

3
B Singular value decomposition and principal component analysis 41
B.1 Matrix decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
B.2 Principal component analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

C Concentration inequalities 43

D Optimization of convex functions 45

4
Notations
• h·, ·i is the usual inner product in Rd . k · k is the Euclidean norm. The elements forming the canonical
basis of Rd are denoted by e0 , . . . , ed−1 . Additionally, the `q -norm of x ∈ Rd is denoted by kxkqq =
Pd q
k=1 xk .

• If A ∈ Rn×d is a matrix, AT ∈ Rd×n is the transpose matrix, ker(A) = {u ∈ Rd : Au = 0}.

Pd
• For any set of vectors (u1 , . . . , ud ) in Rn , span(u1 , . . . , ud ) = { k=1 αk uk : (α1 , . . . , αd ) ∈ Rd }. When
A is a matrix span(A) stands for the linear subspace generated by its columns.

• When A is a square invertible matrix, the inverse is denoted by A−1 . The Moore–Penrose inverse is
denoted by A+ . The trace of A is given by tr(A).
• The identity matrix in Rd×d is Id . The vector 1n ∈ Rn contains n ones.
Pn
• For any sequence z1 , z2 , . . ., the empirical mean over the n first elements is denoted by z n = i=1 zi /n

• When two random variables X and Y have the same distribution we write X ∼ Y .
• When Xn is a sequence of random variables that converges in distribution (resp. in probability) to X,
p
we write Xn X (resp. Xn → X).

5
6
Chapter 1

Definition of ordinary least-squares

and first properties

1.1 Definition
The general goal of regression is to predict with an output variable y ∈ R, also called the explanatory
variable, based on the observation of some input variables x = (x1 , . . . , xp )T ∈ Rp , with p ≥ 1, also called
the covariates. The statistical approach consists in learning or estimating a regression function (also called
link function) that maps any element x of the input space Rp to the output space R. Generally, the estimation
of the regression function is based on the observation of a sample made of examples. These examples are also
called the observations and each example is a pair made of an input variable together with the corresponding
output. Let n ≥ 1 denote the number of observations. Let (xi , yi )i=1,...,n be the observations such that for
each i, xi ∈ Rp and yi ∈ R. Many regression model can be estimated on the basis of these observations.
In this course, we focus on linear regression in which the regression function is defined as a simple linear
function, i.e., the variable y is modeled by θ0 + θ1 x1 + . . . + θp xp where (θ0 , . . . , θp ) are the parameters of the
linear regression function. In what follows we introduce the ordinary least squares (OLS) approach which
basically consists in minimizing the sum of squares of the distance between the observed values yi and the
predicted values at xi under the linear model.
We focus on a regression problem with n ≥ 1 observations and p ≥ 1 covariates. For notational conve-
nience, for i = 1, . . . , n, we consider yi ∈ R and zi = (xi,0 , . . . , xi,p )T = (xi,0 , xTi ) ∈ Rp+1 with xi,0 = 1.
This is only to include the intercept in the same way as the other coefficients. The OLS estimator is any
coefficient vector θ̂n = (θ̂n,0 , . . . , θ̂n,p )T ∈ Rp+1 such that
n
X
θ̂n ∈ argminθ∈Rp+1 (yi − ziT θ)2 . (1.1)
i=1

It is useful to introduce the notations

 T  T    
x1 z1 x1,0 ... x1,p y1
X =  ...  ∈ Rn×p , Z =  ...  =  ... ..  ∈ Rn×(p+1) , Y =  ...  .
      
. 
xTn znT xn,0 ... xn,p yn

The matrix X which contains the covariates is called the design matrix. The matrix Z which contains the
covariates is called the extended design matrix. With the previous notation, (1.1) becomes

θ̂n ∈ argminθ∈Rp+1 kY − Zθk2 ,

7
 140
Raw data and fitted 140
Raw data and fitted
Data Data
120 OLS-sklearn-no-intercept 120 OLS-sklearn-w-intercept
OLS-sm-w-intercept
100 100

80 80
Distance

Distance
60 60

40 40

20 20

0 0

−20 −20

0 5 10 15 20 25 30 0 5 10 15 20 25 30
Speed Speed

Figure 1.1: The dataset is the cars dataset from the R software. We use sklearn to compute OLS. The graph
on the left represents the OLS line without intercept and the graph on the right is the OLS line computed
with intercept.

where k · k stands for the Euclidean norm. As soon as θ̂n is obtained, one can use it to define the estimated
regression function
p
X
ĝ(x) = θ̂n,0 + θ̂n,k xk .
k=1

This function is often called the predictor due to its use to predict the output value based on the observation
of the input x.

1.2 Existence and uniqueness

With the above formulation, the OLS has a nice geometric interpretation : Ŷ = Z θ̂n is the closest point
to Y in the linear subspace span(Z) ⊂ Rn (where span(A) stands for the linear subspace generated by the
columns of A). Using the Hilbert projection theorem (Rn is a Hilbert space, span(Z) is a (closed) linear
subspace of Rn ), Ŷ is unique and is characterized by the fact that the vector Y − Ŷ is orthogonal to span(Z).
This property is equivalent to the so-called normal equation:

Z T (Y − Ŷ ) = 0.

Since Ŷ = Z θ̂n , we obtain that the vector θ̂n must verify

Z T Z θ̂n = Z T Y. (1.2)

Note that in contrast with Ŷ (which is always unique), the vector θ̂n is not uniquely defined without further
assumptions on the data. For instance, take u ∈ ker(Z) then θ̂n + u verifies (1.2) as well as θ̂n . The
uniqueness of the OLS is actually determined by the kernel of Z which is related to the invertibility of the
so called Gram matrix introduce below (see Exercise 1).

Definition 1. The matrix Ĝn = Z T Z/n is called the Gram matrix. Denote by Ĥn,Z ∈ Rn×n the orthogonal
projector1 on span(Z).
1 Recall that P is the orthogonal projector on E, a subspace of Rn , if and only if P 2 = P , P T = P and ker(P ) = E ⊥ .

8
 When the Gram matrix is invertible, the OLS is uniquely defined. When it is not the case, (1.1) has an
infinite number of solutions.
Proposition 1. The OLS estimator always exists and the associated prediction is given by Ŷ = Ĥn,Z Y . It
is either
(i) uniquely defined. This happens if and only if the Gram matrix is invertible, which is equivalent to
ker(Z) = ker(Z T Z) = {0}. In this case, the OLS has the following expression:
θ̂n = (Z T Z)−1 Z T Y.

(ii) or not unique, with an infinite number of solutions. This happens if and only if ker(Z) 6= {0}. In this
case, the set of solution writes θ̂n + ker(Z) where θ̂n is a particular solution.
Proof. The existence has already been shown using the Hilbert projection theorem. The linear system (1.2)
has therefore a unique solution or an infinite number of solutions depending on whether the Gram matrix is
invertible or not. Hence it remains to show that ker(Z) = ker(Z T Z) which follows easily from the identity
kZuk2 = uT Z T Zu.
When the OLS is not unique, the solution traditionally considered is
θ̂n = (Z T Z)+ Z T Y,
where (Z T Z)+ denotes the Moore–Penrose inverse of Z T Z, which always exists. For a demi-definite positive
symmetricPmatrix with eigenvectors ui and corresponding eigenvalues λi ≥ 0, the Moore–Penrose inverse is
given by i λ−1i u T
i i {λi >0} .
u 1
Corollary 1. The set of solution of OLS (1.1) is given by {(Z T Z)+ Z T Y + u : u ∈ ker(Z)}.
Proof. Let u ∈ ker(Z). Verify that (Z T Z)+ Z T Y + u is a solution (see exercise 7). Then assuming that v is
a solution, note that v − (Z T Z)+ Z T Y belongs to ker(Z).

1.3 To centre the data or not to centre the data

We now state the equivalence between this 2 procedures : doing OLS, with the intercept, on (Y, X) (as de-
scribed before) and doing OLS, without the intercept, on centered variables. The later estimation procedure
consists in the following. Define
n
X n
X
−1 −1
Y =n Yi , and X = n xi ,
i=1 i=1
T
Yc = Y − 1n Y and Xc = X − 1n X . Hence the quantities Yc and Xc are just centered version of Y and X,
respectively. Define
θ̂n,c = argminθ∈Rp kYc − Xc θk.
The associated predictor is given by
ĝc (x) = Y + θ̂n,c (x − X).
Proposition 2. It holds that
min kYc − X̃c θ̃k = min kY − Xθk.
θ̃∈Rp θ∈Rp+1

and, assuming that Z has full rank, we have the following relationship between the traditional OLS and the
OLS based on centred data,
T
(θ̂n,1 , . . . , θ̂n,p ) = θ̂n,c .
Moreover, the 2 methods give the same predictor, i.e., for all x ∈ Rp , ĝ(x) = ĝc (x).

9
Proof. See exercise 9.
We conclude with the following proposition which expresses the uniqueness condition given in Proposition
1 in terms of X.
Proposition 3. The following conditions are equivalent:
(i) ker(Z) = {0}
(ii) ker(Xc ) = {0}
(iii) The empirical covariance matrix of X, defined as XcT Xc /n, is invertible
Proof. See Exercise 5.

1.4 The determination coefficient

Pn
To avoid trivial cases, we suppose in the following that i=1 (yi − y n )2 > 0, i.e., that the sequence yi is not
constant. The determination coefficient, denoted by R2 , is defined as the quotient between the explained
sum of squares and the total sum of squares. It is given by
Pn
2 (ŷi − y n )2 kŶ − y n 1n k2
R = Pi=1 n = .
i=1 (yi − y )
n 2
kY − y n 1n k2

Because of the orthogonality between Ŷ − Y and Ŷ and between Ŷ − Y and y n 1n , we have that
Pn
2 (ŷi − yi )2 kŶ − Y k2
R = 1 − ni=1 n 2 =1− . (1.3)
kY − y n 1n k2
P
i=1 (yi − y )

The last expression involves a new quantity, called the residual sum of squares, which is small as soon as the
OLS procedure went well, i.e., as soon as the predicted values are close to the observed values. Hence the
closer to 1 the R2 the better. The following statement justifies the use of the R2 as a score supporting the
quality of the OLS estimation :
• R2 = 1 if and only if Y = Ŷ .
• R2 = 0 if and only if Ŷ = Ĥ1n Y implying that θ̂n = (y n , 0, . . . , 0) is one OLS estimator.

Exercises
Exercise 1. Show that ker(Z T Z) = ker(Z) and that span(Z T ) = span(Z T Z) (for the latter, one might first
note that ker(Z) = span(Z T )⊥ ). Deduce that the normal equations always have at least one solution.

Exercise 2. Give θ̂n ∈ R and Ŷ ∈ Rn in the case where Z = 1n and Y ∈ Rn .

Exercise 3. Show that any invertible transformation on the covariate, i.e. Z is replaced by ZA with A
invertible, does not change the predicted values Ŷ nor the predictor.
Pn
Exercise 4. Show that i=1 ˆi = 0, where ˆ = Y − Ŷ = (I − Ĥn,Z )Y .
Exercise 5. AimPis to express the uniqueness condition of the OLS in terms of the empirical covariance
n
matrix Σ̂n = n−1 i=1 (xi − X)(xi − X)T .
(a) Show that ker(Z) = ker(Z T Z).
Pn
(b) Prove that Z T Z = i=1 zi ziT .

10
(c) Verify that ker(Z) = 0 if and only if the empirical covariance matrix Σ̂n is invertible (hint : one might
work on the condition that Σ̂n is non-invertible, i.e., there exists u ∈ Rp \{0} such that Xc u = 0).

Exercise 6. Aim is to obtain the formula Ĥn,Z = Z(Z T Z)+ Z T .

(a) Verify that for any non-negative symmetric matrix A ∈ Rp×p , show that A+ A = A+ .

(b) Show that Z(Z T Z)+ Z T is idempotent and symmetric (making it an orthogonal projector).

(c) Using that Z(Z T Z)+ Z T writes as U U T for some matrix U that we shall specify, obtain that ker(Ĥn,Z ) =
ker(Z T ).

(d) Conclude showing that span(Ĥn,Z ) = span(Z).

Exercise 7. Show that θ̂n = (Z T Z)+ Z T Y is a solution of the OLS problem.

Exercise 8. Show (1.3).

Exercise 9. Aim is to prove Proposition 2.

(a) Start by obtaining that the inequality ≥ holds true.

Pn and for all w ∈ R, it holds that kW −

(b) Then show that for any collection (wi )i=1,...,n of real numbers,
w1n k ≥ kW − W 1n k, where W = (z1 , . . . , zn ) and W = n−1 i=1 zi .

(c) Find ân such that, for any θ0 ∈ R and θ̃ ∈ Rp , kY − θ0 1n − X θ̃k ≥ kY − ân (θ̃)1n − X θ̃k.

(d) Conclude that minθ∈Rp kYc − Xc θk = minθ∈Rp , θ0 ∈R kY − Z(θ0 , θ T )T k

T
(e) Use the Hilbert projection theorem to conclude that whenever ker(Z) = {0}, (θ̂n,1 , . . . , θ̂n,p ) = θ̂n,c .

Exercise 10 (on-line ols and cross-validation). The goal of this exercise is to show that the OLS estimator
θ̂n associated with design matrix X(n) ∈ Rn×(p+1) and output y(n) ∈ Rn can be easily updated when a new
pair of observation (xTn+1 , yn+1 ) ∈ R(p+1) × R is given. We apply the result to cross validation procedure in
the end.
To clarify the notation:
   
X(n) (n+1)×(p+1) y(n)
X(n+1) = ∈R , and y(n+1) = ∈ Rn+1
xTn+1 yn+1
We assume from now on that X(n) and X(n+1) are full column rank (i.e., the columns of each matrix are
independent vectors).
NB : Some of the questions require some computation (in particular obtaining (1.4) and (1.6)). Even if
you could not prove it, it can be use later.

(a) Let A, B, C, D be matrices with respective sizes (d, d), (d, k), (k, k), (k, d). Show that if A and C are
invertible, then
(A + BCD)−1 = A−1 − A−1 B(DA−1 B + C −1 )−1 DA−1 . (1.4)

(b) Obtain that

T
ζn+1 ζn+1
T
(X(n+1) X(n+1) )−1 = (X(n)
T
X(n) )−1 − (1.5)
1 + bn+1
T
where ζn+1 = (X(n) X(n) )−1 xn+1 and bn+1 = xTn+1 (X(n)
T
X(n) )−1 xn+1 .
T T
(c) Express X(n+1) y(n+1) with respect to X(n) y(n) and yn+1 xn+1 .

11
(d) Show that the OLS estimator θ̂n+1 associated with design matrix X(n+1) and output y(n+1) can be
obtained as follows:
un+1
θ̂n+1 = θ̂n + ζn+1 , (1.6)
1 + bn+1

where un+1 = yn+1 − xTn+1 θ̂n .

(e) Keeping in memory (X(n) T

X(n) )−1 and θ̂n , explain how to update θ̂n+1 using a minimal number of
operations of the kind : matrix (p + 1, p + 1) times vector (p + 1, 1). How many such operation are
needed?
(f ) Using Equation (1.5) above, show that
1
1 + bn+1 =
1 − hn+1

where hn+1 = xTn+1 (X(n+1)

T
X(n+1) )−1 xn+1 .

(g) The prediction of yn+1 given by the model is ŷn+1 := xTn+1 θ̂n+1 . With the following formula

un+1 bn+1
ŷn+1 = xTn+1 θ̂n + .
1 + bn+1

prove that
yn+1 − ŷn+1 = un+1 (1 − hn+1 ).

(h) Given some data (y, X), leave-one-out cross-validation consists in computing the risk
n
1X
Rcv = (yi − xTi θ̂(−i) )2
n i=1

where θ̂(−i) is the OLS estimator based on (y(−i) , X(−i) ), i.e., the data (y, X) without the i-th line.
Applying what have been done so far, show that
n
1X
Rcv = (yi − ŷi )2 /(1 − ĥi )2 ,
n i=1

with ĥi = xTi (X T X)−1 xi and ŷi = xTi θ̂n , θ̂n being the OLS estimator of (y, X).

12
Chapter 2

Statistical model

In the previous section, we have defined the OLS estimator based on the observed data without any assump-
tion on the generating process associated to the data. When assuming that the observations are independent
realizations of some random variables, we can rely on probability theory to further study the behaviour of
the OLS. In the following we describe different probabilistic models : fixed design model, random design
model and the Gaussian noise model.

2.1 The fixed-design model

The fixed design model takes the form:
Yi = xTi θ ? + i , for all i = 1, . . . , n,
where (xi ) is a sequence of deterministic points in Rp+1 and (i ) is a sequence of random variables in R such
that
 
1
E[] = 0, var() = σ 2 In , with  =  ...  .
 

n
For instance, (i ) can be an identically distributed and independent sequence of centred random variables
with variance σ 2 . The level of noise σ of course reflects the difficulty of the problem.
The fixed-design model is appropriate when the sequence (xi ) is chosen by the analyst, e.g., in a physics
laboratory experiment, one can fix some variables such as the temperature, or in a clinical survey one can
give to patients a determined quantity of some serum. In contrast, the random design (see Section 2.3)
model is appropriate when the covariates are unpredictable as for instance the wind speed observed in the
nature or the age of some individuals in a survey.
Based on this model, we can derive some statistical properties that we present in the following. These
properties are concerned with different types of error related to the estimation of θ ? by θ̂n and will be
obtained under the assumption that the dimension of span(X) equals p + 1, implying that ker(X) = {0} and
that θ̂n is unique. We therefore implicitly assume that n ≥ p + 1. We can now state a useful decomposition:
provided that ker(X) = {0}, it holds that
θ̂n − θ ? = (X T X)−1 X T . (2.1)

2.1.1 Bias, variance and risk

The bias, the variance and the risk are important quantities because they are measures of the estimation
quality. For instance, an estimator is accurate when the bias is 0 and the variance is small. The following
notion of bias is related to the whole statistical model (for all θ? , not for a particular one).

13
Definition 2. An estimator θ(X, Y ) is said to be unbiased if for all (X, , θ ? ) used to generate Y according
to the model, it holds that E[θ(X, Y )] = θ ? .
The risk measures the average error associated to an estimation procedure. Different notions of risk can
be defined: the quadratic risk is defined on the regression coefficients β, the prediction risk takes care of the
prediction error, i.e., the error when predicting y. Formal definitions are given below.
Definition 3. The quadratic risk associated to θ̂n estimating θ ? is
Rquad (θ̂n , θ ? ) = E[kθ̂n − θ ? k2 ].
The prediction risk is
Rpred (θ̂n , θ ? ) = E[kY ? − Ŷ k2 ]/n,
where Y ? is the prediction we would make if we knew the true regression vector, i.e., Y ? = Xθ ? .
Proposition 4. When ker(X) = {0}, the following holds:
(i) the OLS estimator is unbiased i.e., it holds that E[θ̂n ] = θ ? .

(ii) Its variance is given by var(θ̂n ) = (X T X)−1 σ 2 .

(iii) Rpred (θ̂n , θ ? ) = (p + 1)σ 2 /n.

(iv) Rquad (θ̂n , θ ? ) = tr((X T X)−1 )σ 2 .

Hence whenever the smallest eigenvalue of Ĝn is larger than b (independently of n), the quadratic risk
of the OLS decreases with the rate 1/n, which is the classical estimation rate in statistics, e.g., empirical
average estimating the expectation.

2.1.2 Best linear unbiased estimator (BLUE)

This section is dedicated to the so called Gauss-Markov theorem which asserts that the OLS is BLUE.
We introduce the following partial order (reflexivity, anti-symmetry and transitivity) on the set of sym-
metric matrices. Let V1 ∈ Rd×d and V2 ∈ Rd×d be two symmetric matrices. We write V1 ≤ V2 whenever
uT V1 u ≤ uT V2 u for every u ∈ Rd . This partial order is particularly useful to compare the covariance matrices
of estimators. Indeed if β̂1 and β̂2 are estimators with respective covariance V1 and V2 . Then, V1 ≤ V2 if
and only if any linear combination of β̂1 has a smaller variance than the same linear combination of β̂2 .
Definition 4. An estimator is said to be linear if, for any dataset (Y, X), it writes as AY , where A ∈
R(p+1)×n depends only on X.
Proposition 5 (Gauss-Markov). Under the fixed design model, among all the unbiased linear estimators
AY , θ̂n is the one with minimal variance, i.e.,
cov(θ̂n ) ≤ cov(AY ),

with equality if and only if A = (X T X)−1 X T .

Proof. First note that AY is unbiased if and only if (A − (X T X)−1 X T )Xθ ? = 0 for all θ ? , equivalently,
BX = 0 with B = (A − (X T X)−1 X T ). Consequently, using that ET = σ 2 In , cov(BY, θ̂n ) = 0. Then, just
write
cov(AY ) = cov(BY + θ̂n )
= cov(BY ) + cov(θ̂n )
= σ 2 BB T + cov(θ̂n ) ≥ cov(θ̂n ).
The previous inequality is an equality if and only if B = 0.

14
2.1.3 Noise estimation
Providing only an estimate θ̂n of θ ? is often not enough as it does not give any clue on the accuracy of the
estimation. When possible, one should also furnish an estimation of the error σ 2 . If one knew the residuals
(i ), one would take the empirical variance of 1 , . . . , n , but this is not possible. Alternatively, one can take
n
X
−1
σ̃n2 =n (Yi − Ŷi )2 .
i=1

Pn
Because of the first normal equations expressed in (1.2), we have i=1 (Yi − Ŷi ) = 0. Consequently, σ̃n2 is
simply the empirical variance estimate of the residual vector Yi − Ŷi . Noting that σ̃n2 = n−1 k(In − Ĥn,X )k2
one can compute the expectation:

E[σ̃n2 ] = σ 2 (n − p − 1)/n.

The unbiased version (which should be used in practice) is then

 
n
σ̂n2 = σ̃n2 ,
n−p−1

where from now on we assume that n > p + 1. In the case when n = p + 1 and X has rank p + 1, we obtain
that Yi = Ŷi for all i = 1, . . . , n.

2.2 The Gaussian model

Here we introduce the Gaussian model as a submodel of the fixed design model where the distribution of
the noise sequence (i ) is supposed to be Gaussian with mean 0 and variance σ 2 . The Gaussian model can
then be formulated as follows:
i.i.d.
Yi ∼ N (xTi θ ? , σ 2 ), for all i = 1, . . . , n,

where (xi ) is non-random sequence of vector in Rp+1 . We keep assuming that ker(X) = {0} in the following.

2.2.1 The Cochran lemma

ThepStudent’s t-distribution with p degrees of freedom is defined as the distribution of the random variable
X/ Z/p, where X (resp. Z) has standard normal distribution (resp. chi-square distribution with p degrees
of freedom).

Proposition 6. Under the Gaussian model, if ker(X) = {0} and n > p + 1, it holds that

• θ̂n and σ̂n2 are independent,

• n1/2 (θ̂n − θ ? ) ∼ N (0, nσ 2 (X T X)−1 ) ,

• (n − p − 1)(σ̂n2 /σ 2 ) ∼ χ2n−p−1 ,

• if ŝ2n,k is the k-th term in the diagonal of Ĝ−1

n , then

(n1/2 /ŝn,k σ̂n )(θ̂n,k − θk? ) ∼ Tn−p−1 ,

where Tn−p−1 is the Student’s t-distribution with n − p − 1 degrees of freedom.

15
Proof. For the first point, remark that X T  and (I − Ĥn,X ) are two independent Gaussian vectors:

cov(X T , (I − Ĥn,X )) = E[X T T (I − Ĥn,X )] = 0.

Then writing
(n − p − 1)σ̂ 2 = kY − Ŷ k2 = k(I − Ĥn,X )Y k2 = k(I − Ĥn,X )k2
θ̂n − θ ? = (X T X)−1 X T ,

we see that θ̂n and σ̂ 2 are measurable transformations of two independent Gaussian vector. They then
are independent. We can use for instance the following characterisation of independence, say for random
variables ξ1 and ξ2 : for any f1 and f2 positive and measurable, E[f1 (ξ1 )f2 (ξ2 )] = E[f1 (ξ1 )]E[f2 (ξ2 )].
For the second point, as  is Gaussian, one just has to compute the variance.
For the third point, let V ∈ Rn×n be an orthogonal matrix such that V = (V1 , V2 ) where V1 is a basis of
span(X), and note that V1T (I − Ĥn,X ) = 0 and V2T (I − Ĥn,X ) = V2T . As the norm is invariant by orthogonal
transformation, one has
(n − p − 1)σ̂ 2 = k(I − Ĥn,X )k2 = kV T (I − Ĥn,X )k2 = kV2T k2 .
Consequently,
n−p−1
X
(n − p − 1)(σ̂ 2 /σ 2 ) = ˜2i ,
i=1

with ˜ = V2T /σ. It remains to show that ˜ is a Gaussian vector with covariance In−p−1 .
For the fourth point, use the second point to obtain that
(n1/2 /ŝn,k σ)(θ̂n,k − θk? ) ∼ N (0, 1).

Then (n1/2 /ŝn,k σ̂n )(θ̂n,k − θk? ) writes as the quotient of two independent random variables: a Gaussian and
the square root of a chi-square. This is a Student’s t-distribution with n − p − 1 degrees of freedom.

A direct application of the previous proposition gives us the following equality, which is informative on
the estimation error, for any k = 0, . . . , p,
P(|θ̂n,k − θk? | ≥ t) = 2STn−p−1 (tn1/2 /ŝn,k σ̂n ),
where STn−p−1 is the survival function of the distribution Tn−p−1 .

2.2.2 A concentration inequality

We now provide an additional guarantee for the OLS estimator under the Gaussian model. It consists of a
concentration inequality : an upper bound on the probability that the estimation error exceeds any given
t > 0. The upper bound unsurprisingly depends on p, n, t, and the smallest eigenvalue of Ĝn .
Proposition 7. Suppose that the Gaussian model is valid. Denote by λ̂n the smallest eigenvalue of Ĝn and
suppose that λ̂n > 0 for all n ≥ 1. Then, for any k ∈ {0, . . . , p}, n ≥ 1 and δ > 0, it holds with probability
1 − δ,
s
2σ 2 ŝ2n,k log(2/δ)
θ̂n,k − θk? ≤ .
n
where ŝ2n,k = eTk Ĝ−1
n ek . Moreover,
s
2σ 2 log(2(p + 1)/δ)
max θ̂n,k − θk? ≤ .
k=0,...,p nλ̂n

16
 Pn
Proof. Let X̃i = (X T X)−1 Xi . Apply Lemma 4 of Appendix C to the sequence i=1 (uT X̃i )i to obtain that
n
! n
!
X X
T 2 2 T 2
P (X̃i u)i > t ≤ 2 exp −t /(2σ (u X̃i ) ) ,
i=1 i=1
Pn
Choosing u = ek , we have i=1 (u
T
X̃i )2 = n−1 ŝ2n,k , and using (2.1), we obtain that
 
P θ̂n,k − θk? > t ≤ 2 exp(−t2 n/(2σ 2 ŝ2n,k )).

Choose t appropriately to obtain the first inequality. The second inequality follows from ŝ2n,k ≤ λ̂−1
n =
maxkuk=1 |uT Ĝ−1
n u| and the union bound:
 
  [ n o
P max θ̂n,k − θk? > t = P |θ̂n,k − θk? | > t 
k=0,...,p
k=0,...,p
X  
≤ P θ̂n,k − θk? > t .
k=0,...,p

Remark 1. The first inequality of Proposition 7 is important as it shows that each coordinate might not
behave similarly depending
p on the associated diagonal element of Ĝn . For instance, for the intercept, the
bound just becomes 2σ 2 log(2/δ)/n. The quantity ŝn,k will play an important role in practice when building
confidence intervals (see section 3).

Remark 2. Proposition 7 suggests that the value of the smallest eigenvalue λ̂n of Ĝn plays a certain role
on the accuracy of the estimation. The smaller λ̂n the worst the estimation accuracy.

2.3 The random design model

In the random design model, we suppose that (Yi , Xi ) is a sequence of independent and identically dis-
tributed random vectors defined on the probability space (Ω, A, P). The aim is to estimate the best linear
approximation of Y1 made up with X1 in terms of L2 -risk, i.e., to find θ that minimizes E[(Y1 − X1T θ ∗ )2 ].
Such a minimizer can be characterized with the help of the normal equation. Recall that X1 ∈ Rp+1 and
X1,0 = 1 almost surely.
2
Proposition 8. Suppose that for all k = 0, . . . , p, E[X1,k ] < ∞ and E[Y12 ] < ∞, then

inf E[(Y1 − X1T θ)2 ] = E[(Y1 − X1T θ ∗ )2 ],

θ

if and only if

E[X1 X1T ]θ ∗ = E[X1 Y1 ].

Proof. Note that the minimization problem of interest is equivalent to

inf E[(Y1 − Z1 )2 ],
Z1 ∈F

where F is the linear subspace of the Hilbert space L2 (Ω, A, P) generated by X1,0 , . . . , X1,p . As F is a closed
linear subspace (because it has a finite dimension), the minimizer is unique and characterized by the normal
equations.

17
 The previous proposition does not imply that θ ∗ is unique. In fact we are facing a similar situation as
in Proposition 1 : either θ∗ is unique, which is equivalent to E[X1 X1T ] is invertible, or θ ∗ is not uniquely
defined. Note that θ ∗ is not unique whenever one variable is a combination of the others. In this case one
might consider any of the solution, e.g., θ ∗ = E[X1 X1T ]+ E[X1 Y1 ]. Some asymptotic properties are available.
They will be useful to run some statistical tests. We consider the following definition, valid for any n ≥ 1,

θ̂n = (X T X)+ X T Y.

Proposition 9. Suppose that E[X1 X1T ] and E[Y12 ] exist and that E[X1 X1T ] is invertible. Then

n1/2 (θ̂n − θ ∗ ) N (0, σ 2 G−1 ),

where σ 2 = var(Y1 − X1T θ ∗ ) and G = E[X1 X1T ]. Moreover

σ̂n2 → σ 2 , in probability.

In particular, (n1/2 /ŝn,k σ̂n )(θ̂n,k − θk? ) N (0, 1).

Proof. Note that

n1/2 (θ̂n − θ ∗ ) = n1/2 (X T X)+ X T  + n1/2 ((X T X)+ (X T X) − Ip+1 )θ ∗ .

It suffices to show that the term in the right converges to 0 in probability and that the term in the left
converges in distribution to the stated limit. The first point is a consequence of the continuity of the
determinant. The second point is a consequence of Slutsky’s theorem using the fact that the Moore-Penrose
inverse is a continuous operation. For more details, see Exercise 11.
The convergence of σ̂n2 is obtained by the decomposition

σ̂n2 = (n − p + 1)−1 k(I − Ĥn,X )k22

= (n − p + 1)−1 kk2 − T X(X T X)+ X T  .

Invoking the law of large number, we only need to show that the term on the right goes to 0 in probability.
We have

n
!T n
!
X X
T T + T −1/2 −1/2
 X(X X) X  = n Xi i Ĝ+
n n Xi i
i=1 i=1

−1
Pn
Because Ĝ+
n →G and n−1/2 i=1 Xi i N (0, G), we get that

T X(X T X)+ X T  kN (0, σ 2 Ip+1 )k2 = σ 2 χ2p+1 .

When divided by (n − p + 1) the previous term goes to 0.

Remark 3. A more general regression problem can be formulated without specifying a linear link : the
regression function f ∗ is any measurable function that minimizes the risk

R(f ) = E[(Y1 − f (X1 ))2 ].

When E[Y12 ] < ∞, the minimizer is unique and coincides, in L2 (Ω, A, P), with the conditional expectation
of Y given X1 : f ∗ (X1 ) = E[Y1 |X1 ], almost surely.

18
Exercises
Exercise 11 (Asymptotics for the OLS in Random design). Let (X1 , Y1 ), (X2 , Y2 ), . . . be an i.i.d. sequence
(1) (p)
of random vectors. Each pair (Xi , Yi ) is valued in Rp × R. Denote by XiT = (Xi , . . . , Xi ). Suppose that
(k) (l)
for all (k, l) ∈ {1, . . . , p}2 , E[|X1 X1 |] < ∞ and G = E[X1 X1T ] is invertible. The goal is to show that

n1/2 (θ̂n − θ ∗ ) N (0, σ 2 G−1 ),

where θ ∗ is defined in Proposition 8. Recall that for each n ∈ N∗+ , the OLS is given by
n
!
X
θ̂n = Ĝ+
n n−1 Xi Yi , (2.2)
i=1
Pn
with Ĝn = n−1 i=1 Xi XiT .

1. Let Sp (R) be the space of symmetric matrices with real coefficients. Let T : Sp (R) → Sp (R) be such
that T (A) = A+ . Show the continuity of T at each point A such that det(A) 6= 0. We recall that for
any A such that det(A) 6= 0, A−1 = (det(A))−1 Com(A)T where Com : Sp (R) → Sp (R) is continuous
(it is called the comatrix).

2. What is the limit of Ĝ+

n ? In which sense?
Pn
3. Show that θ̂n − θ ∗ = Ĝ+ + ∗
n µ̂n + (Ĝn Ĝn − Ip )θ with µ̂n = n
−1
i=1 Xi (Yi − XiT θ ∗ ).
√ L
4. Obtain that nĜ+
n µ̂n → N (0, Σ) where Σ needs to be determined.
√ P
5. Prove that n(Ĝ+
n Ĝn − I)β0 → 0. One can consider the event det(Gn ) 6= 0.
√ L
6. Show that n(θ̂n − θ ∗ ) → N (0, Σ).
√ ∗
θ̂n,k −θk L
7. Let k ∈ {1, . . . , p}, find ŝn,k , depending only on σ 2 and Ĝn , such that n( ŝn,k ) → N (0, 1).

8. Deduce a (1 − α)-confidence interval for θk∗ . Verify it has level 1 − α.

19
20
Chapter 3

Confidence intervals and hypothesis

testing

3.1 Confidence intervals

From a practical perspective, building confidence intervals is often an inevitable step as it permits to evaluate
the quality of the estimation. The construction of confidence intervals follows the estimation step. Intuitively,
a confidence interval is simply a region (based on the observed data) in which the parameter of interest is
most likely to lie. The accuracy/quality of the estimation is then naturally measured by the size of the
underlying confidence interval. As we shall see, the construction of a confidence interval is based on the
estimation of the variance.
We consider a regression model with n observed data points (Y, X) and we focus on the task of building
confidence intervals for the k-th coordinate θk? of the regression vector (where k =∈ {0, . . . , p}).

Definition 5. A confidence interval of level 1 − α is an interval Iˆn (Y, X) ⊂ R satisfying, for all n ≥ 1,

P(θk? ∈ Iˆn (Y, X)) ≥ 1 − α.

3.1.1 Gaussian model

Confidence intervals for the regression coefficients
Confidence intervals can be obtained easily when the assumption on the model allows to know the distribution
of the quantity θ̂n,k −θk? . This is the case for instance in the popular Gaussian model in virtue of Proposition
6. Recall that, when it exists,

ŝ2n,k = eTk Ĝ−1

n ek .

Proposition 10. In the Gaussian model, if ker(X) = {0} and n > p + 1,

     
ŝn,k σ̂n ŝn,k σ̂n
θ̂n,k + − Qn−p−1 (1 − α/2) , Qn−p−1 (1 − α/2) ,
n1/2 n1/2

where Qn−p−1 is the quantile function of the distribution Tn−p−1 , is a confidence interval of level 1 − α.

Confidence intervals for the predicted values

We are now interested in building confidence intervals for the predicted value under the true model at a
single given point x = (1, x1 , . . . , xp ) ∈ Rp . The predicted value at x under the true model is defined as

21
y ∗ = xT θ ∗ . In the Gaussian model, using preservation properties of the Student’s distribution, we find the
following confidence interval CI(x) of level 1 − α. With probability equal to 1 − α,

y ∗ ∈ CI(x),

where
q
CI(x) = xT θ̂n ± Qn−p−1 (1 − α/2)σ̂ xT (X T X)−1 x,
Pn
2
and σ̂n2 = i=1 Yi − xTi θ̂n /(n − p − 1) (it has been introduced in Chapter 2). A related question is to
build a confidence interval on the value of y (not y ∗ ) under the true model. This can be done in a similar
manner as before but one needs to pay a particular attention to the additive noise in the model. Indeed, we
have that y = y ∗ +  where  ∼ N (0, σ 2 ). It follows that

y ∈ PI(x),

with
q
PI(x) = xT θ̂n ± Qn−p−1 (1 − α/2)σ̂ 1 + xT (X T X)−1 x.

For more details on the derivation of those confidence intervals, see Exercise 12.

3.1.2 Nongaussian case

When the noise distribution is not Gaussian, the previous confidence interval has no reason to be valid. In
this case, there are basically two techniques permitting the construction of confidence intervals:
• Concentration inequalities. This usually produces pessimistic (too large) confidence interval.
• Asymptotics. This only produces asymptotically valid confidence interval (often too small).
We start by deriving 2 confidence intervals based, respectively, on two concentration inequalities : the Markov
and the Hoeffding inequalities.
Proposition 11. In the fixed design model, suppose that (for clarity) X T X = nIn and that (i ) is an
identically distributed sequence of centered random variables with variance σ 2 , then for each k ∈ {0, 1, . . . , p},
the interval
h p p i
θ̂n,k + − σ 2 /(nα) , σ 2 /(nα) ,

is a confidence interval of level 1−α. If moreover, |i | ≤ c for all i = 1, . . . , n, then for each k ∈ {0, 1, . . . , p},
the interval
h p p i
θ̂n,k + − 2c log(2/α)c/n , 2c log(2/α)/n ,

is a confidence interval of level 1 − α.

Proof. We have using (2.1), θ̂n − θ ? = X T /n. Applying the Markov inequality

P(|θ̂n,k − θk? | ≥ t) ≤ t−2 E[(θ̂n,k − θk? )2 ]

Xn
2 2
≤σ Xi,k /(t2 n2 )
i=1
= σ 2 /(t2 n),

leading to the first confidence interval.

22
 Applying Hoeffding inequality with the sequence (i ), one has
n
!
X
P(|θ̂n,k − θk? | ≥ t) ≤ 2 exp −2(nt) /2
(bi − ai )2 ,
i=1

where ai ≤ i ≤ bi . Choosing ai = −c, bi = c, we get that

P(|θ̂n,k − θk? | ≥ t) ≤ 2 exp(−t2 n/2c),

leading to the second confidence interval.

Note that the first confidence interval based on Markov inequality is very pessimistic (i.e., very large)
compared to the second one, based on Hoeffding’s inequality. This is because log(1/α) << 1/α when α → 0.
2
Proposition 12. In the random design model, suppose that E[X1,k ] < ∞ and E[Y12 ] < ∞, then
     
ŝn,k σ̂n − ŝn,k σ̂n −
θ̂n,k + − Φ (1 − α/2) , Φ (1 − α/2) ,
n1/2 n1/2

where Φ− is the quantile function of the distribution N (0, 1), is, asymptotically, a confidence interval of level
1 − α, i.e.,

lim inf P(θk? ∈ Iˆn (α)) ≥ 1 − α.

n→∞

Proof. That Xn N (0, 1) means that P (Xn ∈ [−Φ− (1 − α/2), Φ− (1 − α/2)]) → Φ(Φ− (1 − α/2)) −
−
Φ(Φ (α/2)) = 1 − α where Φ is the cumulative distribution function of N (0, 1).

3.2 Hypothesis testing

We start by recalling some definitions and some vocabulary related to statistical testing. Then we consider
no effect tests on the covariates of a regression. These tests play an important role in practice as they might
quantify the importance of each covariate in the regression. As an application, we consider the forward
variable selection method in Section 3.3.

3.2.1 Definitions
Statistical testing aims at answering whether or not an hypothesis H0 is likely. It is usually performed by
constructing a test statistic T̂n and deciding to reject, or not, whenever T̂n is in R, or not. The region R is
called the reject region. As soon as T̂n and R are specified, the process is quite simple:

Reject whenever T̂n ∈ R

Do not reject whenever T̂n ∈
/ R.

The terminology “not to reject” rather than ”to accept” comes from the fact that H0 is often too much thin
and unlikely to be “accepted”, e.g., a simple hypothesis θ1? = 3.14159. There are basically 2 kinds of error
that we wish to control:

Type-1: to reject whereas H0 is true

Type-2: not to reject whereas H0 is not true.

The proportion of Type-1 errors is called the level of the test. One minus the proportion of Type-2 errors is
called the power of the test. The consistency imposes that, for any level 1 − α, asymptotically, the level is
smaller than α while the power is one. To achieve consistency, it is natural to let the reject region depend
on α.

23
Definition 6. A statistical test (T̂n , Rα ) is said to be (asymptotically) consistent whenever for all level
1 − α ∈ (0, 1)
lim sup PH0 (T̂n ∈ Rα ) ≤ α
n→∞

lim PH1 (T̂n ∈ Rα ) = 1.

n→∞

Remark 4. In practice, a standard choice is α = 0.05. Of course when the sample size is too small one
cannot be too demanding and larger values of α might be more reasonable.

3.2.2 Test of no effect

In a linear regression model, a covariate has no effect if and only if its associated regression coefficient is null.
A test of no effect of a covariate, say the k-th, then consists in testing the nullity of its regression coefficient
θk? :
H0 : θk? = 0.

Proposition 13. Under the random design model, if E[X1 X1T ] and E[Y12 ] exist and E[X1 X1T ] is invertible,
the statistic and reject region, respectively given by
 1/2 
n
T̂n,k = |θ̂n,k |,
ŝn,k σ̂n
Rα = (Φ− (1 − α/2), +∞),
produce a consistent test.
Proof. For the level, it is very similar to confidence interval. For the power, suppose that θk? 6= 0. Let
Zn = (n1/2 /ŝn,k σ̂n )(θ̂n,k − θk? ) and q = Φ− (1 − α/2). Then T̂n,k ∈ Rα if and only if

Zn + (n1/2 /ŝn,k σ̂n )θk? < −q or Zn + (n1/2 /ŝn,k σ̂n )θk? > q.
If θk? is positive (resp. negative) one can show that the event on the right (resp. left) has probability going to
1. We consider only the case θk? > 0. It has been shown in the proof of Proposition 9 that ŝn,k σ̂n converges
in probability to a finite value. We can work on the event that ŝn,k σ̂n < M . Let K > 0. For n large enough
q − (n1/2 /ŝn,k σ̂n )θk? < −K. Hence

P (Zn + (n1/2 /ŝn,k σ̂n )θk? > q) ≥ P (Zn > −K).

Hence
lim inf P (Zn + (n1/2 /ŝn,k σ̂n )θk? > q) ≥ 1 − Φ(−K).
n→∞

But K is arbitrary and the result follows.

Remark 5. In practice, the statistic T̂n,k is scale invariant: if D is a positive diagonal matrix, then the
statistic T̂n,k constructed from the sample X is the same as the statistic T̂n,k constructed from the sample
XD.
Remark 6. In the Gaussian case, the test statistic and the reject region are given by
 1/2 
n
T̂n,k = |θ̂n,k |,
ŝn,k σ̂n
Rα = (Qn−p−1 (1 − α/2), ∞).
Such a test has a level exactly equal to 1 − α. To derive that the power goes to 1, one can assume that for
all n ≥ 1, ŝn,k σ̂n is bounded.

24
 age sex bmi bp Serum measurements output
patient x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 y
1 59 2 32.1 101 157 93 38 4 4.9 87 151
2 48 1 21.6 87 183 103 70 3 3.9 69 75
··· ··· ···
··· ··· ···
441 36 1 30.0 95 201 125 42 5 5.1 85 220
442 36 1 19.6 71 250 133 97 3 4.6 92 57

Table 3.1: The dataset is composed of n = 442 patients, p = 10 variables “baseline” body mass index,
bmi), average blood pressure (bp), etc... The output is a score corresponding to the disease evolution. Each
covariate has been standardized Efron et al. (2004).

Remark 7 (test and confidence intervals). Making no effect tests consists in rejecting whenever 0 (or more
generally any tested values) is not lying inside the confidence interval. For instance, in the random design
model, to reject is equivalent to
n1/2
|θ̂n,k | ∈ (Φ− (1 − α/2), +∞),
ŝn,k σ̂n
which is equivalent to
     
ŝn,k σ̂n − ŝn,k σ̂n −
0∈
/ θ̂n,k + − Φ (1 − α/2) , Φ (1 − α/2) .
n1/2 n1/2

3.3 Forward variable selection

The method of forward selection is a stepwise procedure that aims at selecting the most “important” vari-
ables. The method starts with no covariate and add a new one at each step. This kind of methods is
sometimes referred to as greedy methods. The criterion used to select the best covariate follows from the test
statistic for the test of no effect: n1/2 |θ̂n,k |/(ŝn,k σ̂n ). Intuitively, the larger the statistic, the more important
the effect of the k-th variable.
More formally, let X = (1n , X̃1 , . . . , X̃p ). Each (non-constant) covariate X̃k is competing against the
others via 1-dimensional regression submodels Y ' θ0 + Xk θk . For any Y ∈ Rn and X̃k ∈ Rn , define the
OLS
θ̂n (Y, X̃k ) = argmin(θ0 ,θ1 )∈R2 kY − θ0 1n − θ1 X̃k k2 .
Within each submodel, the Gram matrix and the noise level estimate are given by
Ĝn (X̃k ) = n−1 (1n , X̃k )T (1n , X̃k ),
σ̂n2 (Y, X̃k ) = (n − 2)−1 kY − (1n , X̃k )θ̂n (Y, X̃k )k2 .

Another quantity of interest is ŝn (X̃k )2 = eT1 Ĝn (X̃k )−1 e1 . The criterion used to compare the importance of
each variable is the test statistic of the test of no effect, computed within each submodel:

θ̂n (Y, X̃k )

T̂n (Y, X̃k ) = n1/2 .
ŝn (X̃k )σ̂n (Y, X̃k )
For each covariate, such a quantity is compared and the largest value is selected. This criterion has an
interpretation in terms of p-values. When the test is described by (T̂n (Y, X̃k ), Rα ), the p-value is the
smallest value of α for which we still reject. For instance, in the random design model,
inf{α ∈ [0, 1] : T̂n (Y, X̃k ) > Φ− (1 − α/2)} = 2(1 − Φ(T̂n,k )).

25
 values of the t-stat at each steps
40 step 0
step 1
step 2
30

20

10

0
0 2 4 6 8
features

Figure 3.1: The statistics of each selected variable is 0 in the next step. The intercept is the first selected
variable, then X3 , etc...

Hence taking the largest T̂n (Y, X̃k ) is equivalent to take the smallest p-value for the underlying test of no
effect. A stopping rule can be based on the p-value: stop as soon as none of the p-value is smaller than
0.05. As soon as one variable, say X̃k , is selected, one needs to account for the predictive information it
has brought in the modeling of Y . This is to prevent from selecting 2 identical covariates. This is done by
replacing the output Y by the residual Y − (1n , X̃k )θ̂n (Y, X̃k ).
Algorithm 1 (forward variable selection).
Inputs: (Y, X) a threshold pstop . Start with r = Y , S = ∅ ⊂ A = {0, . . . , p}.
(i) For each k ∈ A\S, compute T̂n (r, X̃k ).
(ii) Stop if no p-values are smaller than pstop .
Else compute k ∗ ∈ argmax T̂n (r, X̃k ).
And update S = S ∪ {k ∗ } and r = r − (1n , X̃k∗ )θ̂n (Y, X̃k∗ ).

Figure 3.3 illustrates the procedure described by Algorithm 1 applied to the “diabetes” dataset of sklearn
presented in Table 3.1.
Remark 8. Different stopping rules might be considered. For instance, in Zhang (2009), the authors rec-
ommend to consider the residuals sum of squares and to stop as soon as krk2 < .

Exercises
Exercise 12 (explicit formulas when p = 1 for prediction intervals). Let us consider the following fixed-design
one-dimensional (p = 1) linear regression model:
Yi = β0 + β1 xi + i , i ∼ N (0, σ) i.i.d. , i = 1, ..., n .
Being a particular but simply interpretable Pncase it facilitates intuitive
Pn understanding and enables easy two-
n −1 n −1
dimensional visualization.
Pn Let x = n x
i=1 i and Y = n i=1 Yi . We further assume that xi is not
constant, i.e., that i=1 (xi − xn )2 6= 0.

1. Show that the OLS estimators β̂0 and β̂1 are

Pn n
n (x − xn )(Yi − Y )
β̂0 = Y − β̂1 x n
and β̂1 = Pni
i=1
n 2
i=1 (xi − x )

26
2. Show that
1 xn2  1
eT0 (X T X)−1 e0 = + Pn n 2 and eT1 (X T X)−1 e1 = Pn n 2,
n i=1 (xi − x ) i=1 i − x )
(x

3. Give the distribution of V[β̂0 ]−1/2 (β̂0 − β0 ) and V[β̂1 ]−1/2 (β̂1 − β1 )
1 xn2  σ̂ 2
V[β̂0 ] = σ̂ 2 + Pn n 2 and V[β̂1 ] = Pn n2 2 ,
n i=1 (xi − x ) i=1 (xi − x )

1
Pn
2
where σ̂ 2 = n−2 i=1 Yi − (β̂0 + β̂1 xi ) .
4. Give the reject region for the test H0 : βj = 0.

5. For a new pair (Y, x) observed from the Gaussian model above, the value β̂0 + β̂1 x is called the point
prediction. Show that

(β̂0 + β̂1 x) − (β0 + β1 x) Y − (β̂0 + β̂1 x)

q n )2
∼ t(n − 2) and q n 2
∼ t(n − 2) .
σ̂ n1 + Pn(x−x (xi −xn )2 σ̂ 1+ 1
+ Pn(x−x ) n 2
i=1 n i=1 (xi −x )

6. Build confidence intervals for (β0 + β1 x) and Y . Note that these intervals correspond, respectively, to
CI and PI given in section 3.1.1. The last one is often called prediction interval.

27
28
Chapter 4

Ridge regularization

In this chapter we use the singular-value decomposition (SVD), a matrix decomposition presented in Ap-
pendix B. The SVD of X shall provide useful expression for quantities related to the OLS estimate. The
SVD is also important to understand principal component analysis (PCA), a method that compresses the
data without loosing too much information, also presented in Appendix B.
The ridge estimator is introduced to overcome the issues caused by poorly conditioned Gram matrix
Ĝn , i.e., P
when some of the eigenvalues are too small. As indicated by the singular-value decomposition
r
of X = k=1 si ui viT , where r stands for the rank of X and si (resp. ui and vi ) are the singular-values
Pr
(resp. singular-vector) of X, we have that θ̂n = k=1 si vi uTi y. Consequently, the estimate is numerically
unstable as soon as some of the si are close to 0. As we can see looking at the variance of the OLS or at
Proposition 7, the smallest eigenvalues of Ĝn have a bad influence on the statistical behaviour of the OLS.
The ridge estimator is a solution to control these bad effects due to poor conditioning. Before going through
the definition of the ridge estimator, we discuss another method which consists in doing PCA before running
OLS.
For ease of notation (to avoid working with X̃c as before), in the rest of the chapter, we consider the
formulation of OLS without intercept with centered variable Y ∈ Rn and X ∈ Rn×p . As established in
Proposition 2, this is equivalent to include an intercept in the OLS with non-centered variables.

4.1 PCA before OLS

Many practitioners are familiar with the method of combining PCA and OLS. In addition to visualize and
explore the centered covariates X, aim is to reduce the number of covariates to avoid inverting a possibly too
large matrix X T X. Be careful that after running PCA, due to its definition (see Definition 8 in Appendix
B), the prediction must be operated with respect to the centered covariate X and centered output Y . Hence
the intercept is no longer necessary as explained in Proposition 2. The algorithm is as follows.
Algorithm 2 (PCA before OLS).
Inputs: (Y, X), centered variables, and an integer k (the number of components to keep). output: prediction
of Y at x ∈ Rp .
(i) Do a PCA on X and keep the k first components U1 , . . . , Uk .
Pk
(ii) Let Pk = i=1 Ui UiT . Compute θ̂n,k , the OLS associated with (Y, XPk ).

(iii) Return the prediction xT Pk θ̂n,k .

Trying to legitimize the approach, one can write

kY − X θ̂n k ≤ kY − XPk θ̂n,k k ≤ kY − X θ̂n k + kX(θ̂n − θ̂n,k )k + k(XPk − X)θ̂n,k k

29
in which, by Proposition 26, the last term should be small. However, the second term in the right hand
side might be large. The lack of guarantee for this approach is due to the fact that the PCA used from the
beginning is independent of the output Y . Doing such a process might result in some loss in accuracy.

4.2 Definition of the Ridge estimator

For centered variables (Y, X), the ridge estimator is defined as a solution of the following minimization
problem

kY − Xθk2 + nλkθk2 , (4.1)

where λ > 0, called the regularization parameter, is fixed by the analyst. Before dealing with the choice of
λ, we describe some properties of the ridge estimate. First of all, let us briefly state some simple remarks:

• Intuitively, when λ → 0, we obtain the OLS. When λ → +∞, we estimate 0.

• Doing ridge is adding a regularization term to the square

Pploss of OLS, aiming to penalize for large
coefficients in θ. Other norms might be used such as k=1 |θk | (see the next chapter about the
LASSO).

• As the expression in (4.1) is a Lagrangian with constraint kθk2 ≤ c the Ridge is an OLS under
constraints. The link between c and λ is not explicit.

• To make the ridge estimate scale invariant, one might replace X by XD−1/2 where D is the diagonal
matrix with entries eTk X T Xek . Actually, this normalization permits to justify having 1 single parameter
λ to control the influence of the penalty. The ridge estimate is classically defined without intercept (to
prevent from penalizing the intercept). Hence one needs to first center Y and X so that the intercept
of the OLS is automatically 0.

Proposition 14. The minimizer of (4.1) exists and is unique. It is given by

θ̂n(rdg) = (X T X + nλIp )−1 X T Y.

Proof. Let f denote the objective function of (4.1). Considering the behaviour of f at the limit of the
domain, there exists A such that whenever kθk > A, f (θ) > f (0). But the set kθk ≤ A is compact and so a
minimum exists and is achieved. Note that for any θ,

f (θ) − f (0) = −2 < Y, Xθ > +kXθk2 + nλkθk2

= −2 < Y, Xθ > +kAθk2 ,

with A = ((X T X) + nλIp )1/2 a positive matrix. For uniqueness, note that, for any u and v, we have

ktu + (1 − t)vk2 = tkuk2 + (1 − t)kvk2 − t(1 − t)ku − vk2 . (4.2)

Then suppose that θ1 and θ2 are two distinct minimizers with f ∗ = f (θ1 ) = f (θ2 ). We have, from (4.2),

f (tθ1 + (1 − t)θ2 )
= tf (θ1 ) + (1 − t)f (θ2 ) − t(1 − t)kA(θ1 − θ2 )k2 < f ∗ .
(rdg)
Hence θ̂n is unique. The first order equation is

((X T X) + nλIp )θ = X T Y.

30
4.3 Bias and variance
We have seen that, similarly to the OLS, the ridge estimator is the solution of a linear system of equations.
In the ridge system of equations the matrix that was previously X T X in the OLS is now replaced by
X T X + nλIp . As λ is chosen by the user, it allows us to control the smallest eigenvalue of the underlying
Gram matrix. Such a change of course influence the bias and the variance of the estimate. To express these
quantities, we consider the fixed design model.
Proposition 15. In the fixed-design model :
(rdg)
(i) The bias of the ridge is E[θ̂n ] − θ ? = −λn(X T X + nλIp )−1 θ ?
(rdg)
(ii) The variance of the ridge estimator expresses as var(θ̂n ) = σ 2 (X T X + nλIp )−1 X T X(X T X +
nλIp )−1 .
(rdg)
(iii) We have that var(θ̂n ) < var(θ̂n ), where θ̂n is the OLS solution.
Proof. For the last point, we use the SVD of X to write that
p
X s2k
var(θ̂n(rdg) ) = σ 2 ui uTi .
(s2k + nλ)2
k=1

In terms of eigenvalues λ̂k associated to Ĝn , we have

p
X λk
var(θ̂n(rdg) ) = σ 2 uk uTk .
(λk + nλ)2
k=1

Doing the same for θ̂n and using that λk /(λk + λ)2 < 1/λk , we obtain the result.

4.4 Choice of the regularization parameter

As we have seen before, ridge regression reduces the variance of the OLS but introduces some bias. Actually
this is the parameter λ that decides whether we reduce the bias, λ → 0, or the variance, λ → ∞. As it cannot
be accomplished simultaneously, we are facing a trade-off commonly known under the name of bias-variance
trade-off. In the next few lines, we promote the use of cross validation to select the parameter λ. This
technique of cross validation works in more general context and is of common use as soon as one needs to
choose a parameter to run a method. Examples include the choice of the bandwidth in kernel smoothing
methods, the choice of the scale parameter in RKHS as well as the choice of the cut-off parameter in Huber
regression.
Divide the data (Y, X) according to the lines into K-folds of (approximately) equal size bK/nc. Let
(Y(k) , X(k) ) (resp. (Y−(k) , X−(k) ) denote the observation in the k-th fold (resp. all the observation outside
the k-th fold). Proceed as follows:
(rdg)
(i) Compute θ̂n,k based on each sample (Y−(k) , X−(k) ).
(rdg)
(ii) Compute the (unnormalized) prediction error over each fold Y(k) − X(k) θ̂n,k . The risk is given by
K
(rdg)
X
R̂(λ) = kY(k) − X(k) θ̂n,k k2 .
k=1

The quantity R̂(λ) reflects the prediction risk associated to λ. It is then natural to minimize R̂ over
λ ∈ (0, ∞). In practice, this is usually done by taking a finite grid.
Remark 9. A computational advantage of using the SVD is that even if considering many values of λ the
SVD could be done once for each fold.

31
Exercises
Pr
Exercise 13. Recall the SVD of X = V SU T = k=1 si vk uTk where r = rank(X).
Pr Pr
1. Show that θ̂n = k=1 s−1 T +
k uk vk y = X Y and its variance is var(θ̂n ) = σ
2 −2 T
k=1 si uk uk .

2. Show that

(X T X + nλIp )−1 X T = X T (XX T + nλIn )−1

(hint : one might prefer to use the complete SVD rather than its reduced form)
3. If n  p, give an efficient method that would compute the Ridge estimator and would cost less than
the formula of Proposition 14. Compare the number of operations required.

32
Chapter 5

The LASSO

The LASSO (least absolute shrinkage and selection operator), introduced in Tibshirani (1996) is a regression
technique that consists in minimizing the usual least-squares loss with an `1 -norm regularization. As another
regularization method, it is similar to the Ridge method, presented in the previous chapter, which uses the
`2 -norm to regularize. In contrast with the Ridge, some of the LASSO coefficients are usually equal to 0,
meaning that the corresponding variables are no longer included in the predictive model. The LASSO thus
achieves in the mean time estimation and variable selection.

5.1 Definition
As for the Ridge estimator, the LASSO is usually defined with centered variables so that we can skip
estimating the intercept. We consider the following framework: X ∈ Rn×p denote the covariates vector is
such that 1Tn X = 0 and Y ∈ Rn is the output and satisfies 1Tn Y = 0. In other words, X and Y are supposed
to have (empirical) mean 0. The LASSO estimate is defined by
 
1 2
θ̂LASSO ∈ argminθ∈Rp kY − Xθk2 + λkθk1 , (5.1)
2
where k · kq stands for the `q -norm. Some remarks are of order:
• Note that the solution of (5.1) can be recovered by working with non centered variables X and Y and
adding an intercept.
• As for the Ridge, it is standard to let the LASSO estimate be scale invariant. This is usually done
by an additional standardization step that consists in using XD−1/2 in place of X, where D is the
diagonal matrix with entries eTk X T Xek , k = 1, . . . , p (ek being the k-th element of the canonical basis
of the space Rp ).
• The LASSO is not unique. Some conditions for uniqueness are given and discussed in Tibshirani (2013).
• In contrast with the Ridge approach, the `1 -penalty of the LASSO objective function allows to shrink
to 0 the coefficients in θ̂LASSO associated to the variables that are useless to predict Y .

5.2 Theoretical properties

From a theoretical perspective, the LASSO takes advantage of sparse regression models. A regression model
is sparse whenever many of the coefficients of the parameter vector θ are equal to zero, i.e., many of the
covariates are useless to predict Y . We consider the Gaussian regression model
Y = Xθ ∗ + , with  ∼ N (0, σ 2 In ), (5.2)

33
Figure 5.1: Graphical representation of the Ridge and LASSO penalties with the level set of the quadratic
loss.

and we define the active set S ∗ ⊂ {1, . . . , p} as

S ∗ = {j = 1, . . . , p : θj∗ 6= 0}.

The number of elements in S ∗ , that we denote by s, quantifies the level of sparsity associated to the regression
model. We will see that the generalization bounds for the LASSO improve whenever s becomes small. We
follow the approach presented in Hastie et al. (2015), in which the theoretical analysis of the LASSO is carried
out using the restricted eigenvalue condition. This condition is basically dealing with the smallest eigenvalue
of the matrix X T X. It is called “restricted” because it is concerned only with particular eigenvectors that
are “away” from the not active directions. More precisely, it only considers vectors living in certain cone
that leaves away the direction S ∗c . The collection of cones of interest are now defined. For α > 0 and
S ⊂ {1, . . . , p}, we set

C(α, S) = {u ∈ Rp : kuS c k1 ≤ αkuS k1 }.

The restricted eigenvalue condition (RE) for (γ, α, S) is satisfied whenever

n−1 kXuk22 ≥ γkuk22 , ∀u ∈ C(α, S). (5.3)

The following lemma is crucial to understand the role played by the cone C(3, S ∗ ) in the analysis of the
LASSO.
Lemma 1. Whenever λ ≥ 2kX T k∞ , then

0 ≤ kX ûk22 ≤ λ(3kûk1 − kûS ∗c k1 ). (5.4)

In particular, (θ̂LASSO − θ ∗ ) ∈ C(3, S ∗ ).

Proof. Define

G(u) = kY − X(θ ∗ + u)k2 /2 + λkθ ∗ + uk1 = k − Xuk2 /2 + λkθ ∗ + uk1 .

Let û = θ̂LASSO − θ ∗ . Because G(θ̂LASSO ) ≤ G(0), we have

kX ûk22 /2 ≤ h, X ûi + λ(kθ ∗ k1 − kθ ∗ + ûk1 ).

From the triangle inequality, k(θ ∗ − (−û))S ∗ k1 ≥ |kθS∗ ∗ k1 − kûS ∗ k1 | ≥ kθS∗ ∗ k1 − kûS ∗ k1 , implying that

kθ ∗ k1 − kθ ∗ + ûk1 = kθ ∗ k1 − k(θ ∗ + û)S ∗ k1 − k(θ ∗ + û)S ∗c k1

≤ kθ ∗ k1 − kθS∗ ∗ k1 + kûS ∗ k1 − k(θ ∗ + û)S ∗c k1
= kûS ∗ k1 − kûS ∗c k1 .

34
From Holder inequality, we get h, X ûi ≤ kX T k∞ kûk1 , which leads to
kX ûk22 /2 ≤ kX T k∞ kûk1 + λ(kûS ∗ k1 − kûS ∗c k1 ).

Consequently, because 2kX T k∞ ≤ λ, we obtain that

0 ≤ kX ûk22 /2 ≤ λ(kûk1 /2 + kûS ∗ k1 − kûS ∗c k1 ),
and the conclusion follows.

Now we can state the main result dealing with the analysis of the LASSO error.
Theorem 1. Under the Gaussian model (5.2), p assume that ∀k ∈ {1, . . . , p}, (X T X)k,k ≤ n and that RE for
∗
(γ, 3, S ) is satisfied. Then provided that λ = 2 2nσ 2 log(2p), we have with probability 1 − δ that
64sσ 2 log(2p)
kX(θ̂LASSO − θ ∗ )k22 ≤ .
γ
In addition, we have with probability 1 − δ that
r
∗ 6 2σ 2 s log(2p)
kθ̂LASSO − θ k2 ≤ .
γ n

Proof. Let û = θ̂LASSO −θ ∗ . Suppose for now that λ ≥ 2kX T k∞ (this will be shown to hold with probability
1 − δ at the end of the proof). We have, from (5.4) and Jensen inequality, that
√ √
kX ûk22 ≤ 3λkûS ∗ k1 ≤ 3λ skûS ∗ k2 ≤ 3λ skûk2 . (5.5)
By Lemma 1 and the RE condition, it holds that
kûk22 ≤ (γn)−1 kX ûk22 .
Injecting this in (5.5), we obtain the first and second statement. It remains to show that with probability
1 − δ, λ ≥ 2kX T k∞ . This follows from the use of Lemma 4, a Gaussian concentration result stated in
Appendix C. Let k ∈ {1, . . . , p}, note that (X T )k is distributed as N (0, (X T X)k,k σ 2 ). Applying Lemma 4
gives that
P(|(X T )k | > t) ≤ 2 exp(−t2 /(2(X T X)k,k σ 2 )).
From the union bound, it follows that
P(kX T k∞ > t) ≤ (2p) exp(−t2 /(2nσ 2 )),
or equivalently, that with probability 1 − δ,
p
kX T k∞ ≤ 2nσ 2 log(2p/δ) = λ/2.

Some other (asymptotic) properties of the LASSO are derived in Knight and Fu (2000). The authors
assume the following
n−1 X T X → C, a positive definite matrix, (5.6)
and
(Yi − θ0? − XiT θ ? )i is an iid sequence with mean 0 and variance σ 2 . (5.7)

Theorem 2 (Knight and Fu (2000)). Suppose that (5.6) and (5.7) hold and that λ/n → 0, then θ̂LASSO →
√ √
θ ? , in probability. If moreover λ/ n → λ0 ≥ 0, then n(θ̂LASSO − θ ? ) converges weakly.

35
5.3 Computation
In contrast with the OLS or the Ridge, we have no closed formula for the LASSO solutions. This is due
to the lack of smoothness of the `1 -norm. In particular, the traditional first order conditions are derived
using subgradients (rather than gradients). In Appendix D, some basic definitions and properties are given
concerning subgradients and subdifferentials.
The following remark is helpful to characterize the set of LASSO solutions : if f : Rp → R is a convex
function, the point x? is a minimum if and only if 0 ∈ ∂f (x? ), where ∂f is the subdifferential of f (see
Appendix D for details). This is often refereed to as the Fermat’s rule.

Proposition 16. Denote by Xk the k-th column of X. The LASSO solution satisfies

{sign(θ̂LASSO,k )} if θ̂LASSO,k 6= 0
∀k = {1, . . . , p}, hXk , Y − X θ̂LASSO i ∈
[−1, 1] if θ̂LASSO,k = 0

Actually, the previous set of equations has no explicit solutions. The LASSO problem becomes much
simpler when we fix all coordinates except one, and try to minimize with respect to this coordinate. For this
reason, the LASSO is usually computed using a coordinate descent, i.e., by iteratively solving the first order
conditions (involving subgradients) for each coordinate. For any λ ≥ 0, define the function ηλ : R → R given
by

 z+λ if z < −λ
ηλ (z) = 0 if z ∈ [−λ, λ] (5.8)
z−λ if z > λ


The function ηλ intervenes in solving least-squares with an absolute penalty (see Exercise 14). Note that
ηλ (z) = sign(z)(|z| − λ)+ . As shown in the following development, the function η is useful to update
the LASSO in the coordinate descent algorithm. To avoid trivial cases, we suppose in the following that
kXk k2 6= 0. In practice, one just need to remove the constant variables. Let
X
zk = Yk − Xj θj .
j6=k

Minimizing (5.1) with respect to the k-th coordinates is the same as minimizing
 
1 2 2
(kzk − Xk θk k2 − kzk k2 ) + λ|θk |
2
 
1 2 2
= (−2hzk , Xk iθk + θk kXk k2 ) + λ|θk | ,
2

which is the same as minimizing

(   2 )
1 Xk λ
zk , − θk + |θk | .
2 kXk k22 kXk k22

Consequently, the update is given by

 
Xk
θ̂k = ηλ/kXk k22 zk , .
kXk k22

As it is standard, one can start with a Ridge solution and then update each coordinates with the previous
formula.

36
5.4 Extensions
Among the extension of the LASSO, we have the LSLASSO (Least-Square LASSO) which consists in (i)
running the LASSO to find the support and (ii) applying OLS on the non-zero coefficients. Another extension
of LASSO is Elastic Net, introduced in Zou and Hastie (2005), which computes the regression coefficient
 
1 1
θ̂E-NET ∈ argminθ∈Rp kY − Xθk22 + λ{αkθk1 + (1 − α)kθk22 } . (5.9)
2 2

The previous estimate is unique. Finally, a notable extension is adaptive LASSO, introduced in Zou (2006).
The adaptive LASSO is an iterative strategy that attributes weights to each coefficient |θk | in the penalty
1/2
term. The weights might take the form 1/θ̂OLS,k penalizing mostly the small OLS coefficients. The resulting
estimate recover the support with probability going to 1.

Exercises
Exercise 14. Show that ηλ (z) = argminx∈R {(z − x)2 /2 + λ|x|} where ηλ is defined in (5.8).

Exercise 15. Following the approach given for the LASSO, find the update for the Elastic Net defined in
(5.9).

37
38
Appendix A

Elementary results from linear algebra

The vector space Rd is endowed with the usual inner product

d
X
∀(u, v) ∈ Rd × Rd , hu, vi = uT v = u k vk ,
k=1

where uT stands for the transpose of u. If hu, vi = 0 we say that u and v are orthogonal and we write u ⊥ v.
If E is a set of vector in Rd , we define its orthogonal complement as

E ⊥ = {u ∈ Rd : xT u = 0, ∀x ∈ E}.

Proposition 17. If E is a linear subspace of Rd , then (E ⊥ )⊥ = E.

For any matrix A ∈ Rp×d , define

span(A) = {Ax : x ∈ Rd },
ker(A) = {x ∈ Rd : Ax = 0}.

The set span(A) is called the image of the matrix A. It is the linear space generated by the columns of A.
The set ker(A) is called the kernel of A. Both sets are linked by the following property.

Proposition 18. Let A ∈ Rp×d . Then ker(A) = span(AT )⊥ .

Proposition 19. Let A ∈ Rp×d . Then ker(A) = {0} if and only if span(AT ) = Rd . Consequently, if p < d
then ker(A) 6= {0}.

Let A ∈ Rp×d , b ∈ Rd . Let S be the set of solutions of the linear system Ax = b.

Proposition 20. We have only three possible configurations:

1. S contains only one element,

2. S = ∅,

3. the number of elements in S is infinite.

Note that S is empty if and only if b ∈

/ span(A).

Proposition 21. Suppose that b ∈ span(A) and let x0 ∈ S, then

S = x0 + ker(A).

39
 We now recall a classical result called the spectral decomposition of symmetric matrices or the eigen
decomposition of symmetric matrices.
Proposition 22. Let A ∈ Rd×d be a symmetric matrix. Then there exist λ1 ≥ . . . ≥ λd , called eigenvalues,
and an orthonormal matrix U ∈ Rd×d (i.e., U T U = Id ) of eigenvectors, such that A = U DU T , where
D = diag(λ1 , . . . , λd ).

Definition 7. A linear transformation P ∈ Rd×d is called orthogonal projector if P 2 = P and P T = P .

The next proposition says that an orthogonal projector is characterized by its span and, therefore, by its
kernel from Proposition 17 and 18.
Proposition 23. The eigenvalues of an orthogonal projector are either 1 or 0. Hence any orthogonal
projector can be written as U U T where U ∈ Rp×r forms a basis of span(P ).
Proposition 24. The trace of an orthogonal projector is equal to the dimension of its span.

40
Appendix B

Singular value decomposition and

principal component analysis

Before we present the method of principal component analysis (PCA), it is appropriate to recall some matrix
decomposition results and more particularly the singular value decomposition (SVD).

B.1 Matrix decomposition

The usual eigen decomposition of symmetric matrices can be extended to arbitrary matrices (even not
squared matrix). The price to pay is that the left and right eigenvectors are different. This is called the
SVD.
Proposition 25. Let X ∈ Rn×p . Then there exist two orthogonal matrices : U ∈ Rp×p and V ∈ Rn×n of
singular vectors; and s1 ≥ . . . ≥ smin(n,p) ≥ 0, called singular values, such that

X = V SU T ,

where S ∈ Rn×p contains 0 everywhere except on the diagonal formed by (s1 , . . . , smin(n,p) ) .
Proof. Without loss of generality, we suppose that p ≤ n. Otherwise we apply the result to the X T .
Applying Proposition 22 to X T X, there exists U ∈ Rp×p such that U T (X T X)U is diagonal with r positive
coefficients. Hence U1T (X T X)U1 = D ∈ Rr×r and XU2 = 0. Take V1T = D−1/2 U1T X T (an orthogonal set of
r vectors : V1T V1 = Ir ) to find that V1T XU1 = D1/2 . Consequently, V1T X(U1 , U2 ) = (D1/2 , 0). Remarking
that v orthogonal to V1 means that v T XU1 = 0 implying that v T X(U1 , U2 ) = 0 leading to v T X = 0.
Now taking V2 such that V = (V1 , V2 ) ∈ Rn×p is orthogonal, we obtain the claimed decomposition with
S 2 = diag(d1 , . . . , dp ).
We have the following reduced SVD formula, if r ≥ 1 stands for the dimension of span(X),

X = Ṽr S̃r ŨrT ,

where Ũr = (U1 , . . . , Ur ), Ṽr = (V1 , . . . , Vr ), and S̃r ∈ Rr×r contains only the positive singular-values.
An attractive property of the SVD is that it defines subspaces on which one can project the data X
without loosing too much.
Proposition 26. Let X ∈ Rn×p . For any projector P ∈ Rp×p with rank smaller than k, it holds that

kX − XPk kF ≤ kX − XP kF ,

Ui UiT .
P
where Pk = i≤k

41
Proof. Suppose that 1 ≤ k < r. By Pythagorean identity, kX − XP k2F = kXk2F − kXP k2F . Hence one just
has to show that kXPk k2F ≥ kXP k2F . Considering the reduced SVD X = Ur Sr VrT , we have

kXP k2F = tr (P Ur )S 2 (P Ur )T


 
X
= tr  s2i Wi Wi 
i≤r
X
= s2i kWi k22 ,
i≤r

with Wi = P Ui and the constraints that kWi k22 ≤ 1 and i≤r kWi k22 ≤ k. Note that this corresponds to the
P
optimization problem
X X
max s2i mi u.c. mi ∈ (0, ki ), mi ≤ k,
m1 ,...,mr0
i≤r 0 i≤r 0

0
r 0 with r ≤ r and ki ≥ 1 stands for the multiplicity. We derive the
in which we suppose that s1 < . . . < sP
maximum. Note first that necessarily i≤r0 mi = k. Then if i is the first index such that 0 < mi < ki , the
function cannot achieve its P
maximum. Then we get that the maximizer is achieved when mi is either 0 or
1. Clearly the maximum is i≤k s2i which is achieved when P = i≤k Ui UiT .
P

B.2 Principal component analysis

nT
Definition 8. Let X ∈ Rn×p and define Xc = X − 1n X . The PCA of X of degree k is given by the k first
elements of the SVD of Xc , i.e., the singular values (s1 , . . . , sk ), the principal components U1 , . . . Uk and the
principal axes V1 , . . . , Vk .
Introduce the estimated covariance matrix

Σ̂n = n−1 XcT Xc .

Proposition 27. The principal components U = U1 , . . . Uk forms a set of orthonormal vectors along which
the empirical variance is maximal, i.e.,
X X
UiT Σ̂n Ui ≥ ŨiT Σ̂n Ũi ,
i≤k i≤k

for any (Ũ1 , . . . , Ũk ) orthonormal vectors. The principal components U can be obtained by an eigendecom-
position of Σ̂n .
Proof. Take Ũ and U as define in the statement. Define P̃ = Ũ Ũ T and P = U U T , the associated projectors
of rank k. Write
X
UiT Σ̂n Ui = tr(Σ̂n P ) = n−1 tr(XcT Xc P ) = n−1 kXc P k2F .
i≤k

Using Proposition 26 and the Pythagorean identity, we get that kXc P k2 ≥ kXc P̃ k2F .

Remark 10. As the PCA of X depends on the scale of each covariate, one may prefer in practice to rescale
the matrix X before running the PCA algorithm. This can be done by taking XD−1/2 rather than X, with
D equal to the diagonal matrix whose elements are eTk Σ̂n ek , k = 1, . . . , n. Then each covariate of XD has
the same empirical variance.

42
Appendix C

Concentration inequalities

To derive concentration inequalities for the errors of the estimators, we use the notion of sub-Gaussianity as
defined for instance in (Boucheron et al., 2013, Section 2.3). Recall that the moment generating function of a
Gaussian random variable W with mean µ and variance σ 2 is equal to λ 7→ E[exp(λW )] = exp(µλ + λ2 σ 2 /2).
Definition 9. A centered random variable Y is sub-Gaussian with variance factor τ 2 > 0, notation Y ∈
G(τ 2 ), if log E[exp(λY )] ≤ λ2 τ 2 /2 for all λ ∈ R.
If Y ∈ G(τ 2 ), then necessarily var(Y ) ≤ τ 2 (Boucheron et al., 2013, Exercise 2.16). Centered, bounded
random variables taking values in an interval [a, b] are sub-Gaussian with variance factor at most (b − a)2 /4
(Boucheron et al., 2013, Lemma 2.2). Chernoff’s inequality provides exponential bounds on the tails of
sub-Gaussian random variables.
Lemma 2 (Chernoff). If Y ∈ G(τ 2 ), then P (Y > t) ≤ exp(−t2 /(2τ 2 )).

Proof. For any λ ∈ R, we have 1Y >t ≤ exp(λ(Y − t)). Hence it holds that P (Y > t) ≤ E[exp(λ(Y − t))] ≤
exp(λ2 τ 2 /2 − λt). Minimizing the previous bound in λ gives λ = t/τ 2 , from which we deduce that the stated
bound.
Finally, the sum of independent sub-Gaussian variables is again sub-Gaussian. This is the statement of
the following Lemma, which proof is left as an exercise.
2
Lemma i≥1 is a sequence of independent random variables such that, for all i ≥ 1, Yi ∈ G(τi ),
Pn 3. If (Yi )P n
then i=1 Yi ∈ G( i=1 τi ).
To conclude we state a concentration inequality for the sum of independent sub-Gaussian random vari-
ables. This is just a consequence of the two previous Lemma.

Lemma 4. If (Yi )i≥1 is a sequence of independent random variables such that, for all i ≥ 1, Yi ∈ G(τi2 ),
then, for all n ≥ 1, and t ≥ 0,
n
! n
!
X X
2 2
P Yi > t ≤ 2 exp −t /(2 τi )
i=1 i=1
.

43
44
Appendix D

Optimization of convex functions

We recall here some definitions and basic properties dealing with the minimization of convex functions.
Definition 10 (convex function). A function f : Rp → R is said to be convex if

∀(x, y) ∈ Rp × Rp , ∀α ∈ [0, 1], f (αx + (1 − α)y) ≤ αf (x) + (1 − α)f (y).

A function f : Rp → R is said to be strictly convex if

∀x 6= y, ∀α ∈ (0, 1), f (αx + (1 − α)y) ≤ αf (x) + (1 − α)f (y).

Definition 11. Let f : Rp → R be a convex function. A subgradient of f at x is any vector u ∈ Rp satisfying

∀y ∈ Rp , f (y) − f (x) ≥ uT (y − x).

The subdifferential of f at x, noted ∂f (x), is the set of all subgradients of f at x.

By simply using the definition of the subdifferential, we obtain the following caracterization of minimum
points (often refereed to as the Fermat’s rule). This is useful in deriving the LASSO first-order conditions.
Proposition 28. if f : Rp → R is a convex function, the point x? is a minimum if and only if 0 ∈ ∂f (x? ).
For differentiable (convex) function, the notion of subgradient coincides with the notion of gradient. This
is stated in the following.
Definition 12 (differentiable function). A function f : Rp → R is said to be differentiable at x0 if there
exists a vector u ∈ Rp such that
|f (x0 + h) − f (x0 ) − uT h|
lim exists
h→0 khk

As a consequence of the previous definition (when taking h = tek , t → 0) the partial derivatives (gradient)
exists for differentiable functions. The gradient of a differentiable function f at x ∈ Rp is denoted by ∇f (x).

Figure D.1: Draws of subgradients.

45
Proposition 29 (differentiable function and convexity). A differentiable function f : Rp → R is convex if
and only if

∀(x, y) ∈ Rp × Rp , f (y) − f (x) ≥ ∇f (x)T (y − x).

Moreover, for any x ∈ Rp , the gradient ∇f (x) is the only vector satisfying the previous equation for any y.
Consequently, for differentiable and convex functions ∂f (x) = {∇f (x)}.

46
Bibliography

Boucheron, S., G. Lugosi, and P. Massart (2013). Concentration inequalities: A nonasymptotic theory of
independence. Oxford university press.

Efron, B., T. Hastie, I. Johnstone, R. Tibshirani, et al. (2004). Least angle regression. The Annals of
statistics 32 (2), 407–499.
Hastie, T., R. Tibshirani, and M. Wainwright (2015). Statistical learning with sparsity: the lasso and
generalizations. CRC press.

Knight, K. and W. Fu (2000). Asymptotics for lasso-type estimators. Annals of statistics, 1356–1378.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical
Society. Series B (Methodological), 267–288.
Tibshirani, R. J. (2013). The lasso problem and uniqueness. Electronic Journal of Statistics 7, 1456–1490.

Zhang, T. (2009). Adaptive forward-backward greedy algorithm for sparse learning with linear models. In
Advances in Neural Information Processing Systems, pp. 1921–1928.
Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American statistical associa-
tion 101 (476), 1418–1429.

Zou, H. and T. Hastie (2005). Regularization and variable selection via the elastic net. Journal of the Royal
Statistical Society: Series B (Statistical Methodology) 67 (2), 301–320.

47