III Modern Statistical Methods (Full)

Part III — Modern Statistical Methods

Based on lectures by R. D. Shah

Notes taken by Dexter Chua

Michaelmas 2017

These notes are not endorsed by the lecturers, and I have modified them (often

significantly) after lectures. They are nowhere near accurate representations of what

was actually lectured, and in particular, all errors are almost surely mine.

The remarkable development of computing power and other technology now allows

scientists and businesses to routinely collect datasets of immense size and complexity.

Most classical statistical methods were designed for situations with many observations

and a few, carefully chosen variables. However, we now often gather data where we

have huge numbers of variables, in an attempt to capture as much information as we

can about anything which might conceivably have an influence on the phenomenon

of interest. This dramatic increase in the number variables makes modern datasets

strikingly different, as well-established traditional methods perform either very poorly,

or often do not work at all.

Developing methods that are able to extract meaningful information from these large

and challenging datasets has recently been an area of intense research in statistics,

machine learning and computer science. In this course, we will study some of the

methods that have been developed to analyse such datasets. We aim to cover some of

the following topics.

–

Kernel machines: the kernel trick, the representer theorem, support vector

machines, the hashing trick.

– Penalised regression: Ridge regression, the Lasso and variants.

– Graphical modelling: neighbourhood selection and the graphical Lasso. Causal

inference through structural equation modelling; the PC algorithm.

–

High-dimensional inference: the closed testing procedure and the Benjamini–

Hochberg pro cedure; the debiased Lasso

Pre-requisites

Basic knowledge of statistics, probability, linear algebra and real analysis. Some

background in optimisation would be helpful but is not essential.

Contents

0 Introduction

1 Classical statistics

2 Kernel machines

2.1 Ridge regression

2.2 v-fold cross-validation

2.3 The kernel trick

2.4 Making predictions

2.5 Other kernel machines

2.6 Large-scale kernel machines

3 The Lasso and beyond

3.1 The Lasso estimator

3.2 Basic concentration inequalities

3.3 Convex analysis and optimization theory

3.4 Properties of Lasso solutions

3.5 Variable selection

3.6 Computation of Lasso solutions

3.7 Extensions of the Lasso

4 Graphical modelling

4.1 Conditional independence graphs

4.2 Structural equation modelling

4.3 The PC algorithm

5 High-dimensional inference

5.1 Multiple testing

5.2 Inference in high-dimensional regression

0 Introduction

In recent years, there has been a rather significant change in what sorts of data

we have to handle and what questions we ask about them, witnessed by the

popularity of the buzzwords “big data” and “machine learning”. In classical

statistics, we usually have a small set of parameters, and a very large data set.

We then use the large data set to estimate the parameters.

However, nowadays we often see scenarios where we have a very large number

of parameters, and the data set is relatively small. If we tried to apply our

classical linear regression, then we would just be able to tune the parameters so

that we have a perfect fit, and still have great freedom to change the parameters

without affecting the fitting.

One example is that we might want to test which genomes are responsible for

a particular disease. In this case, there is a huge number of genomes to consider,

and there is good reason to believe that most of them are irrelevant, i.e. the

parameters should be set to zero. Thus, we want to develop methods that find

the “best” fitting model that takes this into account.

Another problem we might encounter is that we just have a large data set,

and doing linear regression seems a bit silly. If we have so much data, we might

as well try to fit more complicated curves, such as polynomial functions and

friends. Perhaps more ambitiously, we might try to find the best continuously

differentiable function that fits the curve, or, as analysts will immediately suggest

as an alternative, weakly differentiable functions.

There are many things we can talk about, and we can’t talk about all of

them. In this course, we are going to cover 4 different topics of different size:

– Kernel machines

– The Lasso and its extensions

– Graphical modeling and causal inference

– Multiple testing and high-dimensional inference

The four topics are rather disjoint, and draw on different mathematical skills.

1 Classical statistics

This is a course on modern statistical methods. Before we study methods, we

give a brief summary of what we are not going to talk about, namely classical

statistics.

So suppose we are doing regression. We have some predictors

∈ R

and

responses

∈ R

, and we hope to find a model that describes

as a function of

x. For convenience, define the vectors

X =













, Y =













The linear model then assumes there is some β

∈ R

such that

Y = Xβ

+ ε,

where

is some (hopefully small) error random variable. Our goal is then to

estimate β

given the data we have.

has full column rank, so that

is invertible, then we can use ordinary

least squares to estimate β

, with estimate

OLS

= argmin

β∈R

kY − xβk

