IAS/Park City Mathematics Series
Volume 00, Pages 000–000
S 1079-5634(XX)0000-0

Four lectures on probabilistic methods for data science

arXiv:1612.06661v2 [math.PR] 4 Nov 2017

Roman Vershynin
Abstract. Methods of high-dimensional probability play a central role in applications for statistics, signal processing, theoretical computer science and related fields. These lectures present a sample of particularly useful tools of highdimensional probability, focusing on the classical and matrix Bernstein’s inequality and the uniform matrix deviation inequality. We illustrate these tools with
applications for dimension reduction, network analysis, covariance estimation,
matrix completion and sparse signal recovery. The lectures are geared towards
beginning graduate students who have taken a rigorous course in probability but
may not have any experience in data science applications.

Contents
1

2

3

4

Lecture 1: Concentration of sums of independent random variables
1.1 Sub-gaussian distributions
1.2 Hoeffding’s inequality
1.3 Sub-exponential distributions
1.4 Bernstein’s inequality
1.5 Sub-gaussian random vectors
1.6 Johnson-Lindenstrauss Lemma
1.7 Notes
Lecture 2: Concentration of sums of independent random matrices
2.1 Matrix calculus
2.2 Matrix Bernstein’s inequality
2.3 Community recovery in networks
2.4 Notes
Lecture 3: Covariance estimation and matrix completion
3.1 Covariance estimation
3.2 Norms of random matrices
3.3 Matrix completion
3.4 Notes
Lecture 4: Matrix deviation inequality
4.1 Gaussian width

2
3
4
5
5
7
7
9
9
9
11
14
18
19
19
23
25
28
29
30

Received by the editors November 7, 2017.
Partially supported by NSF Grant DMS 1265782 and U.S. Air Force Grant FA9550-14-1-0009.
©0000 (copyright holder)

1

2

Four lectures on probabilistic methods for data science

4.2
4.3
4.4
4.5
4.6
4.7

Matrix deviation inequality
Deriving Johnson-Lindenstrauss Lemma
Covariance estimation
Underdetermined linear equations
Sparse recovery
Notes

31
32
33
35
37
38

1. Lecture 1: Concentration of sums of independent random variables
These lectures present a sample of modern methods of high dimensional probability and illustrate these methods with applications in data science. This sample
is not comprehensive by any means, but it could serve as a point of entry into
a branch of modern probability that is motivated by a variety of data-related
problems.
To get the most out of these lectures, you should have taken a graduate course
in probability, have a good command of linear algebra (including the singular
value decomposition) and be familiar with very basic concepts of functional analysis (familiarity with Lp norms should be enough).
All of the material of these lectures is covered more systematically, at a slower
pace, and with a wider range of applications, in my forthcoming textbook [60].
You may also be interested in two similar tutorials: [58] is focused on random
matrices, and a more advanced text [59] discusses high-dimensional inference
problems.
It should be possible to use these lectures for a self-study or group study. You
will find here many places where you are invited to do some work (marked in
the text e.g. by “check this!”), and you are encouraged to do it to get a better
grasp of the material. Each lecture ends with a section called “Notes” where you
will find bibliographic references of the results just discussed, as well asvarious
improvements and extensions.
We are now ready to start.
Probabilistic reasoning has a major impact on modern data science. There are
roughly two ways in which this happens.
• Radnomized algorithms, which perform some operations at random, have
long been developed in computer science and remain very popular. Randomized algorithms are among the most effective methods – and sometimes the only known ones – for many data problems.
• Random models of data form the usual premise of statistical analysis. Even
when the data at hand is deterministic, it is often helpful to think of it as a
random sample drawn from some unknown distribution (“population”).
In these lectures, we will encounter both randomized algorithms and random
models of data.

Roman Vershynin

3

1.1. Sub-gaussian distributions Before we start discussing probabilistic methods, we will introduce an important class of probability distributions that forms a
natural “habitat” for random variables in many theoretical and applied problems.
These are sub-gaussian distributions. As the name suggests, we will be looking
at an extension of the most fundamental distribution in probability theory – the
gaussian, or normal, distribution N(µ, σ).
It is a good exercise to check that the standard normal random variable X ∼
N(0, 1) satisfies the following basic properties:

Tails: P |X| > t 6 2 exp(−t2 /2) for all t > 0.
√
Moments: kXkp := (E |X|p )1/p = O( p) as p → ∞.
MGF of square: 1 E exp(cX2 ) 6 2 for some c > 0.
MGF: E exp(λX) = exp(λ2 ) for all λ ∈ R.
All these properties tell the same story from four different perspectives. It is
not very difficult to show (although we will not do it here) that for any random
variable X, not necessarily Gaussian, these four properties are essentially equivalent.
Proposition 1.1.1 (Sub-gaussian properties). For a random variable X, the following
properties are equivalent.2

Tails: P |X| > t 6 2 exp(−t2 /K21 ) for all t > 0.
√
Moments: kXkp 6 K2 p for all p > 1.
MGF of square: E exp(X2 /K23 ) 6 2.
Moreover, if E X = 0 then these properties are also equivalent to the following one:
MGF: E exp(λX) 6 exp(λ2 K24 ) for all λ ∈ R.
Random variables that satisfy one of the first three properties (and thus all of
them) are called sub-gaussian. The best K3 is called the sub-gaussian norm of X, and
is usually denoted kXkψ2 , that is
kXkψ2 := inf t > 0 : E exp(X2 /t2 ) 6 2 .
One can check that k · kψ2 indeed defines a norm; it is an example of the general
concept of the Orlicz norm. Proposition 1.1.1 states that the numbers Ki in all four
properties are equivalent to kXkψ2 up to absolute constant factors.
Example 1.1.2. As we already noted, the standard normal random variable X ∼
N(0, 1) is sub-gaussian. Similarly, arbitrary normal random variables X ∼ N(µ, σ)
are sub-gaussian. Another example is a Bernoulli random variable X that takes
values 0 and 1 with probabilities 1/2 each. More generally, any bounded random
variable X is sub-gaussian. On the contrary, Poisson, exponential, Pareto and
1MGF stands for moment generation function.
2The parameters K > 0 appearing in these properties can be different. However, they may differ
i

from each other by at most an absolute constant factor. This means that there exists an absolute
constant C such that property 1 implies property 2 with parameter K2 6 CK1 , and similarly for
every other pair or properties.

4

Four lectures on probabilistic methods for data science

Cauchy distributions are not sub-gaussian. (Verify all these claims; this is not
difficult.)
1.2. Hoeffding’s inequality You may remember from a basic course in probability that the normal distribution N(µ, σ) has a remarkable property: the sum of
independent normal random variables is also normal. Here is a version of this
property for sub-gaussian distributions.
Proposition 1.2.1 (Sums of sub-gaussians). Let X1 , . . . , XN be independent, mean
P
zero, sub-gaussian random variables. Then N
i=1 Xi is a sub-gaussian, and
N
X

Xi

i=1

where C is an absolute constant.

2
ψ2

6C

N
X

kXi k2ψ2

i=1

3

Proof. Let us bound the moment generating function of the sum for any λ ∈ R:
E exp λ

N
X



Xi =

i=1

N
Y

E exp(λXi )

(using independence)

i=1

6

N
Y

exp(Cλ2 kXi k2ψ2 )

(by last property in Proposition 1.1.1)

i=1

= exp(λ2 K2 )

where K2 := C

N
X

kXi k2ψ2 .

i=1

Using again the last property in Proposition 1.1.1, we conclude that the sum
P
S= N
i=1 Xi is sub-gaussian, and kSkψ2 6 C1 K where C1 is an absolute constant.
The proof is complete.

Let us rewrite Proposition 1.2.1 in a form that is often more useful in applications, namely as a concentration inequality. To do this, we simply use the first
P
property in Proposition 1.1.1 for the sum N
i=1 Xi . We immediately get the following.
Theorem 1.2.2 (General Hoeffding’s inequality). Let X1 , . . . , XN be independent,
mean zero, sub-gaussian random variables. Then, for every t > 0 we have
P

N
X



ct2
Xi > t 6 2 exp − PN
.
2
i=1 kXi kψ2
i=1

Hoeffding’s inequality controls how far and with what probability a sum of
independent random variables can deviate from its mean, which is zero.
3In the future, we will always denote positive absolute constants by C, c, C , etc. These numbers do
1

not depend on anything. In most cases, one can get good bounds on these constants from the proof,
but the optimal constants for each result are rarely known.

Roman Vershynin

5

1.3. Sub-exponential distributions Sub-gaussian distributions form a sufficiently
wide class of distributions. Many results in probability and data science are
proved nowadays for sub-gaussian random variables. Still, as we noted, there
are some natural random variables that are not sub-gaussian. For example, the
square X2 of a normal random variable X ∼ N(0, 1) is not sub-gaussian. (Check!)
To cover examples like this, we will introduce the similar but weaker notion of
sub-exponential distributions.
Proposition 1.3.1 (Sub-exponential properties). For a random variable X, the following properties are equivalent, in the same sense as in Proposition 1.1.1.

Tails: P |X| > t 6 2 exp(−t/K1 ) for all t > 0.
Moments: kXkp 6 K2 p for all p > 1.
MGF of the square: E exp(|X|/K3 ) 6 2.
Moreover, if E X = 0 then these properties imply the following one:
MGF: E exp(λX) 6 exp(λ2 K24 ) for |λ| 6 1/K4 .
Just like we did for sub-gaussian distributions, we call the best K3 the subexponential norm of X and denote it by kXkψ1 , that is
kXkψ1 := inf {t > 0 : E exp(|X|/t) 6 2} .
All sub-exponential random variables are squares of sub-gaussian random variables. Indeed, inspecting the definitions you will quickly see that
kX2 kψ1 = kXk2ψ2 .

(1.3.2)
(Check!)

1.4. Bernstein’s inequality A version of Hoeffding’s inequality for sub-exponential
random variables is called Bernstein’s inequality. You may naturally expect to see
a sub-exponential tail bound in this result. So it may come as a surprise that
Bernstein’s inequality actually has a mixture of two tails – sub-gaussian and subexponential. Let us state and prove the inequality first, and then we will comment
on the mixture of the two tails.
Theorem 1.4.1 (Bernstein’s inequality). Let X1 , . . . , XN be independent, mean zero,
sub-exponential random variables. Then, for every t > 0 we have
P

N
X

i
h

t
t2
,
.
Xi > t 6 2 exp − c min PN
2
maxi kXi kψ1
i=1 kXi kψ1
i=1

Proof. For simplicity, we will assume that K = 1 and only prove the one-sided
bound (without absolute value); the general case is not much harder. Our approach will be based on bounding the moment generating function of the sum
PN
S :=
i=1 Xi . To see how MGF can be helpful here, choose λ > 0 and use
Markov’s inequality to get


(1.4.2)
P S > t = P exp(λS) > exp(λt) 6 e−λt E exp(λS).

6

Four lectures on probabilistic methods for data science

Recall that S =
as

PN

i=1 Xi

and use independence to express the right side of (1.4.2)
e−λt

N
Y

E exp(λXi ).

i=1

(Check!) It remains to bound the MGF of each term Xi , and this is a much simpler
task. If we choose λ small enough so that
c
,
(1.4.3)
0<λ6
maxi kXi kψ1
then we can use the last property in Proposition 1.3.1 to get

E exp(λXi ) 6 exp Cλ2 kXi k2ψ1 .
Substitute into (1.4.2) and conclude that


P{S > t} 6 exp −λt + Cλ2 σ2

P
2
where σ2 = N
i=1 kXi kψ1 . The left side does not depend on λ while the right side
does. So we can choose λ that minimizes the right side subject to the constraint
(1.4.3). When this is done carefully, we obtain the tail bound stated in Bernstein’s
inequality. (Do this!)

Now, why does Bernstein’s inequality have a mixture of two tails? The subexponential tail should of course be there. Indeed, even if the entire sum consisted
of a single term Xi , the best bound we could hope for would be of the form
exp(−ct/kXi kψ1 ). The sub-gaussian term could be explained by the central limit
theorem, which states that the sum should becomes approximately normal as the
number of terms N increases to infinity.
Remark 1.4.4 (Bernstein’s inequality for bounded random variables). Suppose the
random variables Xi are uniformly bounded, which is a stronger assumption than
being sub-gaussian. Then there is a useful version of Bernstein’s inequality, which
unlike Theorem 1.4.1 is sensitive to the variances of Xi ’s. It states that if K > 0 is
such that |Xi | 6 K almost surely for all i, then, for every t > 0, we have
P

(1.4.5)
σ2

PN

2
i=1 E Xi

N
X
i=1


Xi > t 6 2 exp −

t2 /2 
.
σ2 + CKt

Here
=
is the variance of the sum. This version of Bernstein’s
inequality can be proved in essentially the same way as Theorem 1.4.1. We will
not do it here, but a stronger Theorem 2.2.1, which is valid for matrix-valued
random variables Xi , will be proved in Lecture 2.
To compare this with Theorem 1.4.1, note that σ2 + CKt 6 2 max(σ2 , CKt). So
we can state the probability bound (1.4.5) as
h
 t2 t i
2 exp − c min 2 ,
.
σ K

Roman Vershynin

7

Just like before, here we also have a mixture of two tails, sub-gaussian and
sub-exponential. The sub-gaussian tail is a bit sharper than in Theorem 1.4.1,
since it depends on the variances rather than sub-gaussian norms of Xi . The subexponential tail, on the other hand, is weaker, since it depends on the sup-norms
rather than the sub-exponential norms of Xi .
1.5. Sub-gaussian random vectors The concept of sub-gaussian distributions
can be extended to higher dimensions. Consider a random vector X taking values
in Rn . We call X a sub-gaussian random vector if all one-dimensional marginals of X,
i.e., the random variables hX, xi for x ∈ Rn , are sub-gaussian. The sub-gaussian
norm of X is defined as
kXkψ2 := sup k hX, xi kψ2
x∈Sn−1

where Sn−1 denotes the unit Euclidean sphere in Rn .
Example 1.5.1. Examples of sub-gaussian random distributions in Rn include the
standard normal distribution N(0, In ) (why?), the uniform distribution on the
√
centered Euclidean sphere of radius n, the uniform distribution on the cube
{−1, 1}n , and many others. The last example can be generalized: a random vector
X = (X1 , . . . , Xn ) with independent and sub-gaussian coordinates is sub-gaussian,
with kXkψ2 6 C maxi kXi kψ2 .
1.6. Johnson-Lindenstrauss Lemma Concentration inequalities like Hoeffding’s
and Bernstein’s are successfully used in the analysis of algorithms. Let us give
one example for the problem of dimension reduction. Suppose we have some data
that is represented as a set of N points in Rn . (Think, for example, of n gene
expressions of N patients.)
We would like to compress the data by representing it in a lower dimensional
space Rm instead of Rn with m  n. By how much can we reduce the dimension
without loosing the important features of the data?
The basic result in this direction is the Johnson-Lindenstrauss Lemma. It states
that a remarkably simple dimension reduction method works – a random linear
map from Rn to Rm with
m ∼ log N,
see Figure 1.6.3. The logarithmic function grows very slowly, so we can usually
reduce the dimension dramatically.
What exactly is a random linear map? Several models are possible to use.
Here we will model such a map using a Gaussian random matrix – an m × n
matrix A with independent N(0, 1) entries. More generally, we can consider an
m × n matrix A whose rows are independent, mean zero, isotropic4 and subgaussian random vectors in Rn . For example, the entries of A can be independent
Rademacher entries – those taking values ±1 with equal probabilities.
4A random vector X ∈ Rn is called isotropic if E XXT = I .
n

8

Four lectures on probabilistic methods for data science

Theorem 1.6.1 (Johnson-Lindenstrauss Lemma). Let X be a set of N points in Rn
and ε ∈ (0, 1). Consider an m × n matrix A whose rows are independent, mean zero,
isotropic and sub-gaussian random vectors in Rn . Rescale A by defining the “Gaussian
random projection”5
1
P := √ A.
m
Assume that
m > Cε−2 log N,
where C is an appropriately large constant that depends only on the sub-gaussian norms of
the vectors Xi . Then, with high probability (say, 0.99), the map P preserves the distances
between all points in X with error ε, that is
(1.6.2)

(1 − ε)kx − yk2 6 kPx − Pyk2 6 (1 + ε)kx − yk2

for all x, y ∈ X.

Figure 1.6.3. Johnson-Lindenstrauss Lemma states that a random projection of N data points from dimension n to dimension m ∼ log N
approximately preserves the distances between the points.

Proof. Take a closer look at the desired conclusion (1.6.2). By linearity, Px − Py =
P(x − y). So, dividing the inequality by kx − yk2 , we can rewrite (1.6.2) in the
following way:
(1.6.4)

1 − ε 6 kPzk2 6 1 + ε for all z ∈ T

where
T :=

x−y
: x, y ∈ X distinct points .
kx − yk2

It will be convenient to square the inequality (1.6.4). Using that 1 + ε 6 (1 + ε)2
and 1 − ε > (1 − ε)2 , we see that it is enough to show that
(1.6.5)

1 − ε 6 kPzk22 6 1 + ε for all z ∈ T .

By construction, the coordinates of the vector Pz =
we can restate (1.6.5) as

√1 Az
m

are

√1
m

hXi , zi. Thus

1 X
hXi , zi2 − 1 6 ε for all z ∈ T .
m
m

(1.6.6)

i=1

5Strictly speaking, this P is not a projection since it maps Rn to a different space Rm .

Roman Vershynin

9

Results like (1.6.6) are often proved by combining concentration and a union
bound. In order to use concentration, we first fix z ∈ T . By assumption, the
random variables hXi , zi2 − 1 are independent; they have zero mean (use isotropy
to check this!), and they are sub-exponential (use (1.3.2) to check this). Then
Bernstein’s inequality (Theorem 1.4.1) gives


m
1 X
2
hXi , zi − 1 > ε 6 2 exp(−cε2 m).
P
m
i=1

(Check!)
Finally, we can unfix z by taking a union bound over all possible z ∈ T :




m
m
X
1 X
1 X
2
2
hXi , zi − 1 > ε 6
hXi , zi − 1 > ε
P max
P
m
z∈T m
z∈T

i=1

i=1
2

6 |T | · 2 exp(−cε m).

(1.6.7)

By definition of T , we have |T | 6
So, if we choose m > Cε−2 log N with
appropriately large constant C, we can make (1.6.7) bounded by 0.01. The proof
is complete.

N2 .

1.7. Notes The material presented in Sections 1.1–1.5 is basic and can be found
e.g. in [58] and [60] with all the proofs. Bernstein’s and Hoeffding’s inequalities
that we covered here are two basic examples of concentration inequalities. There
are many other useful concentration inequalities for sums of independent random
variables (e.g. Chernoff’s and Bennett’s) and for more general objects. The textbook [60] is an elementary introduction into concentration; the books [10, 38, 39]
offer more comprehensive and more advanced accounts of this area.
The original version of Johnson-Lindenstrauss Lemma was proved in [31]. The
version we gave here, Theorem 1.6.1, was stated with probability of success 0.99,
but an inspection of the proof gives probability 1 − 2 exp(−cε2 m) which is much
better for large m. A great variety of ramifications and applications of JohnsonLindenstrauss lemma are known, see e.g. [2, 4, 7, 10, 34, 42].

2. Lecture 2: Concentration of sums of independent random matrices
In the previous lecture we proved Bernstein’s inequality, which quantifies how
a sum of independent random variables concentrates about its mean. We will
now study an extension of Bernstein’s inequality to higher dimensions, which
holds for sums of independent random matrices.
2.1. Matrix calculus The key idea of developing a matrix Bernstein’s inequality
will be to use matrix calculus, which allows us to operate with matrices as with
scalars – adding and multiplying them of course, but also comparing matrices
and applying functions to matrices. Let us explain this.

10

Four lectures on probabilistic methods for data science

We can compare matrices to each other using the notion of being positive semidefinite. Let us focus here on n × n symmetric matrices. If A − B is a positive
semidefinite matrix,6 which we denote A − B  0, then we say that A  B (and,
of course, B  A). This defines a partial order on the set of n × n symmetric matrices. The term “partial” indicates that, unlike the real numbers, there exist n × n
symmetric matrices A and B that can not be compared. (Give an example where
neither A  B nor B  A!)
Next, let us guess how to measure the magnitude of a matrix A. The magnitude
of a scalar a ∈ R is measured by the absolute value |a|; it is the smallest nonnegative number t such that
−t 6 a 6 t.
Extending this reasoning to matrices, we can measure the magnitude of an n × n
symmetric matrix A by the smallest non-negative number t such that7
−tIn  A  tIn .
The smallest t is called the operator norm of A and is denoted kAk. Diagonalizing
A, we can see that
(2.1.1)

kAk = max{|λ| : λ is an eigenvalue of A}.

With a little more work (do it!), we can see that kAk is the norm of A acting as a
linear operator on the Euclidean space (Rn , k · k2 ); this is why kAk is called the
operator norm. Thus kAk is the smallest non-negative number M such that
kAxk2 6 Mkxk2

for all x ∈ Rn .

Finally, we will need to be able to take functions of matrices. Let f : R → R
be a function and X be an n × n symmetric matrix. We can define f(X) in two
equivalent ways. The spectral theorem allows us to represent X as
X=

n
X

λi u i u T
i

i=1

where λi are the eigenvalues of X and ui are the corresponding eigenvectors.
Then we can simply define
f(X) :=

n
X

f(λi )ui uT
i.

i=1

Note that f(X) has the same eigenvectors as X, but the eigenvalues change under
the action of f. An equivalent way to define f(X) is using power series. Suppose
the function f has a convergent power series expansion about some point x0 ∈ R,
6Recall that a symmetric real n × n matrix M is called positive semidefinite if xT Mx > 0 for any
vector x ∈ Rn .
7Here and later, I denotes the n × n identity matrix.
n

Roman Vershynin

i.e.

∞
X

f(x) =

11

ak (x − x0 )k .

k=1

Then one can check that the following matrix series converges8 and defines f(X):
f(X) =

∞
X

ak (X − X0 )k .

k=1

(Check!)
2.2. Matrix Bernstein’s inequality We are now ready to state and prove a remarkable generalization of Bernstein’s inequality for random matrices.
Theorem 2.2.1 (Matrix Bernstein’s inequality). Let X1 , . . . , XN be independent, mean
zero, n × n symmetric random matrices, such that kXi k 6 K almost surely for all i.
Then, for every t > 0 we have
P
Here

σ2

=

PN

N
X


Xi > t 6 2n · exp −

i=1

2
i=1 E Xi

t2 /2 
.
σ2 + Kt/3

is the norm of the “matrix variance” of the sum.

The scalar case, where n = 1, is the classical Bernstein’s inequality we stated
in (1.4.5). A remarkable feature of matrix Bernstein’s inequality, which makes
it especially powerful, is that it does not require any independence of the entries (or
the rows or columns) of Xi ; all is needed is that the random matrices Xi be
independent from each other.
In the rest of this section we will prove matrix Bernstein’s inequality, and give
a few applications in this and next lecture.
Our proof will be based on bounding the moment generating function (MGF)
P
E exp(λS) of the sum S = N
i=1 Xi . Note that to exponentiate the matrix λS in
order to define the matrix MGF, we rely on the matrix calculus that we introduced
in Section 2.1.
If the terms Xi were scalars, independence would yield the classical fact that
the MGF of a product is the product of MGF’s, i.e.
(2.2.2)