= (X

−1

This assumes nothing about

itself, but if we assume that

Eε

= 0 and

var

(

) =

I, then this estimate satisfies

– E

OLS

= (X

−1

Xβ

= β

– var

(

OLS

) = (X

−1

var(ε)X(X

−1

= σ

−1

In particular, this is an unbiased estimator. Even better, this is the best linear

unbiased estimator. More precisely, the Gauss–Markov theorem says any other

linear estimator

β = AY has var(

β) − var(

OLS

) positive semi-definite.

Of course, ordinary least squares is not the only way to estimate

. Another

common method for estimating parameters is maximum likelihood estimation,

and this works for more general models than linear regression. For people who are

already sick of meeting likelihoods, this will be the last time we meet likelihoods

in this course.

Suppose we want to estimate a parameter

via knowledge of some data

We assume Y has density f(y; θ). We define the log-likelihood by

`(θ) = log f(Y, θ).

The maximum likelihood estimator then maximize `(θ) over θ to get

θ.

Similar to ordinary least squares, there is a theorem that says maximum

likelihood estimation is the “best”. To do so, we introduce the Fisher information

matrix. This is a family of d × d matrices indexed by θ, defined by

(θ) = −E



∂

∂θ

`(θ)



The relevant theorem is

Theorem (Cram´er–Rao bound). If

is an unbiased estimator for

, then

var(

θ) −I

−1

(θ) is positive semi-definite.

Moreover, asymptotically, as

n → ∞

, the maximum likelihood estimator is

unbiased and achieves the Carm´er–Rao bound.

Another wonderful fact about the maximum likelihood estimator is that

asymptotically, it is normal distributed, and so it is something we understand

well.

This might seem very wonderful, but there are a few problems here. The

results we stated are asymptotic, but what we actually see in real life is that

n → ∞

, the value of

also increases. In these contexts, the asymptotic

property doesn’t tell us much. Another issue is that these all talk about unbiased

estimators. In a lot of situations of interest, it turns out biased methods do

much much better than these methods we have.

Another thing we might be interested is that as

gets large, we might want

to use more complicated models than simple parametric models, as we have

much more data to mess with. This is not something ordinary least squares

provides us with.

2 Kernel machines

We are going to start a little bit slowly, and think about our linear model

Xβ

, where

(

) = 0 and

var

(

) =

. Ordinary least squares is an

unbiased estimator, so let’s look at biased estimators.

For a biased estimator,

, we should not study the variance, but the mean

squared error

E[(

β − β

)(

β − β

)

] = E(

β − E

β + E

β − β

)(

β − E

β + E

β − β

)

= var(

β) + (E

β − β

)(E

β − β

)

The first term is, of course, just the variance, and the second is the squared bias.

So the point is that if we pick a clever biased estimator with a tiny variance,

then this might do better than unbiased estimators with large variance.

2.1 Ridge regression

Ridge regression was introduced in around 1970. The idea of Ridge regression is

to try to shrink our estimator a bit in order to lessen the variance, at the cost of

introducing a bias. We would hope then that this will result in a smaller mean

squared error.

Definition (Ridge regression). Ridge regression solves

(ˆµ

) = argmin

(µ,β)∈R×R

{kY − µ1 − Xβk

+ λkβk

where 1 is a vector of all 1’s. Here

λ ≥

0 is a tuning parameter, and it controls

how much we penalize a large choice of β.

Here we explicitly have an intercept term. Usually, we eliminate this by

adding a column of 1’s in

. But here we want to separate that out, since

we do not want give a penalty for large

. For example, if the parameter is

temperature, and we decide to measure in degrees Celsius rather than Kelvins,

then we don’t want the resulting ˆµ,

β to change.

More precisely, our formulation of Ridge regression is so that if we make the

modification

= c1 + Y,

then we have

ˆµ

) = ˆµ

(Y ) + c.

Note also that the Ridge regression formula makes sense only if each entry in

have the same order of magnitude, or else the penalty will only have a significant

effect on the terms of large magnitude. Standard practice is to subtract from

each column of

its mean, and then scale it to have

norm

√

. The actual

number is not important here, but it will in the case of the Lasso.

By differentiating, one sees that the solution to the optimization problem is

ˆµ

Y =

= (X

X + λI)

−1

Note that we can always pick

such that the matrix (

λI

) is invertible.

In particular, this can work even if we have more parameters than data points.

This will be important when we work with the Lasso later on.

If some god-like being told us a suitable value of

to use, then Ridge

regression always works better than ordinary least squares.

Theorem. Suppose

rank

(

) =

. Then for

λ >

0 sufficiently small (depending

on β

and σ

), the matrix

OLS

− β

)(

OLS

− β

)

− E(

− β

)(

− β

)

(∗)

is positive definite.

Proof.

We know that the first term is just

(

)

−1

. The right-hand-side

has a variance term and a bias term. We first look at the bias:

β − β

] = (X

X + λI)

−1

Xβ

− β

= (X

X + λI)

−1

X + λI − λI)β

− β

= −λ(X

X + λI)

−1

We can also compute the variance

var(

β) = σ

X + λI)

−1

X(X

X + λI)

−1

Note that both terms appearing in the squared error look like

X + λI)

−1

something(X

X + λI)

−1

So let’s try to write σ

−1

in this form. Note that

−1

= (X

X + λI)

−1

X + λI)(X

−1

X + λI)(X

X + λI)

−1

= (X

X + λI)

−1

X + 2λI + λ

−1

)(X

X + λI)

−1

Thus, we can write (∗) as

X + λI)

−1



X + 2λI + λ

−1

)

− σ

X − λ

(β

)



X + λI)

−1

= λ(X

X + λI)

−1



2σ

I + λ(X

−1

− λβ

(β

)



X + λI)

−1

Since λ > 0, this is positive definite iff

2σ

I + σ

λ(X

−1

− λβ

(β

)

is positive definite, which is true for 0 < λ <

2σ

kβ

While this is nice, this is not really telling us much, because we don’t know

how to pick the correct

. It also doesn’t tell us when we should expect a big

improvement from Ridge regression.

To understand this better, we need to use the singular value decomposition.

Theorem (Singular value decomposition). Let

X ∈ R

n×p

be any matrix. Then

it has a singular value decomposition (SVD)

n×p

= U

n×n

n×p

p×p

where

U, V

are orthogonal matrices, and

≥ D

≥ ··· ≥ D

≥

0, where

m = min(n, p), and all other entries of D are zero.

When

n > p

, there is an alternative version where

is an

n × p

matrix

with orthogonal columns, and

is a

p × p

diagonal matrix. This is done by

replacing

by its first

columns and

by its first

rows. This is known as

the thin singular value decomposition. In this case,

but

is not

the identity.

Let’s now apply try to understand Ridge regressions a little better. Suppose

n > p. Then using the thin SVD, the fitted values from ridge regression are

= X(X

X + λI)

−1

= UDV

(V D

+ λI)

−1

V DU

We now note that

V D

λI

(

λI

)

, since

is still orthogonal.

We then have

(V (D

+ λI)V

)

−1

= V (D

+ λI)

−1

Since

and

λI

are both diagonal, it is easy to compute their inverses as well.

We then have

= UD

+ λI)

−1

Y =

j=1

+ λ

Here U

refers to the jth column of U.

Now note that the columns of

form an orthonormal basis for the column

space of

. If

= 0, then this is just a fancy formula for the projection of

onto the column space of

. Thus, what this formula is telling us is that we

look at this projection, look at the coordinates in terms of the basis given by

the columns of U , and scale accordingly.

We can now concretely see the effect of our

. The shrinking depends on the

size of D

, and the larger D

is, the less shrinking we do.

This is not very helpful if we don’t have a concrete interpretation of the

or rather, the columns of

. It turns out the columns of

are what are known

as the normalized principal components of X.

Take u ∈ R

, kuk

= 1. The sample variance of Xu ∈ R

is then

kXuk

Xu =

V D

We write w = V

u. Then kwk

= 1 since V is orthogonal. So we have

kXuk

w =

≤

and the bound is achieved when

, or equivalently,

V e

. Thus,

gives the coefficients of the linear combination of columns of

that has largest

sample variance. We can then write the result as

= UDV

= U

We can extend this to a description of the other columns of

, which is done

in the example sheet. Roughly, the subsequent principle components obey the

same optimality conditions with the added constraint of being orthogonal to all

earlier principle components.

The conclusion is that Ridge regression works best if

lies in the space

spanned by the large principal components of X.

2.2 v-fold cross-validation

In practice, the above analysis is not very useful, since it doesn’t actually tell

us what

we should pick. If we are given a concrete data set, how do we know

what λ we should pick?

One common method is to use

-fold cross-validation. This is a very general

technique that allows us to pick a regression method from a variety of options.

We shall explain the method in terms of Ridge regression, where our regression

methods are parametrized by a single parameter

, but it should be clear that

this is massively more general.

Suppose we have some data set (

, x

)

i=1

, and we are asked to predict the

value of

∗

given a new predictors

∗

. We may want to pick

to minimize the

following quantity:

∗

− (x

∗

)

(X, Y ))

| X, Y

It is difficult to actually do this, and so an easier target to minimize is the

expected prediction error

∗

− (x

∗

)

(

Y ))

One thing to note about this is that we are thinking of

and

as arbitrary

data sets of size

, as opposed to the one we have actually got. This might be a

more tractable problem, since we are not working with our actual data set, but

general data sets.

The method of cross-validation estimates this by splitting the data into

folds.

(1)

, Y

(1)

)

(v)

, Y

(v)

)

is called the observation indices, which is the set of all indices

so that the

jth data point lies in the ith fold.

We let (X

(−k)

, Y

(−k)

) be all data except that in the kth fold. We define

CV (λ) =

k=1

i∈A

− x

(−k)

, Y

(−k)

))

We write

for the minimizer of this, and pick

(

X, Y

) as our estimate.

This tends to work very will in practice.

But we can ask ourselves — why do we have to pick a single

to make the

estimate? We can instead average over a range of different

. Suppose we have

computed

on a grid of

-values

> λ

> ··· > λ

. Our plan is to take a

good weighted average of the λ. Concretely, we want to minimize

k=1

i∈A

−

i=1

(−k)

, Y

(−k)

)

over

w ∈

, ∞

)

. This is known as stacking. This tends to work better than

just doing

-fold cross-validation. Indeed, it must — cross-validation is just a

special case where we restrict

to be zero in all but one entries. Of course, this

method comes with some heavy computational costs.

2.3 The kernel trick

We now come to the actual, main topic of the chapter. We start with a very

simple observation. An alternative way of writing the fitted values from Ridge

regression is

X + λI)X

= X

(XX

+ λI).

One should be careful that the

λI

are different matrices on both sides, as they

have different dimensions. Multiplying inverses on both sides, we have

(XX

+ λI)

−1

= (X

X + λI)

−1

We can multiply the right-hand-side by

to obtain the

from Ridge regression,

and multiply on the left to obtain the fitted values. So we have

(XX

+ λI)

−1

Y = X(X

X + λI)

−1

Y = X

Note that

is a

p × p

matrix, and takes

(

) time to compute. On the

other hand,

is an

n × n

matrix, and takes

(

) time to compute. If

p  n

, then this way of computing the fitted values would be much quicker. The

key point to make is that the fitted values from Ridge regression only depends

on K = XX

(and Y ). Why is this important?

Suppose we believe we have a quadratic signal

= x

β +

k,`

+ ε

Of course, we can do Ridge regression, as long as we add the products

as predictors. But this requires

(

) many predictors. Even if we use the

way, this has a cost of

(

), and if we used the naive method, it would

require O(np

) operations.

We can do better than this. The idea is that we might be able to compute

K directly. Consider

(1 + x

)

= 1 + 2x

k,`

This is equal to the inner product between vectors of the form

(1,

√

, . . . ,

√

, x

, . . . , x

, x

, . . . , x

) (∗)

If we set

= (1+

)

and form

(

λI

)

−1

, this is equivalent to forming

ridge regression with (

∗

) as our predictors. Note that here we don’t scale our

columns to have the same

norm. This is pretty interesting, because computing

this is only

(

). We managed to kill a factor of

in this computation. The

key idea here was that the fitted values depend only on

, and not on the values

of x

itself.

Consider the general scenario where we try to predict the value of

given

a predictor

x ∈ X

. In general, we don’t even assume

has some nice linear

structure where we can do linear regression.

If we want to do Ridge regression, then one thing we can do is that we

can try to construct some map

X → R

for some

, and then run Ridge

regression using

{φ

(

)

}

as our predictors. If we want to fit some complicated,

non-linear model, then

can potentially be very huge, and this can be costly.

The observation above is that instead of trying to do this, we can perhaps find

some magical way of directly computing

= k(x

, x

) = hφ(x

), φ(x

)i.

If we can do so, then we can simply use the above formula to obtain the fitted

values (we shall discuss the problem of making new predictions later).

Since it is only the function

that matters, perhaps we can specify a

“regression method” simply by providing a suitable function

X × X → R

If we are given such a function

, when would it arise from a “feature map”

φ : X → R

More generally, we will find that it is convenient to allow for

to take values

in infinite-dimensional vector spaces instead (since we don’t actually have to

compute φ!).

Definition (Inner product space). an inner product space is a real vector space

H endowed with a map h·, ·i : H × H → R and obeys

– Symmetry: hu, vi = hv, ui

– Linearity: If a, b ∈ R, then hau + bw, vi = ahu, vi + bhw, vi.

– Positive definiteness: hu, ui ≥ 0 with hu, ui = 0 iff u = 0.

If we want

to come from a feature map

, then an immediate necessary

condition is that

has to be symmetric. There is also another condition that

corresponds to the positive-definiteness of the inner product.

Proposition. Given φ : X × X → H, define k : X × X → R by

k(x, x

) = hφ(x), φ(x

)i.

Then for any x

, . . . , x

∈ X, the matrix K ∈ R

× R

with entries

= k(x

, x

)

is positive semi-definite.

Proof. Let x

, . . . , x

∈ X, and α ∈ R

. Then

i,j

k(x

, x

)α

i,j

hφ(x

), φ(x

)iα

φ(x

)

≥ 0

since the inner product is positive definite.

This suggests the following definition:

Definition (Positive-definite kernel). A positive-definite kernel (or simply ker-

nel) is a symmetric map

k : X×X → R

such that for all

n ∈ N

and

, . . . , x

∈ X

the matrix K ∈ R

× R

with entries

= k(x

, x

)

is positive semi-definite.

We will prove that every positive-definite kernel comes from a feature map.

However, before that, let’s look at some examples.

Example. Suppose k

, k

, . . . , k

are kernels. Then

– If α

, α

≥ 0, then α

+ α

is a kernel. Moreover, if

k(x, x

) = lim

m→∞

(x, x

)

exists, then k is a kernel.

– The pointwise product k

is a kernel, where

)(x, x

) = k

(x, x

Example. The linear kernel is

(

x, x

) =

. To see this, we see that this is

given by the feature map φ = id and taking the standard inner product on R

Example. The polynomial kernel is k(x, x

) = (1 + x

)

for all d ∈ N.

We saw this last time with

= 2. This is a kernel since both 1 and

are

kernels, and sums and products preserve kernels.

Example. The Gaussian kernel is

k(x, x

) = exp



−

kx − x

2σ



The quantity σ is known as the bandwidth.

To show that it is a kernel, we decompose

kx − x

= kxk

+ kx

− 2x

We define

(x, x

) = exp



−

kxk

2σ



exp



−

2σ



This is a kernel by taking φ( ·) = exp



−

k· k

2σ



Next, we can define

(x, x

) = exp





∞

r=0





We see that this is the infinite linear combination of powers of kernels, hence is

a kernel. Thus, it follows that k = k

is a kernel.

Note also that any feature map giving this

must take values in an infinite

dimensional inner product space. Roughly speaking, this is because we have

arbitrarily large powers of x and x

in the expansion of k.

Example. The Sobolev kernel is defined as follows: we take

= [0

1], and let

k(x, x

) = min(x, x

A slick proof that this is a kernel is to notice that this is the covariance function

of Brownian motion.

Example. The Jaccard similarity is defined as follows: Take

to be the set of

all subsets of {1, . . . , p}. For x, x

∈ X, define

k(x, x

) =

(

|x∩x

|x∪x

x ∪ x

6= ∅

1 otherwise

Theorem (Moore–Aronszajn theorem). For every kernel

X × X → R

, there

exists an inner product space H and a feature map φ : X → H such that

k(x, x

) = hφ(x), φ(x

)i.

This is not actually the full Moore–Aronszajn theorem, but a simplified

version of it.

Proof. Let H denote the vector space of functions f : X → R of the form

f( ·) =

i=1

k( ·, x

) (∗)

for some n ∈ N, α

, . . . , α

∈ R and x

, . . . , x

∈ X. If

g( ·) =

i=1

k( ·, x

) ∈ H,

then we tentatively define the inner product of f and g by

hf, gi =

i=1

j=1

k(x

, x

We have to check that this is an inner product, but even before that, we need to

check that this is well-defined, since

and

can be represented in the form (

∗

)

in multiple ways. To do so, simply observe that

i=1

j=1

k(x

, x

) =

i=1

g(x

) =

j=1

f(x

). (†)

The first equality shows that the definition of our inner product does not depend

on the representation of

, while the second equality shows that it doesn’t depend

on the representation of f.

To show this is an inner product, note that it is clearly symmetric and bilinear.

To show it is positive definite, note that we have

hf, fi =

i=1

j=1

k(x

, x

)α

≥ 0

since the kernel is positive semi-definite. It remains to check that if

hf, fi

= 0,

then

= 0 as a function. To this end, note the important reproducing property:

by (†), we have

hk( ·, x), f i = f(x).

This says k(·, x) represents the evaluation-at-x linear functional.

In particular, we have

f(x)

= hk( ·, x), f i

≤ hk( ·, x), k( ·, x)ihf, f i = 0.

Here we used the Cauchy–Schwarz inequality, which, if you inspect the proof,

does not require positive definiteness, just positive semi-definiteness. So it follows

that f ≡ 0. Thus, we know that H is indeed an inner product space.

We now have to construct a feature map. Define φ : X → H by

φ(x) = k( ·, x).

Then we have already observed that

hφ(x), φ(x

)i = hk( ·, x), k( ·, x

)i = k(x, x

as desired.

We shall boost this theorem up to its full strength, where we say there is

a “unique” inner product space with such a feature map. Of course, it will not

be unique unless we impose appropriate quantifiers, since we can just throw in

some random elements into H.

Recall (or learn) from functional analysis that any inner product space

a normed space, with norm

kfk

= hf, fi

We say sequence (

) in a normed space

is Cauchy if

− f

→

0 as

n, m → ∞

. An inner product space is complete if every Cauchy sequence has a

limit (in the space), and a complete inner product space is called a Hilbert space.

One important property of a Hilbert space

is that if

is a closed subspace

, then every

f ∈ B

has a unique representation as

, where

v ∈ V

and

w ∈ V

⊥

= {u ∈ B : hu, yi = 0 for all y ∈ V }.

By adding limits of Cauchy sequences to

, we can obtain a Hilbert space.

Indeed, if (f

) is Cauchy in H, then

(x) − f

(x)| ≤ k

1/2

(x, x)kf

− f

→ 0

m, n → ∞

. Since every Cauchy sequence in

converges (i.e.

is complete),

it makes sense to define a limiting function

f(x) = lim

n→∞

(x),

and it can be shown that after augmenting

with such limits, we obtain a

Hilbert space. In fact, it is a special type of Hilbert space.

Definition (Reproudcing kernel Hilbert space (RKHS)). A Hilbert space

functions

X → R

is a reproducing kernel Hilbert space if for each

x ∈ X

there exists a k

∈ B such that

, fi = f (x)

for all x ∈ B.

The function k : X × X → R given by

k(x, x

) = hk

, k

i = k

) = k

(x)

is called the reproducing kernel associated with B.

By the Riesz representation theorem, this condition is equivalent to saying

that pointwise evaluation is continuous.

We know the reproducing kernel associated with an RKHS is a positive-

definite kernel, and the Moore–Aronszajn theorem can be extended to show

that to each kernel

, there is a unique representing kernel Hilbert space

and

feature map φ : X → H such that k(x, x

) = hφ(x), φ(x

)i.

Example. Take the linear kernel

k(x, x

) = x

By definition, we have

H =

(

f : R

→ R | f(x) =

i=1

)

for some n ∈ N, x

, . . . , x

∈ R

. We then see that this is equal to

H = {f : R

→ R | f(x) = β

x for some β ∈ R

and if f (x) = β

x, then kf k

= k(β, β) = kβk

Example. Take the Sobolev kernel, where

= [0

1] and

(

x, x

) =

min

(

x, x

Then

includes all linear combinations of functions of the form

x 7→ min

(

x, x

where x

∈ [0, 1], and their pointwise limits. These functions look like

min(x, x

)

Since we allow arbitrary linear combinations of these things and pointwise limits,

this gives rise to a large class of functions. In particular, this includes all Lipschitz

functions that are 0 at the origin.

In fact, the resulting space is a Sobolev space, with the norm given by

kfk =



(x)



1/2

For example, if we take f = min( ·, x), then by definition we have

kmin( ·, x)k

= min(x, x) = x,

whereas the formula we claimed gives

dt = x.

In general, it is difficult to understand the RKHS, but if we can do so, this

provides a lot of information on what we are regularizing in the kernel trick.

2.4 Making predictions

Let’s now try to understand what exactly we are doing when we do ridge

regression with a kernel

. To do so, we first think carefully what we were doing

in ordinary ridge regression, which corresponds to using the linear kernel. For

the linear kernel, the

norm

kβk

corresponds exactly to the norm in the

RKHS kf k

. Thus, an alternative way of expressing ridge regression is

f = argmin

f∈H

(

i=1

− f(x

))

+ λkfk

)

, (∗)

where

is the RKHS of the linear kernel. Now this way of writing ridge

regression makes sense for an arbitrary RKHS, so we might think this is what

we should solve in general.

But if

is infinite dimensional, then naively, it would be quite difficult to

solve (∗). The solution is provided by the representer theorem.

Theorem (Representer theorem). Let

be an RKHS with reproducing kernel

. Let

be an arbitrary loss function and

: [0

, ∞

)

→ R

any strictly increasing

function. Then the minimizer

f ∈ H of

(f) = c(Y, x

, . . . , x

, f(x

), . . . , f(x

)) + J(kfk

)

lies in the linear span of {k( ·, x

)}

i=1

Given this theorem, we know our optimal solution can be written in the form

f( ·) =

i=1

ˆα

k( ·, x

and thus we can rewrite our optimization problem as looking for the

ˆα ∈ R

that minimizes

(α) = c(Y, x

, . . . , x

, Kα) + J(α

Kα),

over α ∈ R

(with K

= k(x

, x

)).

For an arbitrary

, this can still be quite difficult, but in the case of Ridge

regression, this tells us that (∗) is equivalent to minimizing

kY − Kαk

+ λα

Kα,

and by differentiating, we find that this is given by solving

K ˆα = K(K + λI)

−1

Of course

K ˆα

is also our fitted values, so we see that if we try to minimize

(∗), then the fitted values is indeed given by K(K + λI)

−1

Y , as in our original

motivation.

Proof. Suppose

f minimizes Q

. We can then write

f = u + v

where u ∈ V = span{k( ·, x

) : i = 1, . . . , n} and v ∈ V

⊥

. Then

f(x

) = hf, k( ·, x

)i = hu + v, k( ·, x

)i = hu, k( ·, x

)i = u(x

So we know that

c(Y, x

, . . . , x

, f(x

), . . . , f(x

)) = c(Y, x

, . . . , x

, u(x

), . . . , u(x

)).

Meanwhile,

kfk

= ku + vk

= kuk

+ kvk

using the fact that u and v are orthogonal. So we know

J(kfk

) ≥ J(kuk

)

with equality iff

= 0. Hence

(

)

≥ Q

(

) with equality iff

= 0, and so we

must have v = 0 by optimality.

Thus, we know that the optimizer in fact lies in

, and then the formula of

just expresses Q

in terms of α.

How well does our kernel machine do? Let

be an RKHS with reproducing

kernel k, and f

∈ H. Consider a model

= f

) + ε

for i = 1, . . . , n, and assume Eε = 0, var(ε) = σ

For convenience, We assume

≤

1. There is no loss in generality by

assuming this, since we can always achieve this by scaling the norm.

Let

be the kernel matrix

(

, x

) with eigenvalues

≥ d

≥ ··· ≥

≥ 0. Define

= argmin

f∈H

(

i=1

− f(x

))

+ λkfk

)

Theorem. We have

i=1

E(f

) −

))

≤

i=1

+ λ)

≤

i=1

min



, λ



Given a data set, we can compute the eigenvalues

, and thus we can compute

this error bound.

Proof. We know from the representer theorem that

(

), . . . ,

))

= K(K + λI)

−1

Also, there is some α ∈ R

such that

), . . . , f

))

= Kα.

Moreover, on the example sheet, we show that

1 ≥ kf

≥ α

Kα.

Consider the eigen-decomposition

UDU

, where

is orthogonal,

and D

= 0 for i 6= j. Then we have

i=1

) −

))

= EkKα − K(K + λI)

−1

(Kα + ε)k

Noting that Kα = (K + λI)(K + λI)

−1

Kα, we obtain

= E



λ(K + λI)

−1

Kα − K(K + λI)

−1



= λ



(K + λI)

−1

Kα



| {z }

(I)

+ E



K(K + λI)

−1



| {z }

(II)

At this stage, we can throw in the eigen-decomposition of K to write (I) as

(I) = λ



(UDU

+ λI)

−1

UDU



= λ

kU(D + λI)

−1

| {z }

i=1

+ λ)

Now we have

Kα = α

UDU

α = α

UDD

where D

is diagonal and

(

−1

> 0

0 otherwise

We can then write this as

Kα =

The key thing we know about this is that it is ≤ 1.

Note that by definition of

, we see that if

= 0, then

= 0. So we can

write

(II) =

i:d

≥0

+ λ)

≤ λ max

i=1,...,n

+ λ)

by H¨older’s inequality with (p, q) = (1, ∞). Finally, use the inequality that

(a + b)

≥ 4ab

to see that we have

(I) ≤

So we are left with (II), which is a bit easier. Using the trace trick, we have

(II) = Eε

(K + λI)

−1

(K + λI)

−1

= E tr



K(K + λI)

−1

εε

(K + λI)

−1



= tr



K(K + λI)

−1

E(εε

)(K + λI)

−1



= σ



(K + λI)

−2



= σ



(UDU

+ λI)

−2



= σ



(D + λI)

−2



= σ

i=1

+ λ)

Finally, writing

+λ)

, we have

+ λ)

≤ min



4λ



and we have the second bound.

How can we interpret this? If we look at the formula, then we see that we can

have good benefits if the

’s decay quickly, and the exact values of the

depend

only on the choice of

. So suppose the

are random and idd and independent

of ε. Then the entire analysis is still valid by conditioning on x

, . . . , x

We define ˆµ

, and λ

. Then we can rewrite our result to say

i=1

) −

))

≤

i=1

E min



ˆµ

, λ



≡ Eδ

(λ

We want to relate this more directly to the kernel somehow. Given a density

p(x) on X, Mercer’s theorem guarantees the following eigenexpansion

k(x, x

) =

∞

j=1

(x)e

where the eigenvalues µ

and eigenfunctions e

obey

) =

k(x, x

(x)p(x) dx

and

(x)e

(x)p(x) dx = 1

j=k

One can then show that

i=1

min



ˆµ

, λ



≤

∞

i=1

min



, λ



up to some absolute constant factors. For a particular density

(

) of the input

space, we can try to solve the integral equation to figure out what the

’s are,

and we can find the expected bound.

When

is the Sobolev kernel and

(

) is the uniform density, then it turns

out we have

(2j − 1)

By drawing a picture, we see that we can bound

∞

i=1

min



, λ



∞

i=1

min



, λ



≤

∞

∧

(2j − 1)

dj ≤

√

So we find that

E(δ

(λ

)) = O



nλ

1/2

+ λ



We can then find the λ

that minimizes this, and we see that we should pick

∼





2/3

which gives an error rate of ∼





2/3

2.5 Other kernel machines

Recall that the representer theorem was much more general than what we used

it for. We could apply it to any loss function. In particular, we can apply it

to classification problems. For example, we might want to predict whether an

email is a spam. In this case, we have Y ∈ {−1, 1}

Suppose we are trying to find a hyperplane that separates the two sets

}

and

}

i:Y

=−1

. Consider the really optimistic case where they can

indeed be completely separated by a hyperplane through the origin. So there

exists

β ∈ R

(the normal vector) with

β >

0 for all

. Moreover, it would

be nice if we can maximize the minimum distance between any of the points and

the hyperplane defined by β. This is given by

max

kβk

Thus, we can formulate our problem as

maximize M among β ∈ R

, M ≥ 0 subject to

kβk

≥ M.

This optimization problem gives the hyperplane that maximizes the margin

between the two classes.

Let’s think about the more realistic case where we cannot completely separate

the two sets. We then want to impose a penalty for each point on the wrong

side, and attempt to minimize the distance.

There are two options. We can use the penalty

i=1



M −

kβk



where t

= t1

t≥0

. Another option is to take

i=1



1 −

Mkβk



which is perhaps more natural, since now everything is “dimensionless”. We

will actually use neither, but will start with the second form and massage it to

something that can be attacked with our existing machinery. Starting with the

second option, we can write our optimization problem as

max

M≥0,β∈R

M − λ

i=1



1 −

Mkβk



Since the objective function is invariant to any positive scaling of

, we may

assume kβk

, and rewrite this problem as maximizing

kβk

− λ

i=1

(1 − Y

β)

Replacing

max

kβk

with minimizing

kβk

and adding instead of subtracting the

penalty part, we modify this to say

min

β∈R

kβk

+ λ

i=1

(1 − Y

β)

The final change we make is that we replace

with

, and multiply the whole

equation by λ to get

min

β∈R

λkβk

i=1

(1 − Y

β)

This looks much more like what we’ve seen before, with

λkβk

being the penalty

term and

i=1

(1 − Y

β)

being the loss function.

The final modification is that we want to allow planes that don’t necessarily

pass through the origin. To do this, we allow ourselves to translate all the

’s

by a fixed vector δ ∈ R

. This gives

min

β∈R

,δ∈R

λkβk

i=1



1 − Y

− δ)



Since

always appears together, we can simply replace it with a constant

and write our problem as

(ˆµ,

β) = argmin

(µ,β)∈R×R

i=1



1 − Y

β + µ)



+ λkβk

. (∗)

This gives the support vector classifier .

This is still just fitting a hyperplane. Now we would like to “kernelize” this,

so that we can fit in a surface instead. Let

be the RKHS corresponding to

the linear kernel. We can then write (∗) as

(ˆµ

) = argmin

(µ,f)∈R×H

i=1

(1 − Y

(f(x

) + µ))

+ λkfk

The representer theorem (or rather, a slight variant of it) tells us that the above

optimization problem is equivalent to the support vector machine

(ˆµ

, ˆα

) = argmin

(µ,α)∈R×R

i=1

(1 − Y

α + µ))

+ λα

Kα

where

(

, x

) and

is the reproducing kernel of

. Predictions at a

new x are then given by

sign

ˆµ

i=1

ˆα

λ,i

k(x, x

)

One possible alternative to this is to use logistic regression. Suppose that

log



P(Y

= 1)

P(Y

= −1)



= x

and that

, . . . , Y

are independent. An estimate of

can be obtained through

maximizing the likelihood, or equivalently,

argmin

b∈R

i=1

log(1 + exp(−Y

β)).

To kernelize this, we introduce an error term of

λkβk

, and then kernelize this.

The resulting optimization problem is then given by

argmin

f∈H

i=1

log(1 + exp(−Y

f(x

))) + λkfk

We can then solve this using the representer theorem.

2.6 Large-scale kernel machines

How do we apply kernel machines at large scale? Whenever we want to make a

prediction with a kernel machine, we need

(

) many steps, which is quite a pain

if we work with a large data set, say a few million of them. But even before that,

forming K ∈ R

n×n

and inverting K + λI or equivalent can be computationally

too expensive.

One of the more popular approaches to tackle these difficulties is to form a

randomized feature map

φ : X → R

such that

φ(x)

φ(x

) = k(x, x

)

To increase the quality of the approximation, we can use

x 7→

√

(

(x), . . . ,

(x))

∈ R

where the

(x) are iid with the same distribution as

φ.

Let Φ be the matrix with

th row

√

(

)

, . . . ,

(

)). When performing,

e.g. kernel ridge regression, we can compute

θ = (Φ

Φ + λI)

−1

The computational cost of this is

((

)

(

)

). The key thing of course is

that this is linear in

, and we can choose the size of

so that we have a good

trade-off between the computational cost and the accuracy of the approximation.

Now to predict a new x, we form

√

(

(x), . . . ,

(x))

θ,

and this is O(Lb).

In 2007, Rahimi and Recht proposed a random map for shift-invariant kernels,

i.e. kernels k such that k(x, x

) = h(x − x

) for some h (we work with X = R

A common example is the Gaussian kernel.

One motivation of this comes from a classical result known as Bochner’s

theorem.

Theorem (Bochner’s theorem). Let

× R

→ R

be a continuous kernel.

Then

is shift-invariant if and only if there exists some distribution

and c > 0 such that if W ∼ F , then

k(x, x

) = cE cos((x − x

)

W ).

Our strategy is then to find some

φ depending on W such that

c cos((x − x

)

W ) =

φ(x)

φ(x

We are going to take b = 1, so we don’t need a transpose on the right.

The idea is then to play around with trigonometric identities to try to write

c cos((x − x

)

W ). After some work, we find the following solution:

Let u ∼ U[−π, π], and take x, y ∈ R. Then

E cos(x + u) cos(y + u) = E(cos x cos u − sin x sin u)(cos y cos u − sin y sin u)

Since u has the same distribution as −u, we see that E cos u sin u = 0.

Also, cos

u + sin

u = 1. Since u ranges uniformly in [−π, π], by symmetry,

we have E cos

u = E sin

u =

. So we find that

2E cos(x + u) cos(y + u) = (cos x cos y + sin x sin y) = cos(x − y).

Thus, given a shift-invariant kernel

with associated distribution

, we set

W ∼ F independently of u. Define

φ(x) =

√

2c cos(W

x + u).

Then

φ(x)

φ(x

) = 2cE[E[cos(W

x + u) cos(W

+ u) | W ]]

= cE cos(W

(x − x

))

= k(x, x

In order to get this to work, we must find the appropriate

. In certain

cases, this W is actually not too hard to find:

Example. Take

k(x, x

) = exp



−

2σ

kx − x



the Gaussian kernel. If W ∼ N(0, σ

−2

I), then

= e

−ktk

/(2σ

)

= E cos(t

W ).

3 The Lasso and beyond

We are interested in the situation where we have a very large number of variables

compared to the size of the data set. For example the data set might be all the

genetic and environmental information about patients, and we want to predict if

they have diabetes. The key property is that we expect that most of the data

are not useful, and we want to select a small subset of variables that are relevant,

and form prediction models based on them.

In these cases, ordinary least squares is not a very good tool. If there are

more variables than the number of data points, then it is likely that we can

pick the regression coefficients so that our model fits our data exactly, but this

model is likely to be spurious, since we are fine-tuning our model based on the

random fluctuations of a large number of irrelevant variables. Thus, we want to

find a regression method that penalizes non-zero coefficients. Note that Ridge

regression is not good enough — it encourages small coefficients, but since it

uses the

norm, it’s quite lenient towards small, non-zero coefficients, as the

square of a small number is really small.

There is another reason to favour models that sets a lot of coefficients to

zero. Consider the linear model

Y = Xβ

+ ε,

where as usual

Eε

= 0 and

var

(

) =

. Let

= 0

}

, and let

|S|

As before, we suppose we have some a priori belief that s  p.

For the purposes of this investigation, suppose

has full column rank, so

that we can use ordinary least squares. Then the prediction error is

EkXβ

− X

OLS

− β

)

OLS

− β

)

E tr(

OLS

− β

)(

OLS

− β

)

tr E(

OLS

− β

)(

OLS

− β

)

tr Eσ

−1

= σ

Note that this does not depend on what what β

is, but only on σ

, p and n.

In the situation we are interested in, we expect

s  p

. So if we can find

and find ordinary least squares just on these, then we would have a mean

squared prediction error of σ

, which would be much much smaller.

We first discuss a few classical model methods that does this.

The first approach we may think of is the best subsets method, where we try

to do regression on all possible choices of

and see which is the “best”. For any

set

of indices, we let

be the submatrix of

formed from the columns of

with index in

. We then regress

for every

M ⊆ {

, . . . , p}

, and

then pick the best model via cross-validation, for example. A big problem with

this is that the number of subsets grows exponentially with

, and so becomes

infeasible for, say, p > 50.

Another method might be forward selection. We start with an intercept-only

model, and then add to the existing model the predictor that reduces the RSS

the most, and then keep doing this until a fixed number of predictors have been

added. This is quite a nice approach, and is computationally quite fast even

is large. However, this method is greedy, and if we make a mistake at the

beginning, then the method blows up. In general, this method is rather unstable,

and is not great from a practical perspective.

3.1 The Lasso estimator

The Lasso (Tibshirani, 1996) is a seemingly small modification of Ridge regression

that solves our problems. It solves

(ˆµ

) = argmin

(µ,β)∈R×R

kY − µ1 − Xβk

+ λkβk

The key difference is that we use an `

norm on β rather than the `

norm,

kβk

k=1

|β

This makes it drastically different from Ridge regression. We will see that for

large, it will make all of the entries of

exactly zero, as opposed to being very

close to zero.

We can compare this to best subset regression, where we replace

kβk

with

something of the form

k=1

. But the beautiful property of the Lasso

is that the optimization problem is now continuous, and in fact convex. This

allows us to solve it using standard convex optimization techniques.

Why is the

norm so different from the

norm? Just as in Ridge regression,

we may center and scale

, and center

, so that we can remove

from the

objective. Define

(β) =

kY − Xβk

+ λkβk

Any minimizer

of Q

(β) must also be a solution to

min kY − Xβk

subject to kβk

≤ k

Similarly,

is a solution of

min kY − Xβk

subject to kβk

≤ k

So imagine we are given a value of

, and we try to solve the above

optimization problem with pictures. The region

kβk

≤ k

is given by a

rotated square

On the other hand, the minimum of

kY − Xβk

is at

OLS

, and the contours

are ellipses centered around this point.

To solve the minimization problem, we should pick the smallest contour that hits

the square, and pick the intersection point to be our estimate of β

. The point

is that since the unit ball in the

-norm has these corners, this

is likely to

be on the corners, hence has a lot of zeroes. Compare this to the case of Ridge

regression, where the constraint set is a ball:

Generically, for Ridge regression, we would expect the solution to be non-zero

everywhere.

More generally, we can try to use the `

norm given by

kβk

k=1

1/q

We can plot their unit spheres, and see that they look like

q = 0.5 q = 1 q = 1.5 q = 2

We see that

= 1 is the smallest value of

for which there are corners, and

also the largest value of

for which the constraint set is still convex. Thus,

= 1

is the sweet spot for doing regression.

But is this actually good? Suppose the columns of

are scaled to have

norm

√

n, and, after centering, we have a normal linear model

Y = Xβ

+ ε − ¯ε1,

with ε ∼ N

(0, σ

I).

Theorem. Let

β be the Lasso solution with

λ = Aσ

log p

for some A. Then with probability 1 − 2p

−(A

/2−1)

, we have

kXβ

− X

βk

≤ 4Aσ

log p

kβ

Crucially, this is proportional to

log p

, and not just

. On the other hand,

unlike what we usually see for ordinary least squares, we have

√

, and not

We will later obtain better bounds when we make some assumptions on the

design matrix.

Proof.

We don’t really have a closed form solution for

, and in general it doesn’t

exist. So the only thing we can use is that it in fact minimizes

(

). Thus, by

definition, we have

kY − X

βk

+ λk

βk

≤

kY − Xβ

+ λkβ

We know exactly what Y is. It is Xβ

+ ε − ¯ε1. If we plug this in, we get

kXβ

− X

βk

≤

β − β

) + λkβ

− λk

βk

Here we use the fact that X is centered, and so is orthogonal to 1.

Now H¨older tells us

|ε

β − β

)| ≤ kX

εk

∞

β − β

We’d like to bound

εk

∞

, but it can be arbitrarily large since

is a Gaussian.

However, with high probability, it is small. Precisely, define

Ω =



εk

∞

≤ λ



In a later lemma, we will show later that

(Ω)

≥

−

−(A

/2−1)

. Assuming Ω

holds, we have

kXβ

− X

βk

≤ λk

β − β

k − λk

βk + λkβ

k ≤ 2λkβ

3.2 Basic concentration inequalities

We now have to prove the lemma we left out in the theorem just now. In this

section, we are going to prove a bunch of useful inequalities that we will later

use to prove bounds.

Consider the event Ω as defined before. Then



εk

∞

> λ



= P





[

j=1



ε| > λ







≤

j=1



ε| > λ



Now

ε ∼ N

). So we just have to bound tail probabilities of normal

random variables.

The simplest tail bound we have is Markov’s inequality.

Lemma (Markov’s inequality). Let

be a non-negative random variable. Then

P(W ≥ t) ≤

EW.

Proof. We have

W ≥t

≤ W.

The result then follows from taking expectations and then dividing both sides

by t.

While this is a very simple bound, it is actually quite powerful, because it

assumes almost nothing about

. This allows for the following trick: given any

strictly increasing function

R →

, ∞

) and any random variable

, we have

P(W ≥ t) = P(ϕ(W ) ≥ ϕ(t)) ≤

Eϕ(W )

ϕ(t)

So we get a bound on the tail probability of

for any such function. Even

better, we can minimize the right hand side over a class of functions to get an

even better bound.

In particular, applying this with ϕ(t) = e

αt

gives

Corollary (Chernoff bound). For any random variable W , we have

P(W ≥ t) ≤ inf

α>0

−αt

αW

Note that

αW

is just the moment generating function of

, which we can

compute quite straightforwardly.

We can immediately apply this when

is a normal random variable,

W ∼

N(0, σ

). Then

αW

= e

So we have

P(W ≥ t) ≤ inf

α>0

exp



− αt



= e

−t

/(2σ

)

Observe that in fact this tail bound works for any random variable whose moment

generating function is bounded above by e

Definition (Sub-Gaussian random variable). A random variable

is sub-

Gaussian (with parameter σ) if

α(W −EW )

≤ e

for all α ∈ R.

Corollary. Any sub-Gaussian random variable W with parameter σ satisfies

P(W ≥ t) ≤ e

−t

/2σ

. 

In general, bounded random variables are sub-Gaussian.

Lemma (Hoeffding’s lemma). If

has mean zero and takes values in [

a, b

then W is sub-Gaussian with parameter

b−a

Recall that the sum of two independent Gaussians is still a Gaussian. This

continues to hold for sub-Gaussian random variables.

Proposition. Let (

)

i=1

be independent mean-zero sub-Gaussian random

variables with parameters (

)

i=0

, and let

γ ∈ R

. Then

is sub-Gaussian

with parameter



(γ

)



1/2

Proof. We have

E exp

i=1

E exp (αγ

)

≤

i=1

exp





= exp

i=1

We can now prove our bound for the Lasso, which in fact works for any

sub-Gaussian random variable.

Lemma. Suppose (

)

i=1

are independent, mean-zero sub-Gaussian with com-

mon parameter σ. Let

λ = Aσ

log p

Let X be a matrix whose columns all have norm

√

n. Then



εk

∞

≤ λ



≥ 1 − 2p

−(A

/2−1)

In particular, this includes ε ∼ N

(0, σ

I).

Proof. We have



εk

∞

> λ



≤

j=1



ε| > λ



But ±

ε are both sub-Gaussian with parameter





1/2

√

Then by our previous corollary, we get

j=1



ε|

∞

> λ



≤ 2p exp



−

2σ



Note that we have the factor of 2 since we need to consider the two cases

ε > λ and −

ε > λ.

Plugging in our expression of λ, we write the bound as

2p exp



−

log p



= 2p

1−A

This is all we need for our result on the Lasso. We are going to go a bit

further into this topic of concentration inequalities, because we will need them

later when we impose conditions on the design matrix. In particular, we would

like to bound the tail probabilities of products.

Definition (Bernstein’s condition). We say that a random variable

satisfies

Bernstein’s condition with parameters (σ, b) where a, b > 0 if

E[|W − EW |

] ≤

k! σ

k−2

for k = 2, 3, . . ..

The point is that these bounds on the moments lets us bound the moment

generating function of W .

Proposition (Bernstein’s inequality). Let

, W

, . . . , W

be independent ran-

dom variables with

, and suppose each

satisfies Bernstein’s condition

with parameters (σ, b). Then

α(W

−µ)

≤ exp



1 − b|α|



for all |α| <

i=1

− µ ≥ t

≤ exp



−

2(σ

+ bt)



for all t ≥ 0.

Note that for large t, the bound goes as e

−t

instead of e

−t

Proof. For the first part, we fix i and write W = W

. Let |α| <

. Then

α(W

−µ)

∞

k=0



− µ|



≤ 1 +

∞

k=2

|α|

k−2

= 1 +

1 − |α|b

≤ exp



1 − b|α|



Then

E exp

i=1

α(W

− µ)

i=1

E exp



− µ)



≤ exp





1 − b



assuming



. So it follows that

i=1

− µ ≥ t

≤ e

−αt

exp





1 − b



Setting

bt + σ

∈





gives the result.

Lemma. Let

W, Z

be mean-zero sub-Gaussian random variables with parameters

and

respectively. Then

W Z

satisfies Bernstein’s condition with parameter

(8σ

, 4σ

Proof.

For any random variable

(which we will later take to be

W Z

), for

k > 1, we know

E|Y − EY |

= 2



Y −



≤ 2



|Y | +

|EY |



Note that



|Y | +

|EY |



≤

|Y |

+ |EY |

by Jensen’s inequality. Applying Jensen’s inequality again, we have

|EY |

≤ E|Y |

Putting the whole thing together, we have

E|Y − EY |

≤ 2

E|Y |

Now take Y = W Z. Then

E|W Z − EW Z| ≤ 2

E|W Z|

≤ 2





1/2

(EZ

)

1/2

by the Cauchy–Schwarz inequality.

We know that sub-Gaussians satisfy a bound on the tail probability. We can

then use this to bound the moments of W and Z. First note that

∞

x<W

dx.

Taking expectations of both sides, we get

∞

P(x < W

) dx.

Since we have a tail bound on

instead of

, we substitute

. Then

dx = 2kt

2k−1

dt. So we get

= 2k

∞

2k−1

P(|W | > t) dt

= 4k

∞

exp



−

2σ



dt.

where again we have a factor of 2 to account for both signs. We perform yet

another substitution

x =

2σ

, dx =

dt.

Then we get

= 2

k+1

kσ

∞

k−1

−x

dx = 4 · k!σ

Plugging this back in, we have

E|W Z − EW Z|

≤ 2

k+1

k!σ

k!2

2k+2

k!(8σ

)

(4σ

)

k−2

3.3 Convex analysis and optimization theory

We’ll leave these estimates aside for a bit, and give more background on convex

analysis and convex optimization. Recall the following definition:

Definition (Convex set). A set

A ⊆ R

is convex if for any

x, y ∈ A

and

t ∈ [0, 1], we have

(1 − t)x + ty ∈ A.

non-convex convex

We are actually more interested in convex functions. We shall let our functions

to take value

∞

, so let us define

R ∪ {∞}

. The point is that if we want our

function to be defined in [

a, b

], then it is convenient to extend it to be defined

on all of R by setting the function to be ∞ outside of [a, b].

Definition (Convex function). A function f : R

→

R is convex iff

f((1 − t)x + ty) ≤ (1 − t)f(x) + tf (y)

for all

x, y ∈ R

and

t ∈

1). Moreover, we require that

(

)

< ∞

for at least

one x.

We say it is strictly convex if the inequality is strict for all

x, y

and

t ∈

1).