E exp(λS) = E

N
Y

exp(λXi ) =

i=1

N
Y

E exp(λXi ).

i=1

But for matrices, this reasoning breaks down badly, for in general
eX+Y 6= eX eY
even for 2 × 2 symmetric matrices X and Y. (Give a counterexample!)
Fortunately, there are some trace inequalities that can often serve as proxies for
the missing equality eX+Y = eX eY . One of such proxies is the Golden-Thompson
8The convergence holds in any given metric on the set of matrices, for example in the metric given by
the operator norm. In this series, the terms (X − X0 )k are defined by the usual matrix product.

12

Four lectures on probabilistic methods for data science

inequality, which states that
(2.2.3)

tr(eX+Y ) 6 tr(eX eY )

for any n × n symmetric matrices X and Y. Another result, which we will actually
use in the proof of matrix Bernstein’s inequality, is Lieb’s inequality.
Theorem 2.2.4 (Lieb’s inequality). Let H be an n × n symmetric matrix. Then the
function
f(X) = tr exp(H + log X)
is concave9 on the space on n × n symmetric matrices.
Note that in the scalar case, where n = 1, the function f in Lieb’s inequality is
linear and the result is trivial.
To use Lieb’s inequality in a probabilistic context, we will combine it with
the classical Jensen’s inequality. It states that for any concave function f and a
random matrix X, one has10
(2.2.5)

E f(X) 6 f(E X).

Using this for the function f in Lieb’s inequality, we get
E tr exp(H + log X) 6 tr exp(H + log E X).
And changing variables to X = eZ , we get the following:
Lemma 2.2.6 (Lieb’s inequality for random matrices). Let H be a fixed n × n symmetric matrix and Z be an n × n symmetric random matrix. Then
E tr exp(H + Z) 6 tr exp(H + log E eZ ).
Lieb’s inequality is a perfect tool for bounding the MGF of a sum of indepenP
dent random variables S = N
random
i=1 Xi . To do this, let us condition on the
P
variables X1 , . . . , XN−1 . Apply Lemma 2.2.6 for the fixed matrix H := N−1
i=1 λXi
and the random matrix Z := λXN , and afterwards take the expectation with respect to X1 , . . . , XN−1 . By the law of total expectation, we get
E tr exp(λS) 6 E tr exp

h N−1
X

i
λXi + log E eλXN .

i=1

Next, apply Lemma 2.2.6 in a similar manner for H :=
and Z := λXN−1 , and so on. After N times, we obtain:

PN−2
i=1

λXi + log E eλXN

9Formally, concavity of f means that f(λX + (1 − λ)Y) > λf(X) + (1 − λ)f(Y) for all symmetric
matrices X and Y and all λ ∈ [0, 1].
10Jensen’s inequality is usually stated for a convex function g and a scalar random variable X, and
it reads g(E X) 6 E g(X). From this, inequality (2.2.5) for concave functions and random matrices
easily follows (Check!).

Roman Vershynin

13

Lemma 2.2.7 (MGF of a sum of independent random matrices). Let X1 , . . . , XN be
P
independent n × n symmetric random matrices. Then the sum S = N
i=1 Xi satisfies
E tr exp(λS) 6 tr exp

N
hX

i
log E eλXi .

i=1

Think of this inequality is a matrix version of the scalar identity (2.2.2). The
main difference is that it bounds the trace of the MGF11 rather the MGF itself.
You may recall from a course in probability theory that the quantity log E eλXi
that appears in this bound is called the cumulant generating function of Xi .
Lemma 2.2.7 reduces the complexity of our task significantly, for it is much
easier to bound the cumulant generating function of each single random variable
Xi than to say something about their sum. Here is a simple bound.
Lemma 2.2.8 (Moment generating function). Let X be an n × n symmetric random
matrix. Assume that E X = 0 and kXk 6 K almost surely. Then, for all 0 < λ < 3/K we
have


λ2 /2
.
E exp(λX)  exp g(λ) E X2
where g(λ) =
1 − λK/3
Proof. First, check that the following scalar inequality holds for 0 < λ < 3/K and
|x| 6 K:
eλx 6 1 + λx + g(λ)x2 .
Then extend it to matrices using matrix calculus: if 0 < λ < 3/K and kXk 6 K
then
eλX  I + λX + g(λ)X2 .
(Do these two steps carefully!) Finally, take the expectation and recall that E X = 0
to obtain


E eλX  I + g(λ) E X2  exp g(λ) E X2 .
In the last inequality, we use the matrix version of the scalar inequality 1 + z 6 ez
that holds for all z ∈ R. The lemma is proved.

Proof of Matrix Bernstein’s inequality. We would like to bound the operator norm
P
of the random matrix S = N
i=1 Xi , which, as we know from (2.1.1), is the largest
eigenvalue of S by magnitude. For simplicity of exposition, let us drop the absolute
value from (2.1.1) and just bound the maximal eigenvalue of S, which we denote
λmax (S). (Once this is done, we can repeat the argument for −S to reinstate the
absolute value. Do this!) So, we are to bound

P λmax (S) > t = P eλ·λmax (S) > eλt
(multiply by λ > 0 and exponentiate)
6 e−λt E eλ·λmax (S)

(by Markov’s inequality)

= e−λt E λmax (eλS )

(check!)

6e

−λt

E tr e

λS

(max of eigenvalues is bounded by the sum)

11Note that the order of expectation and trace can be swapped using linearity.

14

Four lectures on probabilistic methods for data science

6 e−λt tr exp

N
hX

log E eλXi

i

(use Lemma 2.2.7)

i=1

6 tr exp [−λt + g(λ)Z]
where
Z :=

(by Lemma 2.2.8)

N
X

E X2i .

i=1

It remains to optimize this bound in λ. The minimum is attained for λ =
t/(σ2 + Kt/3). (Check!) With this value of λ, we conclude


t2 /2 
.
P λmax (S) > t 6 n · exp − 2
σ + Kt/3
This completes the proof of Theorem 2.2.1.

PN
Bernstein’s inequality gives a powerful tail bound for k i=1 Xi k. This easily
implies a useful bound on the expectation:
Corollary 2.2.9 (Expected norm of sum of random matrices). Let X1 , . . . , XN be
independent, mean zero, n × n symmetric random matrices, such that kXi k 6 K almost
surely for all i. Then
N
X
p
Xi . σ log n + K log n
E
where σ =

PN

i=1

2 1/2 .
i=1 E Xi

Proof. The link from tail bounds to expectation is provided by the basic identity
Z∞

(2.2.10)
EZ =
P Z > t dt
0

which is valid for any non-negative random variable Z. (Check it!) Integrating
the tail bound given by matrix Bernstein’s inequality, you will arrive at the expectation bound we claimed. (Check!)

Notice that the bound in this corollary has mild, logarithmic, dependence on
the ambient dimension n. As we will see shortly, this can be an important feature
in some applications.
2.3. Community recovery in networks Matrix Bernstein’s inequality has many
applications. The one we are going to discuss first is for the analysis of networks.
A network can be mathematically represented by a graph, a set of n vertices
with edges connecting some of them. For simplicity, we will consider undirected
graphs where the edges do not have arrows. Real world networks often tend to
have clusters, or communities – subsets of vertices that are connected by unusually
many edges. (Think, for example, about a friendship network where communities
form around some common interests.) An important problem in data science is
to recover communities from a given network.

Roman Vershynin

15

We are going to explain one of the simplest methods for community recovery,
which is called spectral clustering. But before we introduce it, we will first of all
place a probabilistic model on the networks we consider. In other words, it will be
convenient for us to view networks as random graphs whose edges are formed at
random. Although not all real-world networks are truly random, this simplistic
model can motivate us to develop algorithms that may empirically succeed also
for real-world networks.
The basic probabilistic model of random graphs is the Erdös-Rényi model.
Definition 2.3.1 (Erdös-Rényi model). Consider a set of n vertices and connect every
pair of vertices independently and with fixed probability p. The resulting random graph
is said to follow the Erdös-Rényi model G(n, p).
The Erdös-Rényi random model is very simple. But it is not a good choice if
we want to model a network with communities, for every pair of vertices has the
same chance to be connected. So let us introduce a natural generalization of the
Erdös-Rényi random model that does allow for community structure:
Definition 2.3.2 (Stochastic block model). Partition a set of n vertices into two subsets
(“communities”) with n/2 vertices each, and connect every pair of vertices independently
with probability p if they belong to the same community and q < p if not. The resulting
random graph is said to follow the stochastic block model G(n, p, q).
Figure 2.3.3 illustrates a simulation of a stochastic block model.

Figure 2.3.3. A network generated according to the stochastic block
model G(n, p, q) with n = 200 nodes and connection probabilities
p = 1/20 and q = 1/200.

Suppose we are shown one instance of a random graph generated according
to a stochastic block model G(n, p, q). How can we find which vertices belong to
which community?
The spectral clustering algorithm we are going to explain will do precisely this.
It will be based on the spectrum of the adjacency matrix A of the graph, which is

16

Four lectures on probabilistic methods for data science

the n × n symmetric matrix whose entries Aij equal 1 if the vertices i and j are
connected by an edge, and 0 otherwise.12
The adjacency matrix A is a random matrix. Let us compute its expectation
first. This is easy, since the entires of A are Bernoulli random variables. If i and j
belong to the same community then E Aij = p and otherwise E Aij = q. Thus A
has block structure: for example, if n = 4 then A looks like this:


p p q q


 p p q q 


EA = 

 q q p p 


q q p p
(For illustration purposes, we grouped the vertices from each community together.)
You will easily check that A has rank 2, and the non-zero eigenvalues and the
corresponding eigenvectors are
(2.3.4)



1
1



 1 
 1
p + q
p − q



n, v1 (E A) = 
n, v2 (E A) = 
λ1 (E A) =
 ; λ2 (E A) =
 1 
 −1
2
2



1
−1
(Check!)
The eigenvalues and eigenvectors of E A tell us a lot about the community
structure of the underlying graph. Indeed, the first (larger) eigenvalue,
p + q
d :=
n,
2
is the expected degree of any vertex of the graph.13 The second eigenvalue tells
us whether there is any community structure at all (which happens when p 6= q
and thus λ2 (E A) 6= 0). The first eigenvector v1 is not informative of the structure
of the network at all. It is the second eigenvector v2 that tells us exactly how to
separate the vertices into the two communities: the signs of the coefficients of v2
can be used for this purpose.
Thus if we know E A, we can recover the community structure of the network
from the signs of the second eigenvector. The problem is that we do not know
E A. Instead, we know the adjacency matrix A. If, by some chance, A is not far
from E A, we may hope to use the A to approximately recover the community
structure. So is it true that A ≈ E A? The answer is yes, and we can prove it
using matrix Bernstein’s inequality.
12For convenience, we call the vertices of the graph 1, 2, . . . , n.
13The degree of the vertex is the number of edges connected to it.





.



Roman Vershynin

17

Theorem 2.3.5 (Concentration of the stochastic block model). Let A be the adjacency
matrix of a G(n, p, q) random graph. Then
p
E kA − E Ak . d log n + log n.
Here d = (p + q)n/2 is the expected degree.
Proof. Let us sketch the argument. To use matrix Bernstein’s inequality, let us
break A into a sum of independent random matrices
X
A=
Xij ,
i,j: i6j

where each matrix Xij contains a pair of symmetric entries of A, or one diagonal
entry.14 Matrix Bernstein’s inequality obviously applies for the sum
X
A−EA =
(Xij − E Xij ).
i6j
15

Corollary 2.2.9 gives

p
E kA − E Ak . σ log n + K log n

(2.3.6)

P
2 and K = max kX − E X k. It is a good
where σ2 =
ij
ij
ij
i6j E(Xij − E Xij )
exercise to check that
σ2 . d and K 6 2.
(Do it!) Substituting into (2.3.6), we complete the proof.



How useful is Theorem 2.3.5 for community recovery? Suppose that the network is not too sparse, namely
d  log n.
Then
kA − E Ak .

p

d log n while

k E Ak = λ1 (E A) = d,

which implies that
kA − E Ak  k E Ak.
In other words, A nicely approximates E A: the relative error or approximation
is small in the operator norm.
At this point one can apply classical results from the perturbation theory for
matrices, which state that since A and E A are close, their eigenvalues and eigenvectors must also be close. The relevant perturbation results are Weyl’s inequality
for eigenvalues and Davis-Kahan’s inequality for eigenvectors, which we will not
14Precisely, if i 6= j, then X

ij has all zero entries except the (i, j) and (j, i) entries that can potentially
equal 1. If i = j, the only non-zero entry of Xij is the (i, i).
15We will liberally use the notation . to hide constant factors appearing in the inequalities. Thus,
a . b means that a 6 Cb for some constant C.

18

Four lectures on probabilistic methods for data science

reproduce here. Heuristically, what they give us is

1

 1

v2 (A) ≈ v2 (E A) = 
 −1

−1





.



Then we should expect that most of the coefficients of v2 (A) are positive on one
community and negative on the other. So we can use v2 (A) to approximately
recover the communities. This method is called spectral clustering:
Spectral Clustering Algorithm. Compute v2 (A), the eigenvector corresponding to
the second largest eigenvalue of the adjacency matrix A of the network. Use the signs of
the coefficients of v2 (A) to predict the community membership of the vertices.
We saw that spectral clustering should perform well for the stochastic block
model G(n, p, q) if it is not too sparse, namely if the expected degrees satisfy
d = (p + q)n/2  log n.
A more careful analysis along these lines, which you should be able to do
yourself with some work, leads to the following more rigorous result.
Theorem 2.3.7 (Guarantees of spectral clustering). Consider a random graph generated according to the stochastic block model G(n, p, q) with p > q, and set a = pn,
b = qn. Suppose that
(2.3.8)

(a − b)2  log(n)(a + b).

Then, with high probability, the spectral clustering algorithm recovers the communities
up to o(n) misclassified vertices.
Note that condition (2.3.8) implies that the expected degrees are not too small,
namely d = (a + b)/2  log(n) (check!). It also ensures that a and b are sufficiently different: recall that if a = b the network is Erdös-Rényi graph without
any community structure.
2.4. Notes The idea to extend concentration inequalities like Bernstein’s to matrices goes back to R. Ahlswede and A. Winter [3]. They used Golden-Thompson
inequality (2.2.3) and proved a slightly weaker form of matrix Bernstein’s inequality than we gave in Section 2.2. R. Oliveira [48, 49] found a way to improve
this argument and gave a result similar to Theorem 2.2.1. The version of matrix
Bernstein’s inequality we gave here (Theorem 2.2.1) and a proof based on Lieb’s
inequality is due to J. Tropp [52].
The survey [53] contains a comprehensive introduction of matrix calculus, a
proof of Lieb’s inequality (Theorem 2.2.4), a detailed proof of matrix Bernstein’s
inequality (Theorem 2.2.1) and a variety of applications. A proof of GoldenThompson inequality (2.2.3) can be found in [8, Theorem 9.3.7].

Roman Vershynin

19

In Section 2.3 we scratched the surface of an interdisciplinary area of network analysis. For a systematic introduction into networks, refer to the book
[47]. Stochastic block models (Definition 2.3.2) were introduced in [33]. The
community recovery problem in stochastic block models, sometimes also called
community detection problem, has been in the spotlight in the last few years. A
vast and still growing body of literature exists on algorithms and theoretical results for community recovery, see the book [47], the survey [22], papers such as
[9, 29, 30, 32, 37, 46, 61] and the references therein.
A concentration result similar to Theorem 2.3.5 can be found in [48]; the argument there is also based on matrix concentration. This theorem is not quite
optimal. For dense networks, where the expected degree d satisfies d & log n, the
concentration inequality in Theorem 2.3.5 can be improved to
√
(2.4.1)
E kA − E Ak . d.
This improved bound goes back to the original paper [21] which studies the simpler Erdös-Rényi model but the results extend to stochastic block models [17]; it
can also be deduced from [6, 32, 37].
If the network is relatively dense, i.e. d & log n, one can improve the guarantee
(2.3.8) of spectral clustering in Theorem 2.3.7 to
(a − b)2  (a + b).
All one has to do is use the improved concentration inequality (2.4.1) instead of
Theorem 2.3.5. Furthermore, in this case there exist algorithms that can recover
the communities exactly, i.e. without any misclassified vertices, and with high
probability, see e.g. [1, 17, 32, 43].
For sparser networks, where d  log n and possibly even d = O(1), relativelyfew algorithms were known until recently, but now there exist many approaches that provably recover communities in sparse stochastic block models,
see e.g. [9, 17, 29, 30, 37, 46, 61].

3. Lecture 3: Covariance estimation and matrix completion
In the last lecture, we proved matrix Bernstein’s inequality and gave an application for network analysis. We will spend this lecture discussing a couple of
other interesting applications of matrix Bernstein’s inequality. In Section 3.1 we
will work on covariance estimation, a basic problem in high-dimensional statistics.
In Section 3.2, we will derive a useful bound on norms of random matrices, which
unlike Bernstein’s inequality does not require any boundedness assumptions on
the distribution. We will apply this bound in Section 3.3 for a problem of matrix
completion, where we are shown a small sample of the entries of a matrix and
asked to guess the missing entries.
3.1. Covariance estimation Covariance estimation is a problem of fundamental
importance in high-dimensional statistics. Suppose we have a sample of data

20

Four lectures on probabilistic methods for data science

points X1 , . . . , XN in Rn . It is often reasonable to assume that these points are
independently sampled from the same probability distribution (or “population”)
which is unknown. We would like to learn something useful about this distribution.
Denote by X a random vector that has this (unknown) distribution. The most
basic parameter of the distribution is the mean E X. One can estimate E X from
1 PN
the sample by computing the sample mean N
i=1 Xi . The law of large numbers
guarantees that the estimate becomes tight as the sample size N grows to infinity,
i.e.
N
1 X
Xi → E X as N → ∞.
N
i=1

The next most basic parameter of the distribution is the covariance matrix
Σ := E(X − E X)(X − E X)T .
This is a higher-dimensional version of the usual notion of variance of a random
variable Z, which is
Var(Z) = E(Z − E Z)2 .
The eigenvectors of the covariance matrix of Σ are called the principal components.
Principal components that correspond to large eigenvalues of Σ are the directions
in which the distribution of X is most extended, see Figure 3.1.1. These are often
the most interesting directions in the data. Practitioners often visualize the highdimensional data by projecting it onto the span of a few (maybe two or three) of
such principal components; the projection may reveal some hidden structure of
the data. This method is called Principal Component Analysis (PCA).

Figure 3.1.1. Data points X1 , . . . , XN sampled from a distribution in Rn
and the principal components of the covariance matrix.

One can estimate the covariance matrix Σ from the sample by computing the
sample covariance
N
1 X
ΣN :=
(Xi − E Xi )(Xi − E Xi )T .
N
i=1

Again, the law of large numbers guarantees that the estimate becomes tight as
the sample size N grows to infinity, i.e.
ΣN → Σ

as N → ∞.

Roman Vershynin

21

But how large should the sample size N be for covariance estimation? Generally, one can not have N < n for dimension reasons. (Why?) We are going to
show that
N ∼ n log n
is enough. In other words, covariance estimation is possible with just logarithmic
oversampling.
For simplicity, we shall state the covariance estimation bound for mean zero
distributions. (If the mean is not zero, we can estimate it from the sample and
subtract. Check that the mean can be accurately estimated from a sample of size
N = O(n).)
Theorem 3.1.2 (Covariance estimation). Let X be a random vector in Rn with covariance matrix Σ. Suppose that
kXk22 . E kXk22 = tr Σ

(3.1.3)

almost surely.

Then, for every N > 1, we have
r n log n

n log n 
.
N
N
Before we pass to the proof, let us note that Theorem 3.1.2 yields the covariance
estimation result we promised. Let ε ∈ (0, 1). If we take a sample of size
E kΣN − Σk . kΣk

+

N ∼ ε−2 n log n,
then we are guaranteed covariance estimation with a good relative error:
E kΣN − Σk 6 εkΣk.
Proof. Apply matrix Bernstein’s inequality (Corollary 2.2.9) for the sum of independent random matrices Xi XT
i − Σ and get
(3.1.4)

E kΣN − Σk =

N
X

1 p
1
(Xi XT
E
σ log n + K log n
i − Σ) .
N
N
i=1

where
σ2 =

N
X

2
E(Xi XT
= N E(XXT − Σ)2
i − Σ)

i=1

and K is chosen so that
kXXT − Σk 6 K

almost surely.

It remains to bound σ and K. Let us start with σ. We have
E(XXT − Σ)2 = E kXk22 XXT − Σ2
- tr(Σ) · E XXT

(check by expanding the square)

(drop Σ2 and use (3.1.3))

= tr(Σ) · Σ.
Thus
σ2 . N tr(Σ)kΣk.

22

Four lectures on probabilistic methods for data science

Next, to bound K, we have
kXXT − Σk 6 kXk22 + kΣk
. tr Σ + kΣk

(by triangle inequality)
(using (3.1.3))

6 2 tr Σ =: K.
Substitute the bounds on σ and K into (3.1.4) and get
q

1
E kΣN − Σk .
N tr(Σ)kΣk log n + tr(Σ) log n
N
To complete the proof, use that tr Σ 6 nkΣk (check this!) and simplify the bound.

Remark 3.1.5 (Low-dimensional distributions). Much fewer samples are needed
for covariance estimation for low-dimensional, or approximately low-dimensional,
distributions. To measure approximate low-dimensionality we can use the notion
of the stable rank of Σ2 . The stable rank of a matrix A is defined as the square of
the ratio of the Frobenius to operator norms:16
kAk2F
.
kAk2
The stable rank is always bounded by the usual, linear algebraic rank,
r(A) :=

r(A) 6 rank(A),
and it can be much smaller. (Check both claims.)
Our proof of Theorem 3.1.2 actually gives
r r log n r log n 
E kΣN − Σk 6 kΣk
+
.
N
N
where
tr Σ
r = r(Σ1/2 ) =
.
kΣk
(Check this!) Therefore, covariance estimation is possible with
N ∼ r log n
samples.
Remark 3.1.6 (The boundedness condition). It is a good exercise to check that if
we remove the boundedness condition (3.1.3), a nontrivial covariance estimation
is impossible in general. (Show this!) But how do we know whether the boundedness condition holds for data at hand? We may not, but we can enforce this
condition by truncation. All we have to do is to discard 1% of data points with
largest norms. (Check this accurately, assuming that such truncation does not
change the covariance significantly.)
16The Frobenius norm of an n × m matrix, sometimes also called the Hilbert-Schmidt norm, is de-