x y

(1 − t)x + ty

(1 − t)f (x) + tf(y)

Definition (Domain). Define the domain of a function f : R

→

R to be

dom f = {x : f(x) < ∞}.

One sees that the domain of a convex function is always convex.

Proposition.

(i)

Let

, . . . , f

→

be convex with

dom f

∩ ··· ∩ dom f

∅

, and

let c

, . . . , c

≥ 0. Then c

+ ··· + c

is a convex function.

(ii) If f : R

→ R is twice continuously differentiable, then

(a) f is convex iff its Hessian is positive semi-definite everywhere.

(b) f is strictly convex if its Hessian positive definite everywhere.

Note that having a positive definite Hessian is not necessary for strict con-

vexity, e.g. x

is strictly convex but has vanishing Hessian at 0.

Now consider a constrained optimization problem

minimize f (x) subject to g(x) = 0

where x ∈ R

and g : R

→ R

. The Lagrangian for this problem is

L(x, θ) = f(x) + θ

g(x).

Why is this helpful? Suppose

∗

is the minimum of

. Then note that for any

we have

inf

x∈R

L(x, θ) ≤ inf

x∈R

,g(x)=0

L(x, θ) = inf

x∈R

:g(x)=0

f(x) = c

∗

Thus, if we can find some

∗