P
Pm
2 1/2 . Equivalently, for an n × n symmetric matrix, kAk =
fined
= ( n
F
i=1
j=1 Aij )
Pn as kAk2F 1/2
( i=1 λi (A) ) , P
where λi (A) are the eigenvalues of A. Thus the stable rank of A can be ex2
2
pressed as r(A) = n
i=1 λi (A) / maxi λi (A) .

Roman Vershynin

23

3.2. Norms of random matrices We have worked a lot with the operator norm of
matrices, denoted kAk. One may ask if is there exists a formula that expresses kAk
in terms of the entires Aij . Unfortunately, there is no such formula. The operator
norm is a more difficult quantity in this respect than the Frobenius norm, which as
P
we know can be easily expressed in terms of entries: kAkF = ( i,j A2ij )1/2 .
If we can not express kAk in terms of the entires, can we at least get a good estimate? Let us consider n × n symmetric matrices for simplicity. In one direction,
kAk is always bounded below by the largest Euclidean norm of the rows Ai :
X
1/2
(3.2.1)
kAk > max kAi k2 = max
A2ij
.
i

i

j

(Check!) Unfortunately, this bound is sometimes very loose, and the best possible
upper bound is
√
(3.2.2)
kAk 6 n · max kAi k2 .
i

(Show this bound, and give an example where it is sharp.)
Fortunately, for random matrices with independent entries the bound (3.2.2)
can be improved to the point where the upper and lower bounds almost match.
Theorem 3.2.3 (Norms of random matrices without boundedness assumptions).
Let A be an n × n symmetric random matrix whose entries on and above the diagonal are
independent, mean zero random variables. Then
E max kAi k2 6 E kAk 6 C log n · E max kAi k2 ,
i

i

where Ai denote the rows of A.
In words, the operator norm of a random matrix is almost determined by the
norm of the rows.
Our proof of this result will be based on matrix Bernstein’s inequality – more
precisely, Corollary 2.2.9. There is one surprising point. How can we use matrix
Bernstein’s inequality, which applies only for bounded distributions, to prove
a result like Theorem 3.2.3 that does not have any boundedness assumptions?
We will do this using a trick based on conditioning and symmetrization. Let us
introduce this technique first.
Lemma 3.2.4 (Symmetrization). Let X1 , . . . , XN be independent, mean zero random
vectors in a normed space and ε1 , . . . , εN be independent Rademacher random variables.17
Then
N
N
N
X
X
X
1
E
εi Xi 6 E
Xi 6 2 E
εi Xi .
2
i=1

i=1

i=1

(Xi0 )

be an independent copy of the random
Proof. To prove the upper bound, let
vectors (Xi ), i.e. just different random vectors with the same joint distribution as
17This means that random variables ε take values ±1 with probability 1/2 each. We require that all
i

random variables we consider here, i.e. {Xi , εi : i = 1, . . . , N} are jointly independent.

24

Four lectures on probabilistic methods for data science

(Xi ) and independent from (Xi ). Then
X 
X
X
Xi − E
Xi0
Xi = E
E
i

i

6E

X

Xi −

i

=E

X

(since E

X

Xi0 = 0 by assumption)

i

i

Xi0

(by Jensen’s inequality)

i

X
(Xi − Xi0 ) .
i

The distribution of the random vectors Yi := Xi − Xi0 is symmetric, which means
that the distributions of Yi and −Yi0 are the same. (Why?) Thus the distribution
of the random vectors Yi and εi Yi is also the same, for all we do is change the
signs of these vectors at random and independently of the values of the vectors.
Summarizing, we can replace Xi − Xi0 in the sum above with εi (Xi − Xi0 ). Thus
X
X
εi (Xi − Xi0 )
E
Xi 6 E
i

i

6E

X

εi Xi + E

X

X

(using triangle inequality)

i

i

= 2E

εi Xi0

εi Xi

(the two sums have the same distribution).

i

This proves the upper bound in the symmetrization inequality. The lower bound
can be proved by a similar argument. (Do this!)

Proof of Theorem 3.2.3. We already mentioned in (3.2.1) that the bound in Theorem 3.2.3 is trivial. The proof of the upper bound will be based on matrix Bernstein’s inequality.
First, we decompose A in the same way as we did in the proof of Theorem 2.3.5.
Thus we represent A as a sum of independent, mean zero, symmetric random
matrices Zij each of which contains a pair of symmetric entries of A (or one
diagonal entry):
X
A=
Zij .
i,j: i6j

Apply the symmetrization inequality (Lemma 3.2.4) for the random matrices Zij
and get
X
X
(3.2.5)
E kAk = E
Zij 6 2 E
Xij
i6j

i6j

where we set
Xij := εij Zij
and εij are independent Rademacher random variables.
Now we condition on A. The random variables Zij become fixed values and all
randomness remains in the Rademacher random variables εij . Note that Xij are
(conditionally) bounded almost surely, and this is exactly what we have lacked to

Roman Vershynin

25

apply matrix Bernstein’s inequality. Now we can do it. Corollary 2.2.9 gives18
X
p
(3.2.6)
Eε
Xij . σ log n + K log n,
σ2

P

i6j
2
i6j Eε Xij

and K = maxi6j kXij k.
where
=
A good exercise is to check that
σ . max kAi k2
i

and K . max kAi k2 .
i

(Do it!) Substituting into (3.2.6), we get
X
Xij . log n · max kAi k2 .
Eε
i

i6j

Finally, we unfix A by taking expectation of both sides of this inequality with
respect to A and using the law of total expectation. The proof is complete.

We stated Theorem 3.2.3 for symmetric matrices, but it is simple to extend it
to general m × n random matrices A. The bound in this case becomes

(3.2.7)
E kAk 6 C log(m + n) · E max kAi k2 + E max kAj k2
i

j

where Ai and Aj denote the rows and columns of A. To see this, apply Theorem 3.2.3 to the (m + n) × (m + n) symmetric random matrix
"
#
0 A
.
AT 0
(Do this!)
3.3. Matrix completion Consider a fixed, unknown n × n matrix X. Suppose we
are shown m randomly chosen entries of X. Can we guess all the missing entries?
This important problem is called matrix completion. We will analyze it using the
bounds on the norms on random matrices we just obtained.
Obviously, there is no way to guess the missing entries unless we know something extra about the matrix X. So let us assume that X has low rank:
rank(X) =: r  n.
The number of degrees of freedom of an n × n matrix with rank r is O(rn).
(Why?) So we may hope that
(3.3.1)

m ∼ rn

observed entries of X will be enough to determine X completely. But how?
Here we will analyze what is probably the simplest method for matrix completion. Take the matrix Y that consists of the observed entries of X while all
unobserved entries are set to zero. Unlike X, the matrix Y may not have small
18We stick a subscript ε to the expected value to remember that this is a conditional expectation, i.e.

we average only with respect to εi .

26

Four lectures on probabilistic methods for data science

rank. Compute the best rank r approximation19 of Y. The result, as we will show,
will be a good approximation to X.
But before we show this, let us define sampling of entries more rigorously.
Assume each entry of X is shown or hidden independently of others with fixed
probability p. Which entries are shown is decided by independent Bernoulli
random variables
m
δij ∼ Ber(p) with p := 2
n
which are often called selectors in this context. The value of p is chosen so that
among n2 entries of X, the expected number of selected (known) entries is m.
Define the n × n matrix Y with entries
Yij := δij Xij .
We can assume that we are shown Y, for it is a matrix that contains the observed
entries of X while all unobserved entries are replaced with zeros. The following
result shows how to estimate X based on Y.
Theorem 3.3.2 (Matrix completion). Let X̂ be a best rank r approximation to p−1 Y.
Then
r
1
rn
(3.3.3)
E kX̂ − XkF 6 C log(n)
kXk∞ ,
n
m
Here kXk∞ = maxi,j |Xij | denotes the maximum magnitude of the entries of X.
Before we prove this result, let us understand what this bound says about the
quality of matrix completion. The recovery error is measured in the Frobenius
norm, and the left side of (3.3.3) is
n
 1 X
1/2
1
2
|
X̂
−
X
|
.
kX̂ − XkF =
ij
ij
n
n2
i,j=1

Thus Theorem 3.3.2 controls the average error per entry in the mean-squared sense.
To make the error small, let us assume that we have a sample of size
m  rn log2 n,
which is slightly larger than the ideal size we discussed in (3.3.1). This makes
p
C log(n) rn/m = o(1) and forces the recovery error to be bounded by o(1)kXk∞ .
Summarizing, Theorem 3.3.2 says that the expected average error per entry is much
smaller than the maximal magnitude of the entries of X. This is true for a sample of
almost optimal size m. The smaller the rank r of the matrix X, the fewer entries
of X we need to see in order to do matrix completion.
Proof of Theorem 3.3.2. Step 1: The error in the operator norm. Let us first bound
the recovery error in the operator norm. Decompose the error into two parts using
19The best rank r approximation of an n × n matrix A is a matrix B of rank r that minimizes the

operator norm kA − Bk or, alternatively, the Frobenius norm kA − BkF (the minimizer
out to
Pturns
n
be the same). One can compute
B
by
truncating
the
singular
value
decomposition
A
=
s
ui vT
i
i=1
i
P
of A as follows: B = ri=1 si ui vT
,
where
we
assume
that
the
singular
values
s
are
arranged
in
i
i
non-increasing order.

Roman Vershynin

27

triangle inequality:
kX̂ − Xk 6 kX̂ − p−1 Yk + kp−1 Y − Xk.
Recall that X̂ is a best approximation to p−1 Y. Then the first part of the error is
smaller than the second part, i.e. kX̂ − p−1 Yk 6 kp−1 Y − Xk, and we have
2
(3.3.4)
kX̂ − Xk 6 2kp−1 Y − Xk = kY − pXk.
p
The entries of the matrix Y − pX,
(Y − pX)ij = (δij − p)Xij ,
are independent and mean zero random variables. Thus we can apply the bound
(3.2.7) on the norms of random matrices and get

(3.3.5)
E kY − pXk 6 C log n · E max k(Y − pX)i k2 + E max k(Y − pX)j k2 .
i∈[n]

j∈[n]

All that remains is to bound the norms of the rows and columns of Y − pX.
This is not difficult if we note that they can be expressed as sums of independent
random variables:
n
n
X
X
(δij − p)2 X2ij 6
(δij − p)2 · kXk2∞ ,
k(Y − pX)i k22 =
j=1

j=1

and similarly for columns. Taking expectation and noting that E(δij − p)2 =
Var(δij ) = p(1 − p), we get20
√
(3.3.6)
E k(Y − pX)i k2 6 (E k(Y − pX)i k22 )1/2 6 pn kXk∞ .
This is a good bound, but we need something stronger in (3.3.5). Since the maximum appears inside the expectation, we need a uniform bound, which will say
that all rows are bounded simultaneously with high probability.
Such uniform bounds are usually proved by applying concentration inequalities followed by a union bound. Bernstein’s inequality (1.4.5) yields


n
X

P
(δij − p)2 > tpn 6 exp(−ctpn) for t > 3.


j=1

(Check!) This probability can be further bounded by n−ct using the assumption
that m = pn2 > n log n. A union bound over n rows leads to


n


X
P max
(δij − p)2 > tpn 6 n · n−ct for t > 3.
i∈[n]

j=1

Integrating this tail, we conclude using (2.2.10) that
n
X
E max
(δij − p)2 . pn.
i∈[n]

j=1

20The first bound below that compares the L1 and L2 averages follows from Hölder’s inequality.

28

Four lectures on probabilistic methods for data science

(Check!) And this yields the desired bound on the rows,
√
E max k(Y − pX)i k2 . pn,
i∈[n]

which is an improvement of (3.3.6) we wanted. We can do similarly for the
columns. Substituting into (3.3.5), this gives
√
E kY − pXk . log(n) pn kXk∞ .
Then, by (3.3.4), we get
(3.3.7)

E kX̂ − Xk . log(n)

r

n
kXk∞ .
p

Step 2: Passing to Frobenius norm. Now we will need to pass from the
operator to Frobenius norm. This is where we will use for the first (and only)
time the rank of X. We know that rank(X) 6 r by assumption and rank(X̂) 6 r
by construction, so rank(X̂ − X) 6 2r. There is a simple relationship between the
operator and Frobenius norms:
√
kX̂ − XkF 6 2rkX̂ − Xk.
(Check it!) Take expectation of both sides and use (3.3.7); we get
r
√
rn
E kX̂ − XkF 6 2r E kX̂ − Xk . log(n)
kXk∞ .
p
Dividing both sides by n, we can rewrite this bound as
r
1
rn
E kX̂ − XkF . log(n)
kXk∞ .
n
pn2
But pn2 = m by definition of the sampling probability p. This yields the desired
bound (3.3.3).

3.4. Notes Theorem 3.1.2 on covariance estimation is a version of [58, Corollary 5.52], see also [36]. The logarithmic factor is in general necessary. This
theorem is a general-purpose result. If one knows some additional structural information about the covariance matrix (such as sparsity), then fewer samples may
be needed, see e.g. [12, 16, 40].
A version of Theorem 3.2.3 was proved in [51] in a more technical way. Although the logarithmic factor in Theorem 3.2.3 can not be completely removed in
general, it can be improved. Our argument actually gives
p
E kAk 6 C log n · E max kAi k2 + C log n · E max |Aij |,
i

ij

p
Using different methods, one can save an extra log n factor and show that
p
E kAk 6 C E max kAi k2 + C log n · E max |Aij |
i

ij

(see [6]) and
E kAk 6 C

p
log n · log log n · E max kAi k2 ,
i

Roman Vershynin

29

see [55]. (The results in [6, 55] are stated for Gaussian random matrices; the two
bounds above can be deduced by using conditioning and symmetrization.) The
surveys [6, 58] and the textbook [60] present several other useful techniques to
bound the operator norm of random matrices.
The matrix completion problem, which we discussed in Section 3.3, has attracted a lot of recent attention. E. Candes and B. Recht [14] showed that one can
often achieve exact matrix completion, thus computing the precise values of all
missing values of a matrix, from m ∼ rn log2 (n) randomly sampled entries. For
exact matrix completion, one needs an extra incoherence assumption that is not
present in Theorem 3.3.2. This assumption basically excludes matrices that are
simultaneously sparse and low rank (such as a matrix whose all but one entries
are zero – it would be extremely hard to complete it, since sampling will likely
miss the non-zero entry). Many further results on exact matrix completion are
known, e.g. [15, 18, 28, 56].
Theorem 3.3.2 with a simple proof is borrowed from [50]; see also the tutorial
[59]. This result only guarantees approximate matrix completion, but it does not
have any incoherence assumptions on the matrix.

4. Lecture 4: Matrix deviation inequality
In this last lecture we will study a new uniform deviation inequality for random matrices. This result will be a far reaching generalization of the JohnsonLindenstrauss Lemma we proved in Lecture 1.
Consider the same setup as in Theorem 1.6.1, where A is an m × n random matrix whose rows are independent, mean zero, isotropic and sub-gaussian random
vectors in Rn . (If you find it helpful to think in terms of concrete examples, let
the entries of A be independent N(0, 1) random variables.) Like in the JohnsonLindenstrauss Lemma, we will be looking at A as a linear transformation from
Rn to Rm , and we will be interested in what A does to points in some set in Rn .
This time, however, we will allow for infinite sets T ⊂ Rn .
Let us start by analyzing what A does to a single fixed vector x ∈ Rn . We have
E kAxk22 = E
=

m
X

j=1
m
X

Aj , x

E Aj , x

2

2

(where AT
j denote the rows of A)

(by linearity)

j=1

= mkxk22

(using isotropy of Aj ).

Further, if we assume that concentration about the mean holds here (and in fact,
it does), we should expect that
√
(4.0.1)
kAxk2 ≈ m kxk2
with high probability.

30

Four lectures on probabilistic methods for data science

Similarly to Johnson-Lindenstrauss Lemma, our next goal is to make (4.0.1)
hold simultaneously over all vectors x in some fixed set T ⊂ Rn . Precisely, we
may ask – how large is the average uniform deviation:
√
(4.0.2)
E sup kAxk2 − m kxk ?
x∈T

This quantity should clearly depend on some notion of the size of T : the larger
T , the larger should the uniform deviation be. So, how can we quantify the size
of T for this problem? In the next section we will do precisely this – introduce a
convenient, geometric measure of the sizes of sets in Rn , which is called Gaussian
width.
4.1. Gaussian width
Definition 4.1.1. Let T ⊂ Rn be a bounded set, and g be a standard normal random
vector in Rn , i.e. g ∼ N(0, In ). Then the quantities
w(T ) := E sup hg, xi

and γ(T ) := E sup | hg, xi |
x∈T

x∈T

are called the Gaussian width of T and the Gaussian complexity of T , respectively.
Gaussian width and Gaussian complexity are closely related. Indeed,21
(4.1.2)
1
1
1
1
w(T ) = w(T − T ) = E sup hg, x − yi = E sup | hg, x − yi | = γ(T − T ).
2
2 x,y∈T
2 x,y∈T
2
(Check these identities!)
Gaussian width has a natural geometric interpretation. Suppose g is a unit
vector in Rn . Then a moment’s thought reveals that supx,y∈T hg, x − yi is simply
the width of T in the direction of g, i.e. the distance between the two hyperplanes
with normal g that touch T on both sides as shown in Figure 4.1.3. Then 2w(T )
can be obtained by averaging the width of T over all directions g in Rn .

x
width
T
g
y
Figure 4.1.3. The width of a set T in the direction of g.
21The set T − T is defined as {x − y : x, y ∈ T }. More generally, given two sets A and B in the

same vector space, the Minkowski sum of A and B is defined as A + B = {a + b : a ∈ A, b ∈ B}.

Roman Vershynin

31

This reasoning is valid except where we assumed that g is a unit vector. Instead, for g ∼ N(0, In ) we have E kgk22 = n and
√
kgk2 ≈ n with high probability.
(Check both these claims using Bernstein’s inequality.) Thus, we need to scale
√
by the factor n. Ultimately, the geometric interpretation of the Gaussian width
√
becomes the following: w(T ) is approximately n/2 larger than the usual, geometric
width of T averaged over all directions.
A good exercise is to compute the Gaussian width and complexity for some
simple sets, such as the unit balls of the `p norms in Rn , which we denote Bn
p =
{x ∈ Rn : kxkp 6 1}. In particular, we have
p
√
γ(Bn
log n.
(4.1.4)
γ(Bn
2 ) ∼ n,
1)∼
For any finite set T ⊂ Bn
2 , we have
(4.1.5)

γ(T ) .

q

log |T |.

The same holds for Gaussian width w(T ). (Check these facts!)
A look a these examples reveals that the Gaussian width captures some nonobvious geometric qualities of sets. Of course, the fact that the Gaussian width of
√
n is not surprising: the usual, geometric
the unit Euclidean ball Bn
2 is or order
√
width in all directions is 2 and the Gaussian width is about n times that. But
it may be surprising that the Gaussian width of the `1 ball Bn
1 is much smaller,
and so is the width of any finite set T (unless the set has exponentially large
cardinality). As we will see later, Gaussian width nicely captures the geometric
size of “the bulk” of a set.
4.2. Matrix deviation inequality Now we are ready to answer the question we
asked in the beginning of this lecture: what is the magnitude of the uniform deviation (4.0.2)? The answer is surprisingly simple: it is bounded by the Gaussian
complexity of T . The proof is not too simple however, and we will skip it (see the
notes after this lecture for references).
Theorem 4.2.1 (Matrix deviation inequality). Let A be an m × n matrix whose rows
Ai are independent, isotropic and sub-gaussian random vectors in Rn . Let T ⊂ Rn be a
fixed bounded set. Then
√
E sup kAxk2 − mkxk2 6 CK2 γ(T )
x∈T

where K = maxi kAi kψ2 is the maximal sub-gaussian norm22 of the rows of A.
Remark 4.2.2 (Tail bound). It is often useful to have results that hold with high
probability rather than in expectation. There exists a high-probability version of
the matrix deviation inequality, and it states the following. Let u > 0. Then the
22A definition of the sub-gaussian norm of a random vector was given in Section 1.5. For example, if

A is a Gaussian random matrix with independent N(0, 1) entries, then K is an absolute constant.

32

Four lectures on probabilistic methods for data science

event
(4.2.3)

sup kAxk2 −

√
mkxk2 6 CK2 [γ(T ) + u · rad(T )]

x∈T

holds with probability at least 1 − 2 exp(−u2 ). Here rad(T ) is the radius of T ,
defined as
rad(T ) := sup kxk2 .
x∈T

Since rad(T ) . γ(T ) (check!) we can continue the bound (4.2.3) by
. K2 uγ(T )
for all u > 1. This is a weaker but still a useful inequality. For example, we can
use it to bound all higher moments of the deviation:

p 1/p
√
(4.2.4)
E sup kAxk2 − mkxk2
6 Cp K2 γ(T )
x∈T

√
where Cp 6 C p for p > 1. (Check this using Proposition 1.1.1.)
Remark 4.2.5 (Deviation of squares). It is sometimes helpful to bound the deviation of the square kAxk22 rather than kAxk2 itself. We can easily deduce the deviation of squares by using the identity a2 − b2 = (a − b)2 + 2b(a − b) for a = kAxk2
√
and b = mkxk2 . Doing this, we conclude that
√
(4.2.6)
E sup kAxk22 − mkxk22 6 CK4 γ(T )2 + CK2 m rad(T )γ(T ).
x∈T

(Do this calculation using (4.2.4) for p = 2.) We will use this bound in Section 4.4.
Matrix deviation inequality has many consequences. We will explore some of
them now.
4.3. Deriving Johnson-Lindenstrauss Lemma We started this lecture by promising a result that is more general than Johnson-Lindenstrauss Lemma. So let us
show how to quickly derive Johnson-Lindenstrauss from the matrix deviation
inequality. Theorem 1.6.1 from Theorem 4.2.1.
Assume we are in the situation of the Johnson-Lindenstrauss Lemma (Theorem 1.6.1). Given a set X ⊂ R, consider the normalized difference set
x−y
: x, y ∈ X distinct vectors .
T :=
kx − yk2
Then T is a finite subset of the unit sphere of Rn , and thus (4.1.5) gives
q
q
q
γ(T ) . log |T | 6 log |X|2 . log |X|.
Matrix deviation inequality (Theorem 4.2.1) then yields
p
√
kA(x − y)k2 √
− m . log N 6 ε m.
sup
kx − yk2
x,y∈X
with high probability, say 0.99. (To pass from expectation to high probability, we
can use Markov’s inequality. To get the last bound, we use the assumption on m
in Johnson-Lindenstrauss Lemma.)

Roman Vershynin

33