, x

∗

such that

∗

minimizes

(

x, θ

∗

) and

(

∗

) = 0,

then this is indeed the optimal solution.

This gives us a method to solve the optimization problem — for each fixed

solve the unconstrained optimization problem

argmin

(

x, θ

). If we are doing

this analytically, then we would have a formula for

in terms of

. Then seek

for a θ such that g(x) = 0 holds.

Subgradients

Usually, when we have a function to optimize, we take its derivative and set it

to zero. This works well if our function is actually differentiable. However, the

norm is not a differentiable function, since

|x|

is not differentiable at 0. This

is not some exotic case we may hope to avoid most of the time — when solving

the Lasso, we actively want our solutions to have zeroes, so we really want to

get to these non-differentiable points.

Thus, we seek some generalized notion of derivative that works on functions

that are not differentiable.

Definition (Subgradient). A vector

v ∈ R

is a subgradient of a convex function

at x if f (y) ≥ f (x) + v

(y − x) for all y ∈ R

The set of subgradients of

is denoted

∂f

(

), and is called the subdif-

ferential.

f f

This is indeed a generalization of the derivative, since

Proposition. Let

be convex and differentiable at

x ∈ int

(

dom f

). Then

∂f(x) = {∇f(x)}.

The following properties are immediate from definition.

Proposition. Suppose

and

are convex with

int

(

dom f

)

∩ int

(

dom g

)

∅

and α > 0. Then

∂(αf)(x) = α∂f(x) = {αv : v ∈ ∂f (x)}

∂(f + g)(x) = ∂g(x) + ∂f(x).

The condition (for convex differentiable functions) that “

is a minimum iff

(x) = 0” now becomes

Proposition. If f is convex, then

∗

∈ argmin

x∈R

f(x) ⇔ 0 ∈ ∂f(x

∗

Proof.

Both sides are equivalent to the requirement that

(

)

≥ f

(

∗

) for all

We are interested in applying this to the Lasso. So we want to compute the

subdifferential of the `

norm. Let’s first introduce some notation.

Notation. For

x ∈ R

and

A ⊆ {

. . . , d}

, we write

for the sub-vector of

formed by the components of

induced by

. We write

−j

{j}

{1,...,d}\j

Similarly, we write x

−jk

= x

{jk}

etc.

We write

sgn(x

) =











−1 x

< 0

1 x

> 0

0 x

= 0

and sgn(x) = (sgn(x

), . . . , sgn(x

))

First note that k · k

is convex, as it is a norm.

Proposition. For x ∈ R

and A ∈ {j : x

6= 0}, we have

∂kxk

= {v ∈ R

: kvk

∞

≤ 1, v

= sgn(x

)}.

Proof.

It suffices to look at the subdifferential of the absolute value function,

and then add them up.

For

= 1

, . . . , d

, we define

→ R

by sending

. If

= 0, then

is differentiable at

, and so we know

∂g

(

) =

{sgn

(

)

}

, with

the

standard basis vector.

When x

= 0, if v ∈ ∂g

), then

(y) ≥ g

(x) + v

(y − x).

| ≥ v

(y − x).

We claim this holds iff

∈

[

−

1] and

−j

= 0. The

⇐

direction is an immediate

calculation, and to show

⇒

, we pick

−j

, and

= 0. Then we have

0 ≥ v

−j

So we know that v

−j

= 0. Once we know this, the condition says

| ≥ v

for all

. This is then true iff

∈

[

−

1]. Forming the set sum of the

∂g

(

)

gives the result.

3.4 Properties of Lasso solutions

Let’s now apply this to the Lasso. Recall that the Lasso objective was

(β) =

kY − Xβk

+ λkβk

We know this is a convex function in

. So we know

is a minimizer of

(

)

iff 0 ∈ ∂Q

(

β).

Let’s subdifferentiate Q

and see what this amounts to. We have

∂Q

(

β) =



−

(Y − Xβ)



+ λ

ˆν ∈ R

: kˆνk

∞

≤ 1, ˆν|

ˆρ

= sgn(

λ,ˆρ

)

where ˆρ

= {j :

λ,j

6= 0}.

Thus, 0 ∈ ∂Q

(

) is equivalent to

ˆν =

(Y − X

)

satisfying

kˆνk

∞

≤ 1, ˆν

ˆρ

= sgn(

λ,ˆρ

These are known as the KKT conditions for the Lasso.

In principle, there could be several solutions to the Lasso. However, at least

the fitted values are always unique.

Proposition. X

is unique.

Proof.

Fix

λ >

0 and stop writing it. Suppose

(1)

and

(2)

are two Lasso

solutions at λ. Then

(1)

) = Q(

(2)

) = c

∗

As Q is convex, we know

∗

= Q



(1)

(2)



≤

(1)

) +

(2)

) = c

∗

(1)

(2)

is also a minimizer.

Since the two terms in

(

) are individually convex, it must be the case that



(Y − X

(1)

) +

(Y − X

(2)

)



Y − X

(1)



Y − X

(2)



(

(1)

(2)

)



(1)

(2)

Moreover, since

k · k

is strictly convex, we can have equality only if

(1)

(2)

. So we are done.

Definition (Equicorrelation set). Define the equicorrelation set

to be the

set of k such that

(Y − X

)| = λ,

or equivalently, the

with

1, which is well-defined since it depends only

on the fitted values.

By the KKT conditions,

contains the set of non-zeroes of Lasso solution,

but may be strictly bigger than that.

Note that if rk(X

) = |

|, then the Lasso solution must be unique since

(

(1)

−

(2)

) = 0.

(1)

(2)

3.5 Variable selection

So we have seen that the Lasso picks out some “important” variables and discards

the rest. How well does it do the job?

For simplicity, consider a noiseless linear model

Y = Xβ

Our objective is to find the set

S = {k : β

6= 0}.

We may wlog assume

{

, . . . , s}

, and

{

. . . p}\S

(as usual

n ×p

We further assume that rk(X

) = s.

In general, it is difficult to find out

even if we know

|S|

. Under certain

conditions, the Lasso can do the job correctly. This condition is dictated by the

∞

norm of the quantity

∆ = X

)

−1

sgn(β

We can understand this a bit as follows — the

th entry of this is the dot product

sgn

(

) with (

)

−1

. This is the coefficient vector we would obtain

if we tried to regress

. If this is large, then this suggests we expect

to look correlated to

, and so it would be difficult to determine if

is part of

S or not.

Theorem.

(i) If k∆k

∞

≤ 1, or equivalently

max

k∈N

|sgn(β

)

−1

| ≤ 1,

and moreover

|β

| > λ



sgn(β

)





−1



for all

k ∈ S

, then there exists a Lasso solution

with

sgn

(

) =

sgn

(

(ii) If there exists a Lasso solution with sgn(

) = sgn(β

), then k∆k

∞

≤ 1.

We are going to make heavy use of the KKT conditions.

Proof.

Write

, and write

= 0

}

. Then the KKT conditions

are that

(β

−

β) = λˆν,

where kˆνk

∞

≤ 1 and ˆν

= sgn(

We expand this to say





−



= λ



ˆν



Call the top and bottom equations (1) and (2) respectively.

It is easier to prove (ii) first. If there is such a solution, then

= 0. So

from (1), we must have

(β

−

) = λˆν

Inverting

, we learn that

−

= λ





−1

sgn(β

Substitute this into (2) to get





−1

sgn(β

) = λˆν

By the KKT conditions, we know kˆν

∞

≤ 1, and the LHS is exactly λ∆.

To prove (1), we need to exhibit a

that agrees in sign with

and satisfies

the equations (1) and (2). In particular,

= 0. So we try

(

, ˆν

) =

− λ





−1

sgn(β

), sgn(β

)

(

, ν

) = (0, ∆).

This is by construction a solution. We then only need to check that

sgn(

) = sgn(β

which follows from the second condition.

Prediction and estimation

We now consider other question of how good the Lasso functions as a regression

method. Consider the model

Y = Xβ

+ ε − ¯ε1,

where the

are independent and have common sub-Gaussian parameter

. Let

S, s, N be as before.

As before, the Lasso doesn’t always behave well, and whether or not it does

is controlled by the compatibility factor.

Definition (Compatibility factor). Define the compatibility factor to be

= inf

β∈R

kβ

≤3kβ

6=0

kXβk

kβ

= inf

β∈R

kβ

k=1

kβ

≤3

− X

Note that we are free to use a minus sign inside the norm since the problem

is symmetric in β

↔ −β

In some sense, this

measures how well we can approximate

just with

the noise variables.

Definition (Compatibility condition). The compatibility condition is φ

> 0.

Note that if

has minimum eigenvalue

min

0, then we have

≥ c

min

. Indeed,

kβ

= sgn(β

)

≤

√

skβ

≤

√

skβk

and so

≥ inf

β6=0

kXβk

kβk

= c

min

Of course, we don’t expect the minimum eigenvalue to be positive, but we have

the restriction in infimum in the definition of

and we can hope to have a

positive value of φ

Theorem. Assume φ

> 0, and let

β be the Lasso solution with

λ = Aσ

log p/n.

Then with probability at least 1 − 2p

−(A

/8−1)

, we have

kX(β

−

β)k

+ λk

β − β

≤

16λ

16A

log p

sσ

This is actually two bounds. This simultaneously bounds the error in the

fitted values, and a bound on the error in predicting

β − β

Recall that in our previous bound, we had a bound of

∼

√

, and now we have

∼

. Note also that

sσ

is the error we would get if we magically knew which

were the non-zero variables and did ordinary least squares on those variables.

This also tells us that if

has a component that is large, then

must be

non-zero in that component as well. So while the Lasso cannot predict exactly

which variables are non-zero, we can at least get the important ones.

Proof. Start with the basic inequality Q

(

β) ≤ Q

(β

), which gives us

kX(β

−

β)k

+ λk

βk

≤

β − β

) + λkβ

We work on the event

Ω =



εk

∞

≤



where after applying H¨older’s inequality, we get

kX(β

−

β)k

+ 2λk

βk

≤ λk

β − β

+ 2λkβ

We can move 2

λk

βk

to the other side, and applying the triangle inequality, we

have

kX(

β − β

≤ 3λk

β − β

If we manage to bound the RHS from above as well, so that

3λk

β − β

k ≤ cλ

√

kX(

β − β

for some c, then we obtain the bound

kX(β − β

≤ c

Plugging this back into the second bound, we also have

3λk

β − β

≤ c

To obtain these bounds, we want to apply the definition of

β −β

. We thus

need to show that the

β − β

satisfies the conditions required in the infimum

taken.

Write

a =

nλ

kX(

β − β

Then we have

a + 2(k

+ k

) ≤ k

− β

+ k

+ 2kβ

Simplifying, we obtain

a + k

≤ k

− β

+ 2kβ

− 2k

Using the triangle inequality, we write this as

a + k

− β

≤ 3k

− β

So we immediately know we can apply the compatibility condition, which gives

≤

kX(

β − β

− β

Also, we have

kX(

β − β

+ λk

β − β

≤ 4λk

− β

Thus, using the compatibility condition, we have

kX(

β − β

+ λk

β − β

k ≤

4λ

kX(

β − β

Thus, dividing through by

√

kX(

β − β

, we obtain

√

kX(

β − β

≤

4λ

√

. (∗)

So we substitute into the RHS (∗) and obtain

kX(

β − β

+ λk

β − β

≤

16λ

If we want to be impressed by this result, then we should make sure that

the compatibility condition is a reasonable assumption to make on the design

matrix.

The compatibility condition and random design

For any Σ im R

p×p

, define

(S) = inf

β:kβ

≤3kβ

, β

6=0

Σβ

kβ

/|S|

Our original φ

is then the same as φ

(S).

If we want to analyze how

(

) behaves for a “random” Σ, then it would

be convenient to know that this depends continuously on Σ. For our purposes,

the following lemma suffices:

Lemma. Let Θ, Σ ∈ R

p×p

. Suppose φ

(S) > 0 and

max

j,k

|Θ

− Σ

| ≤

(S)

32|S|

Then

(S) ≥

(S).

Proof.

We suppress the dependence on

for notational convenience. Let

|S|

and

t =

(S)

32s

We have

|β

(Σ − Θ)β| ≤ kβk

k(Σ − Θ)βk

∞

≤ tkβk

where we applied H¨older twice.

If kβ

k ≤ 3kβ

, then we have

kβk

≤ 4kβ

≤ 4

Θβ

√

Thus, we have

Θβ −

32s

16β

Θβ

Θβ ≤ β

Σβ.

Define

Σ,s

= min

S:|S|=s

(S).

Note that if

max

|Θ

− Σ

| ≤

Θ,s

32s

then

(S) ≥

(S).

for all S with |S| = s. In particular,

Σ,s

≥

Θ,s

Theorem. Suppose the rows of

are iid and each entry is sub-Gaussian with

parameter

. Suppose

log p/n →

0 as

n → ∞

, and

is bounded away

from 0. Then if Σ

= E

Σ, then we have



Σ,s

≥



→ 1 as n → ∞.

Proof. It is enough to show that

max

− Σ

| ≤

32s

→ 0

as n → ∞.

Let t =

32s

. Then



max

j,k

− Σ

| ≥ t



≤

j,k

P(|

− Σ

| ≥ t).

Recall that

i=1

So we can apply Bernstein’s inequality to bound

P(|

− Σ

) ≤ 2 exp



−

2(64v

+ 4v



since σ = 8v

and b = 4v

. So we can bound



max

j,k

− Σ

| ≥ t



≤ 2p

exp



−



= 2 exp



−



c −

n log p



→ 0

for some constant c.

Corollary. Suppose the rows of

are iid mean-zero multivariate Gaussian

with variance Σ

. Suppose Σ

has minimum eigenvalue bounded from below by

min

0, and suppose the diagonal entries of Σ

are bounded from above. If

log p/n → 0, then



Σ,s

≥

min



→ 1 as n → ∞.

3.6 Computation of Lasso solutions

We have had enough of bounding things. In this section, let’s think about how

we can actually run the Lasso. What we present here is actually a rather general

method to find the minimum of a function, known as coordinate descent.

Suppose we have a function

→ R

. We start with an initial guess

(0)

and repeat for m = 1, 2, . . .

(m)

= argmin

f(x

, x

(m−1)

, . . . , x

(m−1)

)

(m)

= argmin

f(x

(m)

, x

(m−1)

, . . . , x

(m−1)

)

(m)

= argmin

f(x

(m)

, x

(m)

, . . . , x

(m)

d−1

, x

)

until the result stabilizes.

This was proposed for solving the Lasso a long time ago, and a Stanford

group tried this out. However, instead of using

(m)

when computing

(m)

, they

used

(m−1)

instead. This turned out to be pretty useless, and so the group

abandoned the method. After trying some other methods, which weren’t very

good, they decided to revisit this method and fixed their mistake.

For this to work well, of course the coordinatewise minimizations have to be

easy (which is the case for the Lasso, where we even have explicit solutions). This

converges to the global minimizer if the minimizer is unique,

(

)

≤ f

(

(0)

)

}

is compact, and if f has the form

f(x) = g(x) +

where

is convex and differentiable, and each

is convex but not necessarily

differentiable. In the case of the Lasso, the first is the least squared term, and

the h

is the `

term.

There are two things we can do to make this faster. We typically solve the

Lasso on a grid of

values

> λ

> ··· > λ

, and then picking the appropriate

-fold cross-validation. In this case, we can start solving at

, and then

for each

i >

0, we solve the

problem by picking

(0)

to be the optimal

solution to the

i−1

problem. In fact, even if we already have a fixed

value we

want to use, it is often advantageous to solve the Lasso with a larger

-value,

and then use that as a warm start to get to our desired λ value.

Another strategy is an active set strategy. If

is large, then this loop may

take a very long time. Since we know the Lasso should set a lot of things to zero,

for ` = 1, . . . , L, we set

A = {k :

`−1

6= 0}.