√
Multiplying both sides by kx − yk2 / m, we can write the last bound as follows.
With probability at least 0.99, we have
1
(1 − ε)kx − yk2 6 √ kAx − Ayk2 6 (1 + ε)kx − yk2 for all x, y ∈ X.
m
This is exactly the consequence of Johnson-Lindenstrauss lemma.
The argument based on matrix deviation inequality, which we just gave, can
be easily extended for infinite sets. It allows one to state a version of JohnsonLindenstrauss lemma for general, possibly infinite, sets, which depends on the
Gaussian complexity of T rather than cardinality. (Try to do this!)
4.4. Covariance estimation In Section 3.1, we introduced the problem of covariance estimation, and we showed that
N ∼ n log n
samples are enough to estimate the covariance matrix of a general distribution
in Rn . We will now show how to do better if the distribution is sub-gaussian.
(Recall Section 1.5 for the definition of sub-gaussian random vectors.) In this case,
we can get rid of the logarithmic oversampling and the boundedness condition
(3.1.3).
Theorem 4.4.1 (Covariance estimation for sub-gaussian distributions). Let X be a
random vector in Rn with covariance matrix Σ. Suppose X is sub-gaussian, and more
specifically
(4.4.2)

k hX, xi kψ2 . k hX, xi kL2 = kΣ1/2 xk2

for any x ∈ Rn .

Then, for every N > 1, we have
r n

n
.
N N
This result implies that if, for ε ∈ (0, 1), we take a sample of size
E kΣN − Σk . kΣk

+

N ∼ ε−2 n,
then we are guaranteed covariance estimation with a good relative error:
E kΣN − Σk 6 εkΣk.
Proof. Since we are going to use Theorem 4.2.1, we will need to first bring the
random vectors X, X1 , . . . , XN to the isotropic position. This can be done by a
suitable linear transformation. You will easily check that there exists an isotropic
random vector Z such that
X = Σ1/2 Z.
(For example, Σ has full rank, set Z := Σ−1/2 X. Check the general case.) Similarly,
we can find independent and isotropic random vectors Zi such that
Xi = Σ1/2 Zi ,

i = 1, . . . , N.

34

Four lectures on probabilistic methods for data science

The sub-gaussian assumption (4.4.2) then implies that
kZkψ2 . 1.
(Check!) Then
1 X
Zi ZT
i − In .
N
N

kΣN − Σk = kΣ1/2 RN Σ1/2 k

where

RN :=

i=1

The operator norm of a symmetric n × n matrix A can be computed by maximizing the quadratic form over the unit sphere: kAk = maxx∈Sn−1 | hAx, xi |. (To
see this, recall that the operator norm is the biggest eigenvalue of A in magnitude.) Then
D
E
kΣN − Σk = max Σ1/2 RN Σ1/2 x, x = max hRN x, xi
x∈T

x∈Sn−1

where T is the ellipsoid
T := Σ1/2 Sn−1 .
Recalling the definition of RN , we can rewrite this as
1
1 X
hZi , xi2 − kxk22 =
max kAxk22 − Nkxk22 .
N
N x∈T
N

kΣN − Σk = max
x∈T

i=1

Now we apply the matrix deviation inequality for squares (4.2.6) and conclude
that

√
1
kΣN − Σk .
γ(T )2 + N rad(T )γ(T ) .
N
(Do this calculation!) The radius and Gaussian width of the ellipsoid T are easy
to compute:
rad(T ) = kΣk1/2 and γ(T ) 6 (tr Σ)1/2 .
Substituting, we get

p
1
tr Σ + NkΣk tr Σ .
N
To complete the proof, use that tr Σ 6 nkΣk (check this!) and simplify the bound.

kΣN − Σk .

Remark 4.4.3 (Low-dimensional distributions). Similarly to Section 3.1, we can
show that much fewer samples are needed for covariance estimation of lowdimensional sub-gaussian distributions. Indeed, the proof of Theorem 4.4.1 actually yields
r r
r
(4.4.4)
E kΣN − Σk 6 kΣk
+
N N
where
tr Σ
r = r(Σ1/2 ) =
kΣk
is the stable rank of Σ1/2 . This means that covariance estimation is possible with
N∼r

Roman Vershynin

35

samples.
4.5. Underdetermined linear equations We will give one more application of
the matrix deviation inequality – this time, to the area of high dimensional inference. Suppose we need to solve a severely underdetermined system of linear
equations: say, we have m equations in n  m variables. Let us write it in the
matrix form as
y = Ax
where A is a given m × n matrix, y ∈ Rm is a given vector and x ∈ Rn is an
unknown vector. We would like to compute x from A and y.
When the linear system is underdetermined, we can not find x with any accuracy, unless we know something extra about x. So, let us assume that we do
have some a-priori information. We can describe this situation mathematically by
assuming that
x∈K
where K ⊂ Rn is some known set in Rn that describes anything that we know
about x a-priori. (Admittedly, we are operating on a high level of generality here.
If you need a concrete example, we will consider it in Section 4.6.)
Summarizing, here is the problem we are trying to solve. Determine a solution
x = x(A, y, K) to the underdetermined linear equation y = Ax as accurately as
possible, assuming that x ∈ K.
A variety of approaches to this and similar problems were proposed during the
last decade; see the notes after this lecture for pointers to some literature. The one
we will describe here is based on optimization. To do this, it will be convenient
to convert the set K into a function on Rn which is called the Minkowski functional
of K. This is basically a function whose level sets are multiples of K. To define it
formally, assume that K is star-shaped, which means that together with any point
x, the set K must contain the entire interval that connects x with the origin; see
Figure 4.5.1 for illustration. The Minkowski functional of K is defined as

kxkK := inf t > 0 : x/t ∈ K , x ∈ Rn .
If the set K is convex and symmetric about the origin, kxkK is actually a norm on
Rn . (Check this!)

0

0

Figure 4.5.1. The set on the left (whose boundary is shown) is starshaped, the set on the right is not.

36

Four lectures on probabilistic methods for data science

Now we propose the following way to solve the recovery problem: solve the
optimization program
min kx 0 kK

(4.5.2)

subject to

y = Ax 0 .

Note that this is a very natural program: it looks at all solutions to the equation
y = Ax 0 and tries to “shrink” the solution x 0 toward K. (This is what minimization of Minkowski functional is about.)
Also note that if K is convex, this is a convex optimization program, and thus
can be solved effectively by one of the many available numeric algorithms.
The main question we should now be asking is – would the solution to this
program approximate the original vector x? The following result bounds the
approximation error for a probabilistic model of linear equations. Assume that A
is a random matrix as in Theorem 4.2.1, i.e. A is an m × n matrix whose rows Ai
are independent, isotropic and sub-gaussian random vectors in Rn .
Theorem 4.5.3 (Recovery by optimization). The solution x̂ of the optimization program (4.5.2) satisfies23
w(K)
E kx̂ − xk2 . √ ,
m
where w(K) is the Gaussian width of K.
Proof. Both the original vector x and the solution x̂ are feasible vectors for the
optimization program (4.5.2). Then
kx̂kK 6 kxkK
61

(since x̂ minimizes the Minkowski functional)

(since x ∈ K).

Thus both x̂, x ∈ K.
We also know that Ax̂ = Ax = y, which yields
(4.5.4)

A(x̂ − x) = 0.

Let us apply matrix deviation inequality (Theorem 4.2.1) for T := K − K. It
gives
√
E sup kA(u − v)k2 − m ku − vk2 . γ(K − K) = 2w(K),
u,v∈K

where we used (4.1.2) in the last identity. Substitute u = x̂ and v = x here. We
may do this since, as we noted above, both these vectors belong to K. But then
the term kA(u − v)k2 will equal zero by (4.5.4). It disappears from the bound, and
we get
√
E m kx̂ − xk2 . w(K).
√
Dividing both sides by m we complete the proof.

23Here and in other similar results, the notation . will hide possible dependence on the sub-gaussian

norms of the rows of A.

Roman Vershynin

37

Theorem 4.5.3 says that a signal x ∈ K can be efficiently recovered from
m ∼ w(K)2
random linear measurements.
4.6. Sparse recovery Let us illustrate Theorem 4.5.3 with an important specific
example of the feasible set K.
Suppose we know that the signal x is sparse, which means that only a few
coordinates of x are nonzero. As before, our task is to recover x from the random
linear measurements given by the vector
y = Ax,
where A is an m × n random matrix. This is a basic example of sparse recovery
problems, which are ubiquitous in various disciplines.
The number of nonzero coefficients of a vector x ∈ Rn , or the sparsity of x,
is often denoted kxk0 . This is similar to the notation for the `p norm kxkp =
P
p 1/p , and for a reason. You can quickly check that
( n
i=1 |xi | )
(4.6.1)

kxk0 = lim kxkp
p→0

(Do this!) Keep in mind that neither kxk0 nor kxkp for 0 < p < 1 are actually
norms on Rn , since they fail triangle inequality. (Give an example.)
Let us go back to the sparse recovery problem. Our first attempt to recover x
is to try the following optimization problem:
(4.6.2)

min kx 0 k0

subject to y = Ax 0 .

This is sensible because this program selects the sparsest feasible solution. But
there is an implementation caveat: the function f(x) = kxk0 is highly non-convex
and even discontinuous. There is simply no known algorithm to solve the optimization problem (4.6.2) efficiently.
To overcome this difficulty, let us turn to the relation (4.6.1) for an inspiration.
What if we replace kxk0 in the optimization problem (4.6.2) by kxkp with p > 0?
The smallest p for which f(x) = kxkp is a genuine norm (and thus a convex
function on Rn ) is p = 1. So let us try
(4.6.3)

min kx 0 k1

subject to

y = Ax 0 .

This is a convexification of the non-convex program (4.6.2), and a variety of numeric convex optimization methods are available to solve it efficiently.
We will now show that `1 minimization works nicely for sparse recovery. As
before, we assume that A is a random matrix as in Theorem 4.2.1.
Theorem 4.6.4 (Sparse recovery by optimization). Assume that an unknown vector
x ∈ Rn has at most s non-zero coordinates, i.e. kxk0 6 s. The solution x̂ of the

38

Four lectures on probabilistic methods for data science

optimization program (4.6.3) satisfies
r

s log n
kxk2 .
m
Proof. Since kxk0 6 s, Cauchy-Schwarz inequality shows that
√
(4.6.5)
kxk1 6 s kxk2 .
E kx̂ − xk2 .