We then perform coordinate descent only on coordinates in

to obtain a Lasso

solution

with

= 0. This may not be the actual Lasso solution. To check

this, we use the KKT conditions. We set

V =



k ∈ A

(Y − X

β)| > λ



∅

, and we are done. Otherwise, we add

to our active set

, and then

run coordinate descent again on this active set.

3.7 Extensions of the Lasso

There are many ways we can modify the Lasso. The first thing we might want to

change in the Lasso is to replace the least squares loss with other log-likelihoods.

Another way to modify the Lasso is to replace the

penalty with something

else in order to encourage a different form of sparsity.

Example (Group Lasso). Given a partition

∪ ···G

= {1, . . . , p},

the group Lasso penalty is

j=1

kβ

where

}

is some sort of weight to account for the fact that the groups have

different sizes. Typically, we take m

If we take

{i}

, then we recover the original Lasso. If we take

= 1,

then we recover Ridge regression. What this does is that it encourages the entire

group to be all zero, or all non-zero.

Example. Another variation is the fused Lasso. If

j+1

is expected to be close

to β

, then a fused Lasso penalty may be appropriate, with

p−1

j=1

|β

j+1

− β

| + λ

kβk

For example, if

= µ

+ ε

and we believe that (

)

i=1

form a piecewise constant sequence, we can estimate

argmin

(

kY − µk

+ λ

n−1

i=1

|µ

i+1

− µ

)

Example (Elastic net). We can use

λ,α

= argmin



kY − Xβk

+ λ(αkβk

+ (1 − α)kβk

)



for

α ∈

1]. What the

norm does is that it encourages highly positively

correlated variables to have similar estimated coefficients.

For example, if we have duplicate columns, then the

penalty encourages

us to take one of the coefficients to be 0, while the

penalty encourages the

coefficients to be the same.

Another class of variations try to reduce the bias of the Lasso. Although

the bias of the Lasso is a necessary by-product of reducing the variance of the

estimate, it is sometimes desirable to reduce this bias.

The LARS-OLS hybrid takes the

obtained by the Lasso, and then re-

estimate

by OLS. We can also re-estimate using the Lasso on

, and this

is known as the relaxed Lasso.

In the adaptive Lasso, we obtain an initial estimate of

, e.g. with the Lasso,

and then solve

adapt

= argmin

β:β







kY − Xβk

+ λ

k∈

|β







We can also try to use a non-convex penalty. We can attempt to solve

argmin

(

kY − Xβk

k=1

(|β

)

where

: [0

, ∞

)

→ p

, ∞

) is a non-convex function. One common example is

the MCP , given by

(u) =



λ −



where

is an extra tuning parameter. This tends to give estimates even sparser

than the Lasso.

4 Graphical modelling

4.1 Conditional independence graphs

So far, we have been looking at prediction problems. But sometimes we may

want to know more than that. For example, there is a positive correlation

between the wine budget of a college, and the percentage of students getting

firsts. This information allows us to make predictions, in the sense that if we

happen to know the wine budget of a college, but forgot the percentage of

students getting firsts, then we can make a reasonable prediction of the latter

based on the former. However, this does not suggest any causal relation between

the two — increasing the wine budget is probably not a good way to increase

the percentage of students getting firsts!

Of course, it is unlikely that we can actually figure out causation just by

looking at the data. However, there are some things we can try to answer. If

we gather more data about the colleges, then we would probably find that the

colleges that have larger wine budget and more students getting firsts are also

the colleges with larger endowment and longer history. If we condition on all

these other variables, then there is not much correlation left between the wine

budget and the percentage of students getting firsts. This is what we are trying

to capture in conditional independence graphs.

We first introduce some basic graph terminology. For the purpose of condi-

tional independence graphs, we only need undirected graphs. But later, we need

the notion of direct graphs as well, so our definitions will be general enough to

include those.

Definition (Graph). A graph is a pair

= (

V, E

), where

is a set and

E ⊆ (V, V ) such that (v, v) 6∈ E for all v ∈ V .

Definition (Edge). We say there is an edge between

and

and that

and

are adjacent if (j, k) ∈ E or (k, j) ∈ E.

Definition (Undirected edge). An edge (

j, k

) is undirected if also (

k, j

)

∈ E

Otherwise, it is directed and we write

j → k

to represent it. We also say that

is a parent of k, and write pa(k) for the set of all parents of k.

Definition ((Un)directed graph). A graph is (un)directed if all its edges are

(un)directed.

Definition (Skeleton). The skeleton of

is a copy of

with every edge replaced

by an undirected edge.

Definition (Subgraph). A graph

= (

V, E

) is a subgraph of

= (

V, E

) if

⊆ V

and

⊆ E

. A proper subgraph is one where either of the inclusions are

proper inclusions.

As discussed, we want a graph that encodes the conditional dependence of

different variables. We first define what this means. In this section, we only

work with undirected graphs.

Definition (Conditional independence). Let

X, Y, Z

be random vectors with

joint density

XY Z

. We say that

is conditionally independent of

given

written X q Y | Z, if

XY |Z

(x, y | z) = f

X|Z

(x | z)f

Y |Z

(y | z).

Equivalently,

X|Y Z

(x | y, z) = f

X|Z

(x | z)

for all y.

We shall ignore all the technical difficulties, and take as an assumption that

all these conditional densities exist.

Definition (Conditional independence graph (CIG)). Let

be the law of

= (

, . . . , Z

)

. The conditional independent graph (CIG) is the graph whose

vertices are

{

, . . . , p}

, and contains an edge between

and

iff

and

are

conditionally dependent given Z

−jk

More generally, we make the following definition:

Definition (Pairwise Markov property). Let

be the law of

= (

, . . . , Z

)

We say

satisfies the pairwise Markov property with respect to a graph

if for

any distinct, non-adjacent vertices j, k, we have Z

q Z

| Z

−jk

Example. If

is a complete graph, then

satisfies the pairwise Markov

property with respect to G.

The conditional independence graph is thus the minimal graph satisfying the

pairwise Markov property. It turns out that under mild conditions, the pairwise

Markov property implies something much stronger.

Definition (Separates). Given a triple of (disjoint) subsets of nodes

A, B, S

we say

separates

from

if every path from a node in

to a node in

contains a node in S.

Definition (Global Markov property). We say

satisfies the global Markov

property with respect to

if for any triple of disjoint subsets of

(

A, B, S

), if

S separates A and B, then Z

q Z

| Z

Proposition. If

has a positive density, then if it satisfies the pairwise Markov

property with respect to G, then it also satisfies the global Markov property.

This is a really nice result, but we will not prove this. However, we will prove

a special case in the example sheet.

So how do we actually construct the conditional independence graph? To do

so, we need to test our variables for conditional dependence. In general, this is

quite hard. However, in the case where we have Gaussian data, it is much easier,

since independence is the same as vanishing covariance.

Notation (

A,B

). Let

be a matrix. Then

A,B

refers to the submatrix

given by the rows in A and columns in B.

Since we are going to talk about conditional distributions a lot, the following

calculation will be extremely useful.

Proposition. Suppose Z ∼ N

(µ, Σ) and Σ is positive definite. Then

| Z

= z

∼ N

|A|

(µ

+ Σ

A,B

−1

B,B

− µ

), Σ

A,A

− Σ

A,B

−1

B,B

B,A

Proof.

Of course, we can just compute this directly, maybe using moment

generating functions. But for pleasantness, we adopt a different approach. Note

that for any M , we have

= MZ