n
n
(Check!) Denote the unit ball of the `1 norm in Rn by Bn
1 , i.e. B1 := {x ∈ R :
kxk1 6 1}. Then we can rewrite (4.6.5) as the inclusion
√
x ∈ s kxk2 · Bn
1 := K.

Apply Theorem 4.5.3 for this set K. We noted the Gaussian width of Bn
1 in (4.1.4),
so
p
√
√
√
n
w(K) = s kxk2 · w(Bn
log n.
1 ) 6 s kxk2 · γ(B1 ) 6 s kxk2 ·
Substitute this in Theorem 4.5.3 and complete the proof.



Theorem 4.6.4 says that an s-sparse signal x ∈ Rn can be efficiently recovered
from
m ∼ s log n
random linear measurements.
4.7. Notes For a more thorough introduction to Gaussian width and its role in
high-dimensional estimation, refer to the tutorial [59] and the textbook [60]; see
also [5]. Related to Gaussian complexity is the notion of Rademacher complexity
of T , obtained by replacing the coordinates of g by independent Rademacher (i.e.
±1 symmetric) random variables. Rademacher complexity of classes of functions
plays an important role in statistical learning theory, see e.g. [44]
Matrix deviation inequality (Theorem 4.2.1) is borrowed from [41]. In the special case where A is a Gaussian random matrix, this result follows from the work
of G. Schechtman [57] and could be traced back to results of Gordon [24–27].
In the general case of sub-gaussian distributions, earlier variants of Theorem 4.2.1 were proved by B. Klartag and S. Mendelson [35], S. Mendelson, A. Pajor
and N. Tomczak-Jaegermann [45] and S. Dirksen [20].
Theorem 4.4.1 for covariance estimation can be proved alternatively using more
elementary tools (Bernstein’s inequality and ε-nets), see [58]. However, no known
elementary approach exists for the low-rank covariance estimation discussed in
Remark 4.4.3. The bound (4.4.4) was proved by V. Koltchinskii and K. Lounici
[36] by a different method.
In Section 4.5, we scratched the surface of a recently developed area of sparse
signal recovery, which is also called compressed sensing. Our presentation there
essentially follows the tutorial [59]. Theorem 4.6.4 can be improved: if we take
m & s log(n/s)

Roman Vershynin

39

measurements, then with high probability the optimization program (4.6.3) recovers the unknown signal x exactly, i.e.
x̂ = x.
First results of this kind were proved by J. Romberg, E. Candes and T. Tao [13]
and a great number of further developments followed; refer e.g. to the book [23]
and the chapter in [19] for an introduction into this research area.

Acknowledgement
I am grateful to the referees who made a number of useful suggestions, which
led to better presentation of the material in this chapter.

References
[1] E. Abbe, A. S. Bandeira, G. Hall. Exact recovery in the stochastic block model, IEEE Transactions on
Information Theory 62 (2016), 471–487. 19
[2] D. Achlioptas, Database-friendly random projections: Johnson-Lindenstrauss with binary coins, Journal
of Computer and System Sciences, 66 (2003), 671–687. 9
[3] R. Ahlswede, A. Winter, Strong converse for identification via quantum channels, IEEE Trans. Inf.
Theory 48 (2002), 569–579. 18
[4] N. Ailon, B. Chazelle, Approximate nearest neighbors and the fast Johnson-Lindenstrauss transform,
Proceedings of the 38th Annual ACM Symposium on Theory of Computing. New York: ACM
Press, 2006. pp. 557–563. 9
[5] D. Amelunxen, M. Lotz, M. McCoy, J. Tropp, Living on the edge: phase transitions in convex programs
with random data, Inf. Inference 3 (2014), 224–294. 38
[6] A. Bandeira, R. van Handel, Sharp nonasymptotic bounds on the norm of random matrices with independent entries, Ann. Probab. 44 (2016), 2479–2506. 19, 28, 29
[7] R. Baraniuk, M. Davenport, R. DeVore, M. Wakin, A simple proof of the restricted isometry property
for random matrices, Constructive Approximation, 28 (2008), 253–263. 9
[8] R. Bhatia, Matrix Analysis. Graduate Texts in Mathematics, vol. 169. Springer, Berlin, 1997. 18
[9] C. Bordenave, M. Lelarge, L. Massoulie, |em Non-backtracking spectrum of random graphs:
community detection and non-regular Ramanujan graphs, Annals of Probability, to appear. 19
[10] S. Boucheron, G. Lugosi, P. Massart, Concentration inequalities. A nonasymptotic theory of independence. With a foreword by Michel Ledoux. Oxford University Press, Oxford, 2013. 9
[11] O. Bousquet1, S. Boucheron, G. Lugosi, Introduction to statistical learning theory, in: Advanced
Lectures on Machine Learning, Lecture Notes in Computer Science 3176, pp.169–207, Springer
Verlag 2004.
[12] T. Cai, R. Zhao, H. Zhou, Estimating structured high-dimensional covariance and precision matrices:
optimal rates and adaptive estimation, Electron. J. Stat. 10 (2016), 1–59. 28
[13] E. Candes, J. Romberg, T. Tao, Robust uncertainty principles: exact signal reconstruction from highly
incomplete frequency information, IEEE Trans. Inform. Theory 52 (2006), 489–509. 39
[14] E. Candes, B. Recht, Exact Matrix Completion via Convex Optimization, Foundations of Computational Mathematics 9 (2009), 717–772. 29
[15] E. Candes, T. Tao, The power of convex relaxation: near-optimal matrix completion, IEEE Trans. Inform.
Theory 56 (2010), 2053–2080. 29
[16] R. Chen, A. Gittens, J. Tropp, The masked sample covariance estimator: an analysis using matrix concentration inequalities, Inf. Inference 1 (2012), 2–20. 28
[17] P. Chin, A. Rao, and V. Vu, |em Stochastic block model and community detection in the sparse
graphs: A spectral algorithm with optimal rate of recovery, preprint, 2015. 19
[18] M. Davenport, Y. Plan, E. van den Berg, M. Wootters, 1-bit matrix completion, Inf. Inference 3
(2014), 189–223. 29

40

References

[19] M. Davenport, M. Duarte, Yonina C. Eldar, Gitta Kutyniok, Introduction to compressed sensing, in:
Compressed sensing. Edited by Yonina C. Eldar and Gitta Kutyniok. Cambridge University Press,
Cambridge, 2012. 39
[20] S. Dirksen, Tail bounds via generic chaining, Electron. J. Probab. 20 (2015), 1–29. 38
[21] U. Feige, E. Ofek, Spectral techniques applied to sparse random graphs, Random Structures Algorithms 27 (2005), 251–275. 19
[22] S. Fortunato, Santo; D. Hric, Community detection in networks: A user guide. Phys. Rep. 659 (2016),
1–44. 19
[23] S. Foucart, H. Rauhut, A mathematical introduction to compressive sensing. Applied and Numerical
Harmonic Analysis. Birkhäuser/Springer, New York, 2013. 39
[24] Y. Gordon, Some inequalities for Gaussian processes and applications, Israel J. Math. 50 (1985), 265–289.
38
[25] Y. Gordon, Elliptically contoured distributions, Prob. Th. Rel. Fields 76 (1987), 429–438. 38
[26] Y. Gordon, On Milman’s inequality and random subspaces which escape through a mesh in Rn , Geometric aspects of functional analysis (1986/87), Lecture Notes in Math., vol. 1317, pp. 84–106.
38
[27] Y. Gordon, Majorization of Gaussian processes and geometric applications, Prob. Th. Rel. Fields 91
(1992), 251–267. 38
[28] D. Gross, Recovering low-rank matrices from few coefficients in any basis, IEEE Trans. Inform. Theory
57 (2011), 1548–1566. 29
[29] O. Guedon, R. Vershynin, Community detection in sparse networks via Grothendieck’s inequality, Probability Theory and Related Fields 165 (2016), 1025–1049. 19
[30] A. Javanmard, A. Montanari, F. Ricci-Tersenghi, Phase transitions in semidefinite relaxations, PNAS,
April 19, 2016, vol. 113, no.16, E2218–E2223. 19
[31] W. B. Johnson, J. Lindenstrauss, Extensions of Lipschitz mappings into a Hilbert space. In Beals,
Richard; Beck, Anatole; Bellow, Alexandra; et al. Conference in modern analysis and probability
(New Haven, Conn., 1982). Contemporary Mathematics. 26. Providence, RI: American Mathematical Society, 1984. pp. 189–206. 9
[32] B. Hajek, Y. Wu, J. Xu, Achieving exact cluster recovery threshold via semidefinite programming, IEEE
Transactions on Information Theory 62 (2016), 2788–2797. 19
[33] P. W. Holland, K. B. Laskey, S. Leinhardt, |em Stochastic blockmodels: first steps, Social Networks
5 (1983), 109–137. 19
[34] D. Kane, J. Nelson, Sparser Johnson-Lindenstrauss Transforms, Journal of the ACM 61 (2014): 1. 9
[35] B. Klartag, S. Mendelson, Empirical processes and random projections, J. Funct. Anal. 225 (2005),
229–245. 38
[36] V. Koltchinskii, K. Lounici, Concentration inequalities and moment bounds for sample covariance operators, Bernoulli 23 (2017), 110–133. 28, 38
[37] C. Le, E. Levina, R. Vershynin, Concentration and regularization of random graphs, Random Structures and Algorithms, to appear. 19
[38] M. Ledoux, The concentration of measure phenomenon. American Mathematical Society, Providence,
RI, 2001. 9
[39] M. Ledoux, M. Talagrand, Probability in Banach spaces. Isoperimetry and processes. Springer-Verlag,
Berlin, 1991. 9
[40] E. Levina, R. Vershynin, Partial estimation of covariance matrices, Probability Theory and Related
Fields 153 (2012), 405–419. 28
[41] C. Liaw, A. Mehrabian, Y. Plan, R. Vershynin, A simple tool for bounding the deviation of random matrices on geometric sets, Geometric Aspects of Functional Analysis, Lecture Notes in Mathematics,
Springer, Berlin, to appear. 38
[42] J. Matouĺek, Lectures on discrete geometry. Graduate Texts in Mathematics, 212. Springer-Verlag,
New York, 2002. 9
[43] F. McSherry, Spectral partitioning of random graphs, Proc. 42nd FOCS (2001), 529–537. 19
[44] S.Mendelson,S.Mendelson, A few notes on statistical learning theory, in: Advanced Lectures on Machine Learning, S. Mendelson, A.J. Smola (Eds.) LNAI 2600, pp. 1–40, 2003. 38
[45] S. Mendelson, A. Pajor, N. Tomczak-Jaegermann, Reconstruction and subgaussian operators in asymptotic geometric analysis. Geom. Funct. Anal. 17 (2007), 1248–1282. 38
[46] E. Mossel, J. Neeman, A. Sly, Belief propagation, robust reconstruction and optimal recovery of block
models. Ann. Appl. Probab. 26 (2016), 2211–2256. 19

References

41

[47] M. E. Newman, Networks. An introduction. Oxford University Press, Oxford, 2010. 19
[48] R. I. Oliveira, Concentration of the adjacency matrix and of the Laplacian in random graphs with independent edges, unpublished (2010), arXiv:0911.0600. 18, 19
[49] R. I. Oliveira, Sums of random Hermitian matrices and an inequality by Rudelson, Electron. Commun.
Probab. 15 (2010), 203–212. 18
[50] Y. Plan, R. Vershynin, E. Yudovina, High-dimensional estimation with geometric constraints, Information and Inference 0 (2016), 1–40. 29
[51] S. Riemer, C. Schütt, On the expectation of the norm of random matrices with non-identically distributed
entries, Electron. J. Probab. 18 (2013), no. 29, 13 pp. 28
[52] J. Tropp, User-friendly tail bounds for sums of random matrices. Found. Comput. Math. 12 (2012),
389–434. 18
[53] J. Tropp, An introduction to matrix concentration inequalities. Found. Trends Mach. Learning 8 (2015),
10-230. 18
[54] R. van Handel, Structured random matrices. in: IMA Volume “Discrete Structures: Analysis and
Applications”, Springer, to appear.
[55] R. van Handel, On the spectral norm of Gaussian random matrices, Trans. Amer. Math. Soc., to appear.
29
[56] B. Recht, A simpler approach to matrix completion, J. Mach. Learn. Res. 12 (2011), 3413–3430. 29
[57] G. Schechtman, Two observations regarding embedding subsets of Euclidean spaces in normed spaces,
Adv. Math. 200 (2006), 125–135. 38
[58] R. Vershynin, Introduction to the non-asymptotic analysis of random matrices. Compressed sensing,
210–268, Cambridge University Press, Cambridge, 2012. 2, 9, 28, 29, 38
[59] R. Vershynin, Estimation in high dimensions: a geometric perspective. Sampling Theory, a Renaissance,
3–66, Birkhauser Basel, 2015. 2, 29, 38
[60] R. Vershynin, High-Dimensional Probability. An Introduction with Applications in Data Science. Cambridge University Press, to appear. 2, 9, 29, 38
[61] H. Zhou, A. Zhang, Minimax Rates of Community Detection in Stochastic Block Models, Annals of
Statistics, to appear. 19
Department of Mathematics, University of Michigan, 530 Church Street, Ann Arbor, MI 48109,
U.S.A.
E-mail address: romanv@umich.edu