+ (Z

− MZ

We shall pick M such that Z

− MZ

is independent of Z

, i.e. such that

0 = cov(Z

, Z

− MZ

) = Σ

B,A

− Σ

B,B

So we should take

M = (Σ

−1

B,B

B,A

)

= Σ

A,B

−1

B,B

We already know that

−MZ

is Gaussian, so to understand it, we only need

to know its mean and variance. We have

E[Z

− MZ

] = µ

− Mµ

= µ

− Σ

−1

var(Z

− MZ

) = Σ

A,A

− 2Σ

A,B

−1

B,B

B,A

+ Σ

A,B

−1

B,B

−1

B,B

B,A

= Σ

A,A

− Σ

A,B

−1

B,B

B,A

Then we are done.

Neighbourhood selection

We are going to specialize to

{k}

and

{

, . . . , n} \ {k}

. Then we can

separate out the “mean” part and write

= M

+ Z

−k

−1

−k,−k

−k,k

+ ε

where

= µ

− µ

−k

−1

k,−k

−k,k

| Z

−k

∼ N(0, Σ

k,k

− Σ

k,−k

−1

−k,−k

−k,k

This now looks like we are doing regression.

We observe that

Lemma. Given

, let

be such that (

−k

)

. This

is either

+ 1,

depending on whether it comes after or before k.

If the jth component of Σ

−1

−k,−k

−k,k

is 0, then Z

q Z

| Z

−kj

Proof.

If the

th component of Σ

−1

−k,−k

−k,k

is 0, then the distribution of

| Z

−k

will not depend on (

−k

)

(here

is either

+ 1, depending

on where k is). So we know

| Z

−k

= Z

| Z

−kj

This is the same as saying Z

q Z

| Z

−kj

Neighbourhood selection exploits this fact. Given

, . . . , x

which are iid

∼ Z and

X = (x

, ··· , x

)

we can estimate Σ

−1

−k,−k

−k,k

by regressing

−k

using the Lasso (with

an intercept term). We then obtain selected sets

. There are two ways of

estimating the CIG based on these:

– OR rule: We add the edge (j, k) if j ∈

or k ∈

– AND rule: We add the edge (j, k) if j ∈

and k ∈

The graphical Lasso

Another way of finding the conditional independence graph is to compute

var(Z

| Z

−jk

) directly. The following lemma will be useful:

Lemma. Let M ∈ R

p×p

be positive definite, and write

M =



P Q



where P and Q are square. The Schur complement of R is

S = P − QR

−1

Note that this has the same size as P . Then

(i) S is positive definite.

(ii)

−1



−1

−S

−1

−R

−1

+ R

−1



(iii) det(M) = det(S) det(R)

We have seen this Schur complement when we looked at

var

(

| Z

)

previously, where we got

var(Z

| Z

) = Σ

A,A

− Σ

A,A

−1

= Ω

−1

A,A

where Ω = Σ

−1

is the precision matrix.

Specializing to the case where A = {j, k}, we have

var(Z

{j,k}

| Z

−jk

) =

det(Ω

A,A

)



Ω

k,k

−Ω

j,k

−Ω

j,k

Ω

j,j



This tells us

| Z

−kj

iff Ω

= 0. Thus, we can approximate the conditional

independence graph by computing the precision matrix Ω.

Our method to estimate the precision matrix is similar to the Lasso. Recall

that the density of N

(µ, Σ) is

P (z) =

(2π)

p/2

(det Σ)

1/2

exp



−

(z − µ)

−1

(z − µ)



The log-likelihood of (

µ,

Ω) based on an iid sample (

, . . . , X

) is (after dropping

a constant)

`(µ, Ω) =

log det Ω −

i=1

− µ)

Ω(x

− µ).

To simplify this, we let

¯x =

i=1

, S =

i=1

− ¯x)(x

− ¯x)

Then

− µ)

Ω(x

− µ) =

− ¯x + ¯x − µ)

Ω(x

− ¯x + ¯x − µ)

− ¯x)

Ω(x

− ¯x) + n(

X − µ)

Ω(

X − µ).

We have

− ¯x)

Ω(x

− ¯x) =



− ¯x)

Ω(x

− ¯x)



= n tr(SΩ).

So we now have

`(µ, Ω) = −



tr(SΩ) − log det Ω + (

X − µ)

Ω(

X − µ)



We are really interested in estimating Ω. So we should try to maximize this over

, but that is easy, since this is the same as minimizing (

X − µ

)

Ω(

X − µ

), and

we know Ω is positive-definite. So we should set µ =

X. Thus, we have

`(Ω) = max

µ∈R

`(µ, Ω) = −



tr(SΩ) − log det ω



So we can solve for the MLE of Ω by solving

min

Ω:Ω0



− log det Ω + tr(SΩ)



One can show that this is convex, and to find the MLE, we can just differentiate

∂

∂Ω

log det Ω = (Ω

−1

)

∂

∂Ω

tr(SΩ) = S

using that

and Ω are symmetric. So provided that

is positive definite, the

maximum likelihood estimate is just

Ω = S

−1

But we are interested in the high dimensional situation, where we have loads of

variables, and

cannot be positive definite. To solve this, we use the graphical

Lasso.

The graphical Lasso solves

argmin

Ω:Ω0



− log det Ω + tr(SΩ) + λkΩk



where

kΩk

Ω

Often, people don’t sum over the diagonal elements, as we want to know if off-

diagonal elements ought to be zero. Similar to the case of the Lasso, this gives a

sparse estimate of Ω from which we may estimate the conditional independence

graph.

4.2 Structural equation modelling

The conditional independence graph only tells us which variables are related to

one another. However, it doesn’t tell us any causal relation between the different

variables. We first need to explain what we mean by a causal model. For this,

we need the notion of a directed acyclic graph.

Definition (Path). A path from

is a sequence

, j

, . . . , j

(at least two) distinct vertices such that j

and j

`+1

are adjacent.

A path is directed if j

→ j

`+1

for all `.

Definition (Directed acyclic graph (DAG)). A directed cycle is (almost) a

directed path but with the start and end points the same.

A directed acyclic graph (DAG) is a directed graph containing no directed

cycles.

DAG

not DAG

We will use directed acyclic graphs to encode causal structures, where we

have a path from a to b if a “affects” b.

Definition (Structural equation model (SEM)). A structural equation model

for a random vector Z ∈ R

is a collection of equations

= h

, ε

where

= 1

, . . . , p

and

, . . . , ε

are independent, and

⊆ {

, . . . , p}\{k}

and

such that the graph with pa(k) = p

is a directed acyclic graph.

Example. Consider three random variables:

– Z

= 1 if a student is taking a course, 0 otherwise

– Z

= 1 if a student is attending catch up lectures, 0 otherwise

– Z

= 1 if a student heard about machine learning before attending the

course, 0 otherwise.

Suppose

= ε

∼ Bernoulli(0.25)

= 1

{ε

(1+Z

}

, ε

∼ U[0, 1]

= 1

{ε

}

, ε

∼ U[0, 1].

This is then an example of a structural equation modelling, and the corresponding

DAG is

Note that the structural equation model for

determines its distribution,

but the converse is not true. For example, the following two distinct structural

equation give rise to the same distribution:

= ε Z

= Z

= ε

Indeed, if we have two variables that are just always the same, it is hard to tell

which is the cause and which is the effect.

It would be convenient if we could order our variables in a way that

depends only on Z

for j < k. This is known as a topological ordering:

Definition (Descendant). We say

is a descendant of

if there is a directed

path from j to k. The set of descendant of j will be denoted de(j).

Definition (Topological ordering). Given a DAG

with

{

, . . . , p}

we say

that a permutation

V → V

is a topological ordering if

(

)

< π

(

) whenever

k ∈ de(j).

Thus, given a topological ordering

, we can write

as a function of

−1

(1)

, . . . , ε

−1

(π(k))

How do we understand structural equational models? They give us informa-

tion that are not encoded in the distribution itself. One way to think about them

is via interventions. We can modify a structural equation model by replacing

the equation for

by setting, e.g.

. In real life, this may correspond

to forcing all students to go to catch up workshops. This is called a perfect

intervention. The modified SEM gives a new joint distribution for

. Expecta-

tions or probabilities with respect to the new distribution are written by adding

“do(Z

= a)”. For example, we write

E(Z

| do(Z

= a)).

In general, this is different to

(

| Z

), since, for example, if we conditioned

on Z

= a in our example, then that would tell us something about Z

Example. After the intervention

(

= 1), i.e. we force everyone to go to the

catch-up lectures, we have a new SEM with

= ε

∼ Bernoulli(0.25)

= 1

(1+Z

, ε

∼ U[0, 1].

Then, for example, we can compute

P(Z

= 1 | do(Z

= 1)) =

and by high school probability, we also have

P(Z

= 1 | Z

= 1) =

To understand the DAGs associated to structural equation models, we would

like to come up with Markov properties analogous to what we had for undirected

graphs. This, in particular, requires the correct notion of “separation”, which we

call d-separation. Our notion should be such that if

d-separates

and

the DAG, then

and

are conditionally independent given

. Let’s think

about some examples. For example, we might have a DAG that looks like this:

Then we expect that

(i) Z

and Z

are not independent;

(ii) Z

and Z

are independent given Z

;

(iii) Z

and Z

are independent;

(iv) Z

and Z

are not independent given Z

We can explain a bit more about the last one. For example, the structural

equation model might telling us

. In this case, if we know that

is large but

is small, then chances are,

is also large (in the opposite

sign). The point is that both

and

both contribute to

, and if we know

one of the contributions and the result, we can say something about the other

contribution as well.

Similarly, if we have a DAG that looks like

then as above, we know that Z

and Z

are not independent given Z

Another example is

Here we expect

– Z

and Z

are not independent.

– Z

and Z

are independent given Z

To see (i), we observe that if we know about

, then this allows us to predict

some information about

, which would in turn let us say something about

Definition (Blocked). In a DAG, we say a path (

, . . . , j

) between

and

is blocked by a set of nodes

(with neither

nor

) if there is some

∈ S and the path is not of the form

`−1

`+1

or there is some

such that the path is of this form, but neither

nor any of

its descendants are in S.

Definition (d-separate). If

is a DAG, given a triple of (disjoint) subsets of

nodes

A, B, S

, we say

d-separates

from

blocks every path from

For convenience, we define

Definition (

-structure). A set of three nodes is called a

-structure if one node

is a child of the two other nodes, and these two nodes are not adjacent.

It is now natural to define

Definition (Markov properties). Let

be the distribution of

and let

the density. Given a DAG G, we say P satisfies

(i) the Markov factorization criterion if

f(z

, . . . , z

) =

k=1

f(z

| z

pa(k)

(ii)

the global Markov property if for all disjoint

A, B, S

such that

A, B

d-separated by S, then Z

q Z

| Z

Proposition. If

has a density with respect to a product measure, then (i)

and (ii) are equivalent.

How does this fit into the structural equation model?

Proposition. Let

be the structural equation model with DAG

. Then

obeys the Markov factorization property.

Proof.

We assume

is topologically ordered (i.e. the identity map is a topological

ordering). Then we can always write

f(z

, . . . , z

) = f(z

)f(z

| z

) ···z(z

| z

, z

, . . . , z

p−1

By definition of a topological order, we know

(

)

⊆ {

, . . . , k −

}

. Since

is a function of Z

pa(k)

and independent noise ε

. So we know

q Z

{1,...,p}\{k∪pa(k)}

| Z

pa(k)

Thus,

f(z

| z

, . . . , z

k−1

) = f(z

| z

pa(k)

4.3 The PC algorithm

We now want to try to find out the structural equation model given some data,

and in particular, determine the causal structure. As we previously saw, there is

no hope of determining this completely, even if we know the distribution of the

Z completely. Let’s consider the different obstacles to this problem.

Causal minimality

is generated by an SEM with DAG

, then from the above, we know that

is Markov with respect to

. The converse is also true: if

is Markov with

respect to a DAG

, then there exists a SEM with DAG

that generates

This immediately implies that

will be Markov with respect to many DAGs.

For example, a DAG whose skeleton is complete will always work. This suggests

the following definition:

Definition (Causal minimality). A distribution

satisfies causal minimality

with respect to G but not any proper subgraph of G.

Markov equivalent DAGs

It is natural to aim for finding a causally minimal DAG. However, this does

not give a unique solution, as we saw previously with the two variables that are

always the same.

In general, two different DAGs may satisfy the same set of d-separations,

and then a distribution is Markov with respect to one iff its Markov with respect

to the other, and we cannot distinguish between the two.

Definition (Markov equivalence). For a DAG G, we let

M(G) = {distributions P such that P is Markov with respect to G}.

We say two DAGs G

, G

are are Markov equivalent if M(G

) = M(G

What is nice is that there is a rather easy way of determining when two

DAGs are Markov equivalent.

Proposition. Two DAGs are Markov equivalent iff they have the same skeleton

and same set of v-structure.

The set of all DAGs that are Markov equivalent to a given DAG can be

represented by a CPDAG (completed partial DAG), which contains an edge

(j, k) iff some member of the equivalence class does.

Faithfulness

To describe the final issue, consider the SEM

= ε

, Z

= αZ

+ ε

, Z

= βZ

+ γZ

+ ε

We take ε ∼ N

(0, I). Then we have Z = (Z

, Z

) ∼ N(0, Σ), where

Σ =





1 α β + αγ

α α

+ 1 α(β + αγ) + γ

β + αγ α(β + αγ) + γ β

+ γ

(α

+ 1) + 1 + 2αβγ





Now if we picked values of α, β, γ such that

β + αγ = 0,

then we obtained an extra independence relation

q Z

in our system. For

example, if we pick β = −1 and α, γ = 1, then

Σ =





1 1 0

1 2 1

0 1 2





While there is an extra independence relation, we cannot remove any edge while

still satisfying the Markov property. Indeed:

–

If we remove 1

→

2, then this would require

q Z

, but this is not true.

– If we remove 2 → 3, then this would require Z

q Z

| Z

, but we have

var((Z

, Z

) | Z

) =



2 1

1 2



−







1 0





1 1

1 2



and this is not diagonal.

– If we remove 1 → 3, then this would require Z

q Z

| Z

, but

var((Z

, Z

) | Z

) =



1 0

0 2



−







1 1



which is again non-diagonal.

So this DAG satisfies causal minimality. However,

can also be generated by

the structural equation model

= ˜ε

+ ˜ε

= ˜ε

where the ˜ε

are independent with

˜ε

∼ N(0, 1), ˜ε

∼ N(0, 2), ˜ε

∼ N(0,

Then this has the DAG

This is a strictly smaller DAG in terms of the number of edges involved. It is

easy to see that this satisfies causal minimality.

Definition (Faithfulness). We say

is faithful to a DAG

if it is Markov

with respect to

and for all

A, B, S

disjoint,

q Z

| Z

implies

A, B

are

d-separated by S.

Determining the DAG

We shall assume our distribution is faithful to some

, and see if we can figure

out G

from P , or even better, from data.

To find G, the following proposition helps us compute the skeleton:

Proposition. If nodes

and

are adjacent in a DAG

, then no set can

d-separate them.

If they are not adjacent, and

is a topological order for

with

(

)

< π

(

then they are d-separated by pa(k).

Proof.

Only the last part requires proof. Consider a path

, . . . , j

Start reading the path from

and go backwards. If it starts as

m−1

→ k

, then

m−1

is a parent of k and blocks the path. Otherwise, it looks like k → j

m−1

We keep going down along the path until we first see something of the form

···

Thus must exist, since j is not a descendant of k by topological ordering. So it

suffices to show that

does not have a descendant in

(

), but if it did, then

this would form a closed loop.

Finding the

-structures is harder, and at best we can do so up to Markov

equivalence. To do that, observe the following:

Proposition. Suppose we have j − ` − k in the skeleton of a DAG.

(i) If j → ` ← k, then no S that d-separates j can have ` ∈ S.

(ii)

If there exists

that d-separates

and

` 6∈ S

, then

j → ` ← k

Denote the set of nodes adjacent to the vertex

in the graph

adj

(

G, k

We can now describe the first part of the PC algorithm, which finds the

skeleton of the “true DAG”:

(i) Set

G to be the complete undirected graph. Set ` = −1.

(ii) Repeat the following steps:

(a) Set ` = ` + 1:

(b) Repeat the following steps:

Select a new ordered pair of nodes

j, k

that are adjacent in

and

such that |adj(

G, j) \{k}| ≥ `.

ii. Repeat the following steps:

A. Choose a new S ⊆ adj(

G, j) \{k} with |S| = `.

q Z

| Z

, then delete the edge

, and store

(

k, j

) =

S(j, k) = S

C. Repeat until j − k is deleted or all S chosen.

iii. Repeat until all pairs of adjacent nodes are inspected.

Suppose

is faithful to a DAG

. At each stage of the algorithm, the skeleton

will be a subgraph of

. On the other hand, edges (

j, k

) remaining at

termination will have

q Z

| Z

for all S ⊆ (

G, k), S ⊆ (

G, j).

So they must be adjacent in G

. Thus,

G and G

have the same skeleton.

To find the v-structures, we perform:

(i) For all j − l − k in

G, do:

(a) If ` 6∈ S(j, k), then orient j → ` ← k.

This gives us the Markov equivalence class, and we may orient the other edges

using other properties like acyclicity.

If we want to apply this to data sets, then we need to apply some conditional

independence tests instead of querying our oracle to figure out if things are

conditional dependence. However, errors in the algorithm propagate, and the

whole process may be derailed by early errors. Moreover, the result of the

algorithm may depend on how we iterate through the nodes. People have tried

many ways to fix these problems, but in general, this method is rather unstable.

Yet, if we have large data sets, this can produce quite decent results.

5 High-dimensional inference

5.1 Multiple testing

Finally, we talk about high-dimensional inference. Suppose we have come up

with a large number of potential drugs, and want to see if they are effective in

killing bacteria. Naively, we might try to run a hypothesis test on each of them,

using a

p <

05 threshold. But this is a terrible idea, since each test has a 0

chance of giving a false positive, so even if all the drugs are actually useless, we

would have incorrectly believed that a lot of them are useful, which is not the

case.

In general, suppose we have some null hypothesis

, . . . , H

. By definition,

a p value p

for H

is a random variable such that

≤ α) ≤ α

for all α ∈ [0, 1].

Let

be the number of true null hypothesis. Given a procedure for

rejecting hypothesis (a multiple testing procedure), we let

be the number of

false rejections (false positives), and

the total number of rejections. One can

also think about the number of false negatives, but we shall not do that here.

Traditionally, multiple-testing procedures sought to control the family-wise

error rate (FWER), defined by

(

N ≥

1). The simplest way to minimize this is

to use the Bonferroni correction, which rejects

≤

. Usually, we might

have

α ∼

05, and so this would be very small if we have lots of hypothesis (e.g.

a million). Unsurprisingly, we have

Theorem. When using the Bonferroni correction, we have

FWER ≤ E(N) ≤

≤ α.

Proof.

The first inequality is Markov’s inequality, and the last is obvious. The

second follows from

E(N) = E

i∈I

≤α/m

i∈I



≤



≤

The Bonferroni is a rather conservative procedure, since all these inequalities

can be quite loose. When we have a large number of hypotheses, the criterion

for rejection is very very strict. Can we do better?

A more sophisticated approach is the closed testing procedure. For each

non-empty subset

I ⊆ {

, . . . , m}

, we let

be the null hypothesis that

true for all

i ∈ I

. This is known as an intersection hypothesis. Suppose for each

I ⊆ {

, . . . , m}

non-empty, we have an

-level test

for

(a local test), so

that

(φ

= 1) ≤ α.

Here Φ

takes values in

{

}

, and

= 1 means rejection. The closed testing

procedure then rejects H

iff for all J ⊇ I, we have φ

= 1.

Example. Consider the tests, where the red ones are the rejected one:

1234

134

124

134

234

In this case, we reject

but not

by closed testing. While

is rejected,

we cannot tell if it is H

or H

that should be rejected.

This might seem like a very difficult procedure to analyze, but it turns out it

is extremely simple.

Theorem. Closed testing makes no false rejections with probability

≥

− α

In particular, FWER ≤ α.

Proof.

In order for there to be a false rejection, we must have falsely rejected

with the local test, which has probability at most α.

But of course this doesn’t immediately give us an algorithm we can apply to

data. Different choices for the local test give rise to different multiple testing

procedures. One famous example is Holm’s procedure. This takes

to be the

Bonferroni test, where φ

= 1 iff p

≤

|I|

for some i ∈ I.

When

is large, then we don’t want to compute all

, since there are 2

computations to do. So we might want to find a shortcut. With a moment of

thought, we see that Holm’s procedure amounts to the following:

–

Let (

) be the index of the

th smallest

-value, so

(1)

≤ p

(2)

≤ ··· ≤ p

(m)

–

Step 1: If

(1)

≤

, then reject

(1)

and go to step 2. Otherwise, accept

all null hypothesis.

–

Step

: If

(i)

≤

m−i+1

, then reject

(i)

and go to step

+ 1. Otherwise,

accept H

(i)

, H

(i+1)

, . . . , H

(m)

– Step m: If p

(m)

≤ α, then reject H

(m)

. Otherwise, accept H

(m)

The interesting thing about this is that it has the same bound on FWER as the

Bonferroni correction, but the conditions here are less lenient.

But if

is very large, then the criterion for accepting

(1)

is still quite strict.

The problem is that controlling FWER is a very strong condition. Instead of

controlling the probability that there is one false rejection, when

is large, it

might be more reasonable to control the proportion of false discoveries. Many

modern multiple testing procedures aim to control the false discovery rate

FDR = E



max{R, 1}



The funny maximum in the denominator is just to avoid division by zero. When

= 0, then we must have

= 0 as well, so what is put in the denominator

doesn’t really matter.

The Benjamini–Hochberg procedure attempts to control the FDR at level

and works as follows:

–

Let

max



i : p

(i)

≤

αi



. Then reject

(1)

, . . . , H

(

, or accept all

hypothesis if

k is not defined.

Under certain conditions, this does control the false discovery rate.

Theorem. Suppose that for each

i ∈ I

is independent of

j 6

. Then

using the Benjamini–Hochberg procedure, the false discovery rate

F DR = E

max(R, 1)

≤

αM

≤ α.

Curiously, while the proof requires

to be independent of the others, in

simulations, even when there is no hope that the

are independent, it appears

that the Benjamini–Hochberg still works very well, and people are still trying to

understand what it is that makes Benjamini–Hochberg work so well in general.

Proof. The false discovery rate is

max(R, 1)

r=1

R=r

r=1

i∈I

≤αr/M

R=r

i∈I

r=1



≤

αr

, R = r



The brilliant idea is, for each

i ∈ I

, let

be the number of rejections when

applying a modified Benjamini–Hochberg procedure to

, . . . , p

}\{p

}

with cutoff

= max



j : p

(j)

≤

α(j + 1)



We observe that for i ∈ I

and r ≥ 1, we have

≤

αr

, R = r

≤

, p

(r)

≤

αr

, p

(s)

αs

for all s ≥ r

≤

αr

, p

(r−1)

≤

αr

, p

(s−1)

αs

for all s > r

≤

αr

, R

= r − 1

The key point is that

r −

1 depends only on the other

-values. So the

FDR is equal to

FDR =

i∈I

r=1



≤

αr

, R

= r − 1



i∈I

r=1



≤

αr



P(R

= r − 1)

Using that P(p

≤

αr

) ≤

αr

by definition, this is

≤

i∈I

r=1

P(R

= r − 1)

i∈I

P(R

∈ {0, . . . , m − 1})

αM

This is one of the most used procedures in modern statistics, especially in

the biological sciences.

5.2 Inference in high-dimensional regression

We have more-or-less some answer as to how to do hypothesis testing, given that

we know how to obtain these

-values. But how do we obtain these in the first

place?

For example, we might be trying to do regression, and are trying figure out

which coefficients are non-zero. The the low dimension setting, with the normal

linear model

Xβ

, where

ε ∼ N

, σ

). In the low-dimensional setting,

we have

√

(

OLS

− β

)

∼ N

, σ

(

)

−1

). Since this does not depend on

, we can use this to form confidence intervals and hypothesis tests.

However, if we have more coefficients than there are data points, then we

can’t do ordinary least squares. So we need to look for something else. For

example, we might want to replace the OLS estimate with the Lasso estimate.

However,

√

(

− β

) has an intractable distribution. In particular, since

has a bias, the distribution will depend on β

in a complicated way.

The recently introduced debiased Lasso tries to overcomes these issues. See

van de Geer, B¨uhlmann, Ritov, Dezeure (2014) for more details. Let

be the

Lasso solution at λ > 0. Recall the KKT conditions that says ˆν defined by

(Y − X

β) = λˆν

satisfies kˆνk

∞

≤ 1 and ˆν

= sgn(

), where

S = {k :

6= 0}.

We set

Σ =

X. Then we can rewrite the KKT conditions as

Σ(

β − β

) + λˆν =

ε.

What we are trying to understand is

β − β

. So it would be nice if we can find

some sort of inverse to

. If so, then

β − β

plus some correction term involving

ˆv would then be equal to a Gaussian.

Of course, the problem is that in the high dimensional setting, that

has

no hope of being invertible. So we want to find some approximate inverse

that the error we make is not too large. If we are equipped with such a

, then

we have

√

β + λ

Θˆν − β

) =

√

ΘX

ε + ∆,

where

∆ =

√

Σ − I)(β

−

β).

We hope we can choose

Θ so that δ is small. We can then use the quantity

b =

β + λ

Θˆν =

β +

ΘX

(Y − X

β)

as our modified estimator, called the debiased Lasso.

How do we bound ∆? We already know that (under compatibility and

sparsity conditions), we can make the

norm of

kβ

−

βk

small with high

probability. So if the

∞

norm of each of the rows of

Σ −

1 is small, then

H¨older allows us to bound ∆.

Write

for the jth row of

Θ. Then

)

− Ik

∞

≤ η

is equivalent to

(

)

| ≤ η

for

k 6

and

(

)

−

| ≤ η

. Using the

definition of

Σ, these are equivalent to

| ≤ η,



− 1



≤ η.

The first is the same as saying

−j

∞

≤ η.

This is quite reminiscent of the KKT conditions for the Lasso. So let us define

ˆγ

(j)

= argmin

γ∈R

p−1



− X

−j

γk

+ λ

kγk



ˆτ

− X

−j

ˆγ

(j)

) =

− X

−j

ˆγ

(j)

+ λ

kˆγ

(j)

The second equality is an exercise on the example sheet.

We can then set

= −

ˆτ

(ˆγ

(j)

, . . . , ˆγ

(j)

j−1

, −1, ˆγ

(j)

, . . . , ˆγ

(j)

p−1

)

The factor is there so that the second inequality holds.

Then by construction, we have

− X

−j

ˆγ

(j)

(X − X

−j

ˆγ

(j)

)/n

Thus, we have

= 1, and by the KKT conditions for the Lasso, we have

ˆτ

−j

∞

≤ λ

Thus, with the choice of

Θ above, we have

k∆k

∞

≤

√

β − β

max

ˆτ

Now this is good as long as we can ensure

ˆτ

to be small. When is this true?

We can consider a random design setting, where each row of

is iid

Σ)

for some positive-definite Σ. Write Ω = Σ

−1

Then by our study of the neighbourhood selection procedure, we know that

for each j, we can write

= X

−j

(j)

+ ε

(j)

where

(j)

| X

−j

∼ N

Ω

−1

) are iid and

(j)

−

Ω

−1

Ω

−j,j

. To apply our

results, we need to ensure that γ

(j)

are sparse. Let use therefore define

k6=j

Ω

6=0

and set s

max

= max(max

, s).

Theorem. Suppose the maximum eigenvalue of Σ is always at least

min

and

max

≤

1. Suppose further that

max

log(p)/n →

0. Then there exists

constants A

, A

such that setting λ = λ

= A

log(p)/n, we have

√

b − β

) = W + ∆

W | X ∼ N

(0, σ

and as n, p → ∞,



k∆k

∞

> A

log(p)

√



→ 0.

Note that here X is not centered and scaled.

We see that in particular,

√

(

− β

)

∼ N

, σ

(

)

). In fact, one

can show that

− X

−j

ˆγ

(j)

}

ˆτ

This suggests an approximate (1 − α)-level confidence interval for β

CI =



− Z

α/2

/n,

+ Z

α/2



where

is the upper

point of

1). Note that here we are getting confidence

intervals of width

∼

1/n

. In particular, there is no

log p

dependence, if we are

only trying to estimate β

only.

Proof. Consider the sequence of events Λ

defined by the following properties:

– φ

Σ,s

≥ c

min

/2 and φ

−j,−j

≥ c

min

/2 for all j

–

Σk

∞

≤ λ and

−j

(j)

∞

≤ λ.

–

(j)

≥ (Ω

)

−1

(1 − 4

(log p)/n)

Question 13 on example sheet 4 shows that

(Λ

)

→

1 for

sufficiently

large. So we will work on the event Λ

By our results on the Lasso, we know

kβ

−

βk

≤ c

log p/n.

for some constant

. We now seek a lower bound for

ˆτ

. Consider linear models

= X

−j

(j)

+ ε

(j)

where the sparsity of γ

(j)

is s

, and ε

(j)

−j

∼ N(0, Ω

−1

). Note that

Ω

−1

= var(X

| X

i,−j

) ≤ var(X

) = Σ

≤ A.

Also, the maximum eigenvalue of Ω is at most

−1

min

. So Ω

≤ c

−1

min

. So

Ω

−1

≥ c

min

. So by Lasso theory, we know

kγ

(j)

− ˆγ

(j)

≤ c

log p

for some constant c

. Then we have

ˆτ

− X

−j

ˆγ

(j)

+ λkˆγ

(j)

≥

kε

(j)

+ X

−j

(γ

(j)

− ˆγ

(j)

≥

kε

(j)

−

−j

(j)

∞

kγ

(j)

− ˆγ

(j)

≥ Ω

−1

1 − 4

log p

− c

log p

+ A

log p

In the limit, this tends to Ω

−1

. So for large n, this is ≥

Ω

−1

≥

min

Thus, we have

k∆k

∞

≤ 2λ

√

log p

−1

min

= A

log p

√