An Introduction to Matrix Concentration Inequalities

Chapter 0 Preface

In recent years, random matrices have come to play a major role in computational mathematics, but most of the classical areas of random matrix theory remain the province of experts. Over the last decade, with the advent of matrix concentration inequalities, research has advanced to the point where we can conquer many (formerly) challenging problems with a page or two of arithmetic.

My aim is to describe the most successful methods from this area along with some interesting examples that these techniques can illuminate. I hope that the results in these pages will inspire future work on applications of random matrices as well as refinements of the matrix concentration inequalities discussed herein.

I have chosen to present a coherent body of results based on a generalization of the Laplace transform method for establishing scalar concentration inequalities. In the last two years, Lester Mackey and I, together with our coauthors, have developed an alternative approach to matrix concentration using exchangeable pairs and Markov chain couplings. With some regret, I have chosen to omit this theory because the ideas seem less accessible to a broad audience of researchers. The interested reader will find pointers to these articles in the annotated bibliography.

The work described in these notes reflects the influence of many researchers. These include Rudolf Ahlswede, Rajendra Bhatia, Eric Carlen, Sourav Chatterjee, Edward Effros, Elliott Lieb, Roberto Imbuzeiro Oliveira, Dénes Petz, Gilles Pisier, Mark Rudelson, Roman Vershynin, and Andreas Winter. I have also learned a great deal from other colleagues and friends along the way.

I would like to thank some people who have helped me improve this work. Several readers informed me about errors in the initial version of this manuscript; these include Serg Bogdanov, Peter Forrester, Nikos Karampatziakis, and Guido Lagos. The anonymous reviewers tendered many useful suggestions, and they pointed out a number of errors. Sid Barman gave me feedback on the final revisions to the monograph. Last, I want to thank Léon Nijensohn for his continuing encouragement.

I gratefully acknowledge financial support from the Office of Naval Research under awards N00014-08-1-0883 and N00014-11-1002, the Air Force Office of Strategic Research under award FA9550-09-1-0643, and an Alfred P. Sloan Fellowship. Some of this research was completed at the Institute of Pure and Applied Mathematics at UCLA. I would also like to thank the California Institute of Technology and the Moore Foundation.

Joel A. Tropp Pasadena, CA December 2012 Revised, March 2014 and December 2014

Chapter 1 Introduction

Random matrix theory has grown into a vital area of probability, and it has found applications in many other fields. To motivate the results in this monograph, we begin with an overview of the connections between random matrix theory and computational mathematics. We introduce the basic ideas underlying our approach, and we state one of our main results on the behavior of random matrices. As an application, we examine the properties of the sample covariance estimator, a random matrix that arises in statistics. Afterward, we summarize the other types of results that appear in these notes, and we assess the novelties in this presentation.

Random matrix theory sprang from several different sources in the first half of the 20th century.

Peter Forrester [For10, p. v] traces the field of random matrix theory to work of Hurwitz, who defined the invariant integral over a Lie group. Specializing this analysis to the orthogonal group, we can reinterpret this integral as the expectation of a function of a uniformly random orthogonal matrix.

Another early example of a random matrix appeared in the work of John Wishart [Wis28]. Wishart was studying the behavior of the sample covariance estimator for the covariance matrix of a multivariate normal random vector. He showed that the estimator, which is a random matrix, has the distribution that now bears his name. Statisticians have often used random matrices as models for multivariate data [MKB79, Mui82].

In their remarkable work [vNG47, GvN51] on computational methods for solving systems of linear equations, von Neumann and Goldstine considered a random matrix model for the floating-point errors that arise from an LU decomposition. von Neumann and Goldstine invented and analyzed this algorithm before they had any digital computer on which to implement it! See [Grc11] for a historical account. They obtained a high-probability bound for the norm of the random matrix, which they took as an estimate for the error the procedure might typically incur. Curiously, in subsequent years, numerical linear algebraists became very suspicious of probabilistic techniques, and only in recent years have randomized algorithms reappeared in this field. See the surveys [Mah11, HMT11, Woo14] for more details and references.

In the early 1950s, physicists had reached the limits of deterministic analytical techniques for studying the energy spectra of heavy atoms undergoing slow nuclear reactions. Eugene Wigner was the first researcher to surmise that a random matrix with appropriate symmetries might serve as a suitable model for the Hamiltonian of the quantum mechanical system that describes the reaction. The eigenvalues of this random matrix model the possible energy levels of the system. See Mehta’s book [Meh04, §1.1] for an account of all this.

In each area, the motivation was quite different and led to distinct sets of questions. Later, random matrices began to percolate into other fields such as graph theory (the Erdős–Rényi model [ER60] for a random graph) and number theory (as a model for the spacing of zeros of the Riemann zeta function [Mon73]).

The Modern Random Matrix

By now, random matrices are ubiquitous. They arise throughout modern mathematics and statistics, as well as in many branches of science and engineering. Random matrices have several different purposes that we may wish to distinguish. They can be used within randomized computer algorithms; they serve as models for data and for physical phenomena; and they are subjects of mathematical inquiry. This section offers a taste of these applications. Note that the ideas and references here reflect the author’s interests, and they are far from comprehensive!

The striking mathematical properties of random matrices can be harnessed to develop algorithms for solving many different problems.

Random matrices can be used to develop fast algorithms for computing a truncated singular-value decomposition. In this application, we multiply a large input matrix by a smaller random matrix to extract information about the dominant singular vectors of the input matrix. The seed of this idea appears in [FKV98, DFK+99]. The survey [HMT11] explains how to implement this method in practice, while the two monographs [Mah11, Woo14] cover more theoretical aspects.

One way to accelerate spectral computations on large matrices is to replace the original matrix by a sparse proxy that has similar spectral properties. An elegant way to produce the sparse proxy is to zero out entries of the original matrix at random while rescaling the entries that remain. This approach was proposed in [AM01, AM07], and the papers [AKL13, KD14] contain recent innovations. Related ideas play an important role in Spielman and Teng’s work [ST04] on fast algorithms for solving linear systems.

In large-scale machine learning, one may need to subsample data randomly to reduce the computational costs of fitting a model. For instance, we can combine random sampling with the Nyström decomposition to obtain a randomized approximation of a kernel matrix. This method was introduced by Williams & Seeger [WS01]. The paper [DM05] provides the first theoretical analysis, and the survey [GM14] contains more complete results.

A basic template in the theory of algorithms invokes randomized projection to reduce the dimension of a computational problem. Many types of dimension reduction are based on properties of random matrices. The two papers [JL84, Bou85] established the mathematical foundations of this approach. The earliest applications in computer science appear in the work [LLR95]. Many contemporary variants depend on ideas from [AC09] and [CW13].

One approach to solving a computationally difficult optimization problem is to relax (i.e., enlarge) the constraint set so the problem becomes tractable, to solve the relaxed problem, and then to use a randomized procedure to map the solution back to the original constraint set [BTN01, §4.3]. This technique is called relaxation and rounding. For hard optimization problems involving a matrix variable, the analysis of the rounding procedure often involves ideas from random matrix theory [So09, NRV13].

When acquiring data about an object with relatively few degrees of freedom as compared with the ambient dimension, we may be able to sieve out the important information from the object by taking a small number of random measurements, where the number of measurements is comparable to the number of degrees of freedom [GGI+02, CRT06, Don06]. This observation is now referred to as compressed sensing. Random matrices play a central role in the design and analysis of measurement procedures. For example, see [FR13, CRPW12, ALMT14, Tro14].

Modeling

Random matrices also appear as models for multivariate data or multivariate phenomena. By studying the properties of these models, we may hope to understand the typical behavior of a data-analysis algorithm or a physical system.

Sparse approximation has become an important problem in statistics, signal processing, machine learning and other areas. One model for a “typical” sparse signal poses the assumption that the nonzero coefficients that generate the signal are chosen at random. When analyzing methods for identifying the sparse set of coefficients, we must study the behavior of a random column submatrix drawn from the model matrix [Tro08a, Tro08b].

In data analysis, it is common to encounter a mixture of two structured signals, and the goal is to extract the two signals using prior information about the structures. A common model for this problem assumes that the signals are randomly oriented with respect to each other, which means that it is usually possible to discriminate the underlying structures. Random orthogonal matrices arise in the analysis of estimation techniques for this problem [MT14, ALMT14, MT13].

One probabilistic framework for describing community structure in a network assumes that each pair of individuals in the same community has a relationship with high probability, while each pair of individuals drawn from different communities has a relationship with lower probability. This is referred to as the stochastic block model [HLL83]. It is quite common to analyze algorithms for extracting community structure from data by positing that this model holds. See [ABH14] for a recent contribution, as well as a summary of the extensive literature.

More generally, random models are pervasive in the analysis of statistical estimation procedures for high-dimensional data. Random matrix theory plays a key role in this field [MKB79, Mui82, Kol11, BvdG11].

Random matrices are commonly used as models for wireless channels. See the book of Tulino and Verdú for more information [TV04].

In these examples, it is important to recognize that random models may not coincide very well with reality, but they allow us to get a sense of what might be possible in some generic cases.

Theoretical Aspects

Random matrices are frequently studied for their intrinsic mathematical interest. In some fields, they provide examples of striking phenomena. In other areas, they furnish counterexamples to “intuitive” conjectures. Here are a few disparate problems where random matrices play a role.

An expander graph has the property that every small set of vertices has edges linking it to a large proportion of the vertices. The expansion property is closely related to the spectral behavior of the adjacency matrix of the graph. The easiest construction of an expander involves a random matrix argument [AS00, §9.2].

For worst-case examples, the Gaussian elimination method for solving a linear system is not numerically stable. In practice, however, stability problems rarely arise. One explanation for this phenomenon is that, with high probability, a small random perturbation of any fixed matrix is well conditioned. As a consequence, it can be shown that Gaussian elimination is stable for most matrices [SST06].

Dvoretzky’s Theorem states that, when $N$ is large, the unit ball of each $N$-dimensional Banach space has a slice of dimension $n\approx\log N$ that is close to a Euclidean ball with dimension $n$. It turns out that a random slice of dimension $n$ realizes this property [Mil71]. This result can be framed as a statement about spectral properties of a random matrix [Gor85].

Random matrices appear as counterexamples for a number of conjectures in quantum information theory. Here is one instance. In classical information theory, the total amount of information that we can transmit through a pair of channels equals the sum of the information we can send through each channel separately. It was conjectured that the same property holds for quantum channels. In fact, a pair of quantum channels can have strictly larger capacity than a single channel. This result depends on a random matrix construction [Has09]. See [HW08] for related work.

Random Matrices for the People

Historically, random matrix theory has been regarded as a very challenging field. Even now, many well-established methods are only comprehensible to researchers with significant experience, and it may take months of intensive effort to prove new results. There are a small number of classes of random matrices that have been studied so completely that we know almost everything about them. Yet, moving beyond this terra firma, one quickly encounters examples where classical methods are brittle.

We hope to democratize random matrix theory. These notes describe tools that deliver useful information about a wide range of random matrices. In many cases, a modest amount of straightforward arithmetic leads to strong results. The methods here should be accessible to computational scientists working in a variety of fields. Indeed, the techniques in this work have already found an extensive number of applications.

Basic Questions in Random Matrix Theory

Random matrices merit special attention because they have spectral properties that are quite different from familiar deterministic matrices. Here are some of the questions we might want to investigate.

What is the expectation of the maximum eigenvalue of a random Hermitian matrix? What about the minimum eigenvalue?

How is the maximum eigenvalue of a random Hermitian matrix distributed? What is the probability that it takes values substantially different from its mean? What about the minimum eigenvalue?

What is the expected spectral norm of a random matrix? What is the probability that the norm takes a value substantially different from its mean?

What about the other eigenvalues or singular values? Can we say something about the “typical” spectrum of a random matrix?

Can we say anything about the eigenvectors or singular vectors? For instance, is each one distributed almost uniformly on the sphere?

We can also ask questions about the operator norm of a random matrix acting as a map between two normed linear spaces. In this case, the geometry of the domain and codomain play a role.

In this work, we focus on the first three questions above. We study the expectation of the extreme eigenvalues of a random Hermitian matrix, and we attempt to provide bounds on the probability that they take an unusual value. As an application of these results, we can control the expected spectral norm of a general matrix and bound the probability of a large deviation. These are the most relevant problems in many (but not all!) applications. The remaining questions are also important, but we will not touch on them here. We recommend the book [Tao12] for a friendly introduction to other branches of random matrix theory.

Random Matrices as Independent Sums

Our approach to random matrices depends on a fundamental principle:

In applications, it is common that a random matrix can be expressed as a sum of independent random matrices.

The examples that appear in these notes should provide ample evidence for this claim. For now, let us describe a specific problem that will serve as an illustration throughout the Introduction. We hope this example is complicated enough to be interesting but simple enough to elucidate the main points.

The star ∗ refers to the conjugate transpose operation, and the standard basis matrix $\mathbf{E}_{jk}$ has a one in the $(j,k)$ position and zeros elsewhere. In other words, the $(j,k)$ entry of the sample covariance matrix $\bm{A}$ records the covariance between the $j$th and $k$th entry of the vector $\bm{x}$.

One basic problem in statistical practice is to estimate the covariance matrix from data. Imagine that we have access to $n$ independent samples $\bm{x}_{1},\dots,\bm{x}_{n}$, each distributed the same way as $\bm{x}$. The sample covariance estimator $\bm{Y}$ is the random matrix

The sample covariance estimator can be expressed as a sum of independent random matrices.

This is precisely the type of decomposition that allows us to apply the tools in these notes.

Exponential Concentration Inequalities for Matrices

An important challenge in probability theory is to study the probability that a real random variable $Z$ takes a value substantially different from its mean. That is, we seek a bound of the form

for a positive parameter $t$. When $Z$ is expressed as a sum of independent random variables, the literature contains many tools for addressing this problem. See [BLM13] for an overview.

For a random matrix $\bm{Z}$, a variant of (1) is the question of whether $\bm{Z}$ deviates substantially from its mean value. We might frame this question as

Here and elsewhere, $\left\|{\cdot}\right\|$ denotes the spectral norm of a matrix. As noted, it is frequently possible to decompose $\bm{Z}$ as a sum of independent random matrices. We might even dream that established methods for studying the scalar concentration problem (1) extend to (2).

To control $Z$, we rely on two types of information: global properties of the sum (such as its mean and variance) and local properties of the summands (such as their maximum fluctuation). These pieces of data are usually easy to obtain. Together, they determine how well $Z$ concentrates around zero, its mean value.

Let $S_{1},\dots,S_{n}$ be independent, centered, real random variables, and assume that each one is uniformly bounded:

Introduce the sum $Z=\sum_{k=1}^{n}S_{k}$, and let $v(Z)$ denote the variance of the sum:

See [BLM13, §2.8] for a proof of this result. We refer to Theorem 6.1 as an exponential concentration inequality because it yields exponentially decaying bounds on the probability that $Z$ deviates substantially from its mean.

The Matrix Bernstein Inequality

What is truly astonishing is that the scalar Bernstein inequality, Theorem 6.1, lifts directly to matrices. Let us emphasize this remarkable fact:

There are exponential concentration inequalities for the spectral norm of a sum of independent random matrices.

As a consequence, once we decompose a random matrix as an independent sum, we can harness global properties (such as the mean and the variance) and local properties (such as a uniform bound on the summands) to obtain detailed information about the norm of the sum. As in the scalar case, it is usually easy to acquire the input data for the inequality. But the output of the inequality is highly nontrivial.

To illustrate these claims, we will state one of the major results from this monograph. This theorem is a matrix extension of Bernstein’s inequality that was developed independently in the two papers [Oli10a, Tro11c]. After presenting the result, we give some more details about its interpretation. In the next section, we apply this result to study the covariance estimation problem.

Let $\bm{S}_{1},\dots,\bm{S}_{n}$ be independent, centered random matrices with common dimension $d_{1}\times d_{2}$, and assume that each one is uniformly bounded

and let $v(\bm{Z})$ denote the matrix variance statistic of the sum:

The proof of this result appears in Chapter 6.

To appreciate what Theorem 6.2 means, it is valuable to make a direct comparison with the scalar version, Theorem 6.1. In both cases, we express the object of interest as an independent sum, and we instate a uniform bound on the summands. There are three salient changes:

The variance $v(\bm{Z})$ in the result for matrices can be interpreted as the magnitude of the expected squared deviation of $\bm{Z}$ from its mean. The formula reflects the fact that a general matrix $\bm{B}$ has two different squares $\bm{BB}^{*}$ and $\bm{B}^{*}\bm{B}$. For an Hermitian matrix, the two squares coincide.

The tail bound has a dimensional factor $d_{1}+d_{2}$ that depends on the size of the matrix. This factor reduces to two in the scalar setting. In the matrix case, it limits the range of $t$ where the tail bound is informative.

The expectation bound (6) is the most important aspect of the matrix Bernstein inequality.

For further discussion of this result, turn to Chapter 6. Chapters 4 and 7 contain related results and interpretations.

Example: The Sample Covariance Estimator

Our goal is to study how the spectral-norm distance $\left\|{\bm{Y}-\bm{A}}\right\|$ between the sample covariance and the true covariance depends on the number $n$ of samples.

We are in a situation where it is quite easy to see how the matrix Bernstein inequality applies. Define the random deviation $\bm{Z}$ of the estimator $\bm{Y}$ from the true covariance matrix $\bm{A}$:

The random matrices $\bm{S}_{k}$ are independent, identically distributed, and centered. To apply Theorem 6.2, we need to find a uniform bound $L$ for the summands, and we need to control the matrix variance statistic $v(\bm{Z})$.

First, let us develop a uniform bound on the spectral norm of each summand. We may calculate that

The first relation is the triangle inequality. The second follows from the assumption that $\bm{x}$ is bounded and the observation that

This expression depends on Jensen’s inequality and the hypothesis that $\bm{x}$ is bounded.

Second, we need to bound the matrix variance statistic $v(\bm{Z})$ defined in (4). The matrix $\bm{Z}$ is Hermitian, so the two squares in this formula coincide with each other:

We need to determine the variance of each summand. By direct calculation,

The expression $\bm{H}\preccurlyeq\bm{T}$ means that $\bm{T}-\bm{H}$ is positive semidefinite. We used the norm bound for the random vector $\bm{x}$ and the fact that expectation preserves the semidefinite order. In the last step, we dropped the negative-semidefinite term $-\bm{A}^{2}$. Summing this relation over $k$, we reach

The matrix is positive-semidefinite because it is a sum of squares of Hermitian matrices. Extract the spectral norm to arrive at

We have now collected the information we need to analyze the sample covariance estimator.

We can invoke the estimate (6) from the matrix Bernstein inequality, Theorem 6.2, with the uniform bound $L=2B/n$ and the variance bound $v(\bm{Z})\leq B\left\|{\bm{A}}\right\|/n$. We attain

In other words, the error in approximating the sample covariance matrix is not too large when we have a sufficient number of samples. If we wish to obtain a relative error on the order of $\varepsilon$, we may take

The analysis in this section applies to many other examples. We encapsulate the argument in Corollary 2.1, which we use to study several more problems.

History of this Example

Covariance estimation may be the earliest application of matrix concentration tools in random matrix theory. Rudelson [Rud99], building on a suggestion of Pisier, showed how to use the noncommutative Khintchine inequality [LP86, LPP91, Buc01, Buc05] to obtain essentially optimal bounds on the sample covariance estimator of a bounded random vector. The tutorial [Ver12] of Roman Vershynin offers an overview of this problem as well as many results and references. The analysis of the sample covariance matrix here is adapted from the technical [GT14]. It leads to a result similar with the one Rudelson obtained in [Rud99].

Optimality of the Matrix Bernstein Inequality

Theorem 6.2 can be sharpened very little because it applies to every random matrix $\bm{Z}$ of the form (3). Let us say a few words about optimality now, postponing the details to §2.

Suppose that $\bm{Z}$ is a random matrix of the form (3). To make the comparison simpler, we also insist that each summand $\bm{S}_{k}$ is a symmetric random variable; that is, $\bm{S}_{k}$ and $-\bm{S}_{k}$ have the same distribution for each index $k$. Introduce the quantity

In §2, we will argue that these assumptions imply

The Arsenal of Results

The Bernstein inequality is probably the most familiar exponential tail bound for a sum of independent random variables, but there are many more. It turns out that essentially all of these scalar results admit extensions that hold for random matrices. In fact, many of the established techniques for scalar concentration have analogs in the matrix setting.

This monograph focuses on a few key exponential concentration inequalities for a sum of independent random matrices, and it describes some specific applications of these results.

A matrix Gaussian series is a random matrix that can be expressed as a sum of fixed matrices, each weighted by an independent standard normal random variable. This formulation includes a surprising number of examples. The most important are undoubtedly Wigner matrices and rectangular Gaussian matrices. Another interesting case is a Toeplitz matrix with Gaussian entries. The analysis of matrix Gaussian series appears in Chapter 4.

A matrix Rademacher series is a random matrix that can be written as a sum of fixed matrices, each weighted by an independent Rademacher random variable. A Rademacher random variable takes the two values $\pm 1$ with equal probability. This construction includes things like random sign matrices, as well as a fixed matrix whose entries are modulated by random signs. There are also interesting examples that arise in combinatorial optimization. We treat these problems in Chapter 4.

The matrix Chernoff bounds apply to a random matrix that can be decomposed as a sum of independent, random positive-semidefinite matrices whose maximum eigenvalues are subject to a uniform bound. These results allow us to obtain information about the norm of a random submatrix drawn from a fixed matrix. They are also appropriate for studying the Laplacian matrix of a random graph. See Chapter 5.

The matrix Bernstein inequality concerns a random matrix that can be expressed as a sum of independent, centered random matrices that admit a uniform spectral-norm bound. This result has many applications, including the analysis of randomized algorithms for matrix sparsification and matrix multiplication. It can also be used to study the random features paradigm for approximating a kernel matrix. Chapter 6 contains this material.

Some matrix concentration inequalities can be improved when the random matrix has limited spectral content in most dimensions. In this situation, we may be able to obtain bounds that do not depend on the ambient dimension. See Chapter 7 for details.

We have chosen to present these results because they are illustrative, and they have already found concrete applications.

What’s Not Here…

The program of extending scalar concentration results to the matrix setting has been quite fruitful, and there are many useful results beyond the ones that we detail. Let us mention some of the other tools that are available. For further information, see the annotated bibliography.

First, there are additional exponential concentration inequalities for a sum of independent random matrices. All of the following results can be established within the framework of this monograph.

Matrix Hoeffding. This result concerns a sum of independent random matrices whose squares are subject to semidefinite upper bounds [Tro11c, §7].

Matrix Bennett. This estimate sharpens the tail bound from the matrix Bernstein inequality [Tro11c, §6].

Matrix Bernstein, Unbounded Case. The matrix Bernstein inequality extends to the case where the moments of the summands grow at a controlled rate. See [Tro11c, §6] or [Kol11].

Matrix Bernstein, Nonnegative Summands. The lower tail of the Bernstein inequality can be improved when the summands are positive semidefinite [Mau03]; this result extends to the matrix setting. By a different argument, the dimensional factor can be removed from this bound for a class of interesting examples [Oli13, Thm. 3.1].

The approach in this monograph can be adapted to obtain exponential concentration for matrix-valued martingales. Here are a few results from this category:

Matrix Azuma. This is the martingale version of the matrix Hoeffding bound [Tro11c, §7].

Matrix Bounded Differences. The matrix Azuma inequality gives bounds for the spectral norm of a matrix-valued function of independent random variables [Tro11c, §7].

Matrix Freedman. This result can be viewed as the martingale extension of the matrix Bernstein inequality [Oli10a, Tro11a].

The technical report [Tro11b] explains how to extend other bounds for a sum of independent random matrices to the martingale setting.

Polynomial moment inequalities provide bounds for the expected trace of a power of a random matrix. Moment inequalities for a sum of independent random matrices can provide useful information when the summands have heavy tails or else a uniform bound does not reflect the typical size of the summands.

Matrix Khintchine. The matrix Khintchine inequality is the polynomial version of the exponential bounds for matrix Gaussian series and matrix Rademacher series. This result is presented in (1). See the papers [LP86, Buc01, Buc05] or [MJC+14, Cor. 7.3] for proofs.

Matrix Moment Inequalities. The matrix Chernoff inequality admits a polynomial variant; the simplest form appears in (9). The matrix Bernstein inequality also has a polynomial variant, stated in (6). These bounds are drawn from [CGT12a, App.].

The methods that lead to polynomial moment inequalities differ substantially from the techniques in this monograph, so we cannot include the proofs here. The annotated bibliography includes references to the large literature on moment inequalities for random matrices.

Recently, Lester Mackey and the author, in collaboration with Daniel Paulin and several other researchers [MJC+14, PMT14], have developed another framework for establishing matrix concentration. This approach extends a scalar argument, introduced by Chatterjee [Cha05, Cha07], that depends on exchangeable pairs and Markov chain couplings. The method of exchangeable pairs delivers both exponential concentration inequalities and polynomial moment inequalities for random matrices, and it can reproduce many of the bounds mentioned above. It also leads to new results:

Polynomial Efron–Stein Inequality for Matrices. This bound is a matrix version of the polynomial Efron–Stein inequality [BBLM05, Thm. 1]. It controls the polynomial moments of a centered random matrix that is a function of independent random variables [PMT14, Thm. 4.2].

Exponential Efron–Stein Inequality for Matrices. This bound is the matrix extension of the exponential Efron–Stein inequality [BLM03, Thm. 1]. It leads to exponential concentration inequalities for a centered random matrix constructed from independent random variables [PMT14, Thm. 4.3].

Another significant advantage is that the method of exchangeable pairs can sometimes handle random matrices built from dependent random variables. Although the simplest version of the exchangeable pairs argument is more elementary than the approach in this monograph, it takes a lot of effort to establish the more useful inequalities. With some regret, we have chosen not to include this material because the method and results are accessible to a narrower audience.

Finally, we remark that the modified logarithmic Sobolev inequalities of [BLM03, BBLM05] also extend to the matrix setting [CT14]. Unfortunately, the matrix variants do not seem to be as useful as the scalar results.

About This Monograph

This monograph is intended for graduate students and researchers in computational mathematics who want to learn some modern techniques for analyzing random matrices. The preparation required is minimal. We assume familiarity with calculus, applied linear algebra, the basic theory of normed spaces, and classical probability theory up through the elementary concentration inequalities (such as Markov and Bernstein). Beyond the basics, which can be gleaned from any good textbook, we include all the required background in Chapter 2.

The material here is based primarily on the paper “User-Friendly Tail Bounds for Sums of Random Matrices” by the present author [Tro11c]. There are several significant revisions to this earlier work:

Many of the papers on matrix concentration give limited information about how the results can be used to solve problems of interest. A major part of these notes consists of worked examples and applications that indicate how matrix concentration inequalities apply to practical questions.

This work collects bounds for the expected value of the spectral norm of a random matrix and bounds for the expectation of the smallest and largest eigenvalues of a random symmetric matrix. Some of these useful results have appeared piecemeal in the literature [CGT12a, MJC+14], but they have not been included in a unified presentation.

We explain why each matrix concentration inequality is (nearly) optimal. This presentation includes examples to show that each term in each bound is necessary to describe some particular phenomenon.

Over the last few years, there have been some refinements to the basic matrix concentration bounds that improve the dependence on dimension [HKZ12, Min11]. We describe a new framework that allows us to prove these results with ease.

The matrix concentration inequalities in this monograph depend on a deep theorem [Lie73, Thm. 6] from matrix analysis due to Elliott Lieb. We provide a complete proof of this result, along with all the background required to understand the argument.

We have included a list of the major works on matrix concentration, including a short summary of the main contributions of these papers. We hope this catalog will be a valuable guide for further reading.

The organization of the notes is straightforward. Chapter 2 contains background material that is needed for the analysis. Chapter 3 describes the framework for developing exponential concentration inequalities for matrices. Chapter 4 presents the first set of results and examples, concerning matrix Gaussian and Rademacher series. Chapter 5 introduces the matrix Chernoff bounds and their applications, and Chapter 6 expands on our discussion of the matrix Bernstein inequality. Chapter 7 shows how to sharpen some of the results so that they depend on an intrinsic dimension parameter. Chapter 8 contains the proof of Lieb’s theorem. We conclude with resources on matrix concentration and a bibliography.

To make the presentation smoother, we have not followed all of the conventions for scholarly articles in journals. In particular, almost all the citations appear in the notes at the end of each chapter. Our aim has been to explain the ideas as clearly as possible, rather than to interrupt the narrative with an elaborate genealogy of results.

Chapter 2 Matrix Functions & Probability with Matrices

We begin the main development with a short overview of the background material that is required to understand the proofs and, to a lesser extent, the statements of matrix concentration inequalities. We have been careful to provide cross-references to these foundational results, so most readers will be able to proceed directly to the main theoretical development in Chapter 3 or the discussion of specific random matrix inequalities in Chapters 4, 5, and 6.

Section 1 covers material from matrix theory concerning the behavior of matrix functions. Section 2 reviews relevant results from probability, especially the parts involving matrices.

Matrix Theory Background

Let us begin with the results we require from the field of matrix analysis.

Spaces of Vectors

Spaces of Matrices

Topology & Convergence

We can endow the space of matrices with the Frobenius norm:

Basic Vectors and Matrices

A square matrix $\bm{Q}$ that satisfies $\bm{QQ}^{*}=\mathbf{I}=\bm{Q}^{*}\bm{Q}$ is called a unitary matrix. We reserve the letter $\bm{Q}$ for a unitary matrix. Readers who prefer the real setting may prefer to regard $\bm{Q}$ as an orthogonal matrix.

Hermitian Matrices and Eigenvalues

An Hermitian matrix $\bm{A}$ is a square matrix that satisfies $\bm{A}=\bm{A}^{*}$. A useful intuition from operator theory is that Hermitian matrices are analogous with real numbers, while general square matrices are analogous with complex numbers.

The diagonal entries of $\bm{\Lambda}$ are real numbers, which are referred to as the eigenvalues of $\bm{A}$. The unitary matrix $\bm{Q}$ in the eigenvalue decomposition is not determined completely, but the list of eigenvalues is unique modulo permutations. The eigenvalues of an Hermitian matrix are often referred to as its spectrum.

We denote the algebraic minimum and maximum eigenvalues of an Hermitian matrix $\bm{A}$ by $\lambda_{\min}(\bm{A})$ and $\lambda_{\max}(\bm{A})$. The extreme eigenvalue maps are positive homogeneous:

There is an important relationship between minimum and maximum eigenvalues:

The fact (5) warns us that we must be careful passing scalars through an eigenvalue map.

This work rarely requires any eigenvalues of an Hermitian matrix aside from the minimum and maximum. When they do arise, we usually order the other eigenvalues in the weakly decreasing sense:

On occasion, it is more natural to arrange eigenvalues in the weakly increasing sense:

To prevent confusion, we will accompany this notation with a reminder.

Readers who prefer the real setting may read “symmetric” in place of “Hermitian.” In this case, the eigenvalue decomposition involves an orthogonal matrix $\bm{Q}$. Note, however, that the term “symmetric” has a different meaning when applied to random variables!

The Trace of a Square Matrix

The trace of a square matrix, denoted by $\operatorname{tr}$, is the sum of its diagonal entries.

In particular, the existence of an eigenvalue decomposition (3) shows that the trace of an Hermitian matrix equals the sum of its eigenvalues. This fact also holds true for a general square matrix.

Another valuable relation connects the trace with the Frobenius norm:

This expression follows from the definitions (2) and (6) and a short calculation.

The Semidefinite Partial Order

Equivalently, $\bm{A}$ is positive definite when it is Hermitian and its eigenvalues are all positive.

Positive-semidefinite and positive-definite matrices play a special role in matrix theory, analogous with the role of nonnegative and positive numbers in real analysis. In particular, observe that the square of an Hermitian matrix is always positive semidefinite. The square of a nonsingular Hermitian matrix is always positive definite.

In particular, we write $\bm{A}\succcurlyeq\bm{0}$ to indicate that $\bm{A}$ is positive semidefinite and $\bm{A}\succ\bm{0}$ to indicate that $\bm{A}$ is positive definite. For a diagonal matrix $\bm{\Lambda}$, the expression $\bm{\Lambda}\succcurlyeq\bm{0}$ means that each entry of $\bm{\Lambda}$ is nonnegative.

The semidefinite order is preserved by conjugation, a simple fact whose importance cannot be overstated.

Let $\bm{A}$ and $\bm{H}$ be Hermitian matrices of the same dimension, and let $\bm{B}$ be a general matrix with compatible dimensions. Then

Finally, we remark that the trace of a positive-semidefinite matrix is at least as large as its maximum eigenvalue:

This property follows from the definition of a positive-semidefinite matrix and the fact that the trace of $\bm{A}$ equals the sum of the eigenvalues.

Standard Matrix Functions

Let us describe the most direct method for extending a function on the real numbers to a function on Hermitian matrices. The basic idea is to apply the function to each eigenvalue of the matrix to construct a new matrix.

In particular, we can apply $f$ to a real diagonal matrix by applying the function to each diagonal entry.

It can be verified that the definition of $f(\bm{A})$ does not depend on which eigenvalue decomposition $\bm{A}=\bm{Q\Lambda Q}^{*}$ that we choose. Any matrix function that arises in this fashion is called a standard matrix function.

To confirm that this definition is sensible, consider the power function $f(t)=t^{q}$ for a natural number $q$. When $\bm{A}$ is Hermitian, the power function $f(\bm{A})=\bm{A}^{q}$, where $\bm{A}^{q}$ is the $q$-fold product of $\bm{A}$.

The following result is an immediate, but important, consequence of the definition of a standard matrix function.

This formula can be verified using an eigenvalue decomposition of $\bm{A}$ and the definition of a standard matrix function.

The Transfer Rule

In most cases, the “obvious” generalization of an inequality for real-valued functions fails to hold in the semidefinite order. Nevertheless, there is one class of inequalities for real functions that extends to give semidefinite relationships for standard matrix functions.

Let $f$ and $g$ be real-valued functions defined on an interval $I$ of the real line, and let $\bm{A}$ be an Hermitian matrix whose eigenvalues are contained in $I$. Then

Decompose $\bm{A}=\bm{Q\Lambda Q}^{*}$. It is immediate that $f(\bm{\Lambda})\preccurlyeq g(\bm{\Lambda})$. The Conjugation Rule (12) allows us to conjugate this relation by $\bm{Q}$. Finally, we invoke Definition 1.2, of a standard matrix function, to complete the argument. ∎

The Matrix Exponential

The Spectral Mapping Theorem, Proposition 1.3, implies that the exponential of an Hermitian matrix is always positive definite.

We often work with the trace of the matrix exponential:

This function has a monotonicity property that we use extensively. For Hermitian matrices $\bm{A}$ and $\bm{H}$ with the same dimension,

The Matrix Logarithm

We can define the matrix logarithm as a standard matrix function. The matrix logarithm is also the functional inverse of the matrix exponential:

A valuable fact about the matrix logarithm is that it preserves the semidefinite order. For positive-definite matrices $\bm{A}$ and $\bm{H}$ with the same dimension,

We establish this result in §4. Let us stress that the matrix exponential does not have any operator monotonicity property analogous with (18)!

Singular Values of Rectangular Matrices

A general matrix does not have an eigenvalue decomposition, but it admits a different representation that is just as useful. Every $d_{1}\times d_{2}$ matrix $\bm{B}$ has a singular value decomposition

The unitary matrices $\bm{Q}_{1}$ and $\bm{Q}_{2}$ have dimensions $d_{1}\times d_{1}$ and $d_{2}\times d_{2}$, respectively. The inner matrix $\bm{\Sigma}$ has dimension $d_{1}\times d_{2}$, and we use the term diagonal in the sense that only the diagonal entries $(\bm{\Sigma})_{jj}$ may be nonzero.

The diagonal entries of $\bm{\Sigma}$ are called the singular values of $\bm{B}$, and they are denoted as $\sigma_{j}(\bm{B})$. The singular values are determined completely modulo permutations, and it is conventional to arrange them in weakly decreasing order:

There is an important relationship between singular values and eigenvalues. A general matrix has two squares associated with it, $\bm{BB}^{*}$ and $\bm{B}^{*}\bm{B}$, both of which are positive semidefinite. We can use a singular value decomposition of $\bm{B}$ to construct eigenvalue decompositions of the two squares:

The two squares of $\bm{\Sigma}$ are square, diagonal matrices with nonnegative entries. Conversely, we can always extract a singular value decomposition from the eigenvalue decompositions of the two squares.

We can write the Frobenius norm of a matrix in terms of the singular values:

This expression follows from the expression (8) for the Frobenius norm, the property (20) of the singular value decomposition, and the unitary invariance (7) of the trace.

The Spectral Norm

The spectral norm of an Hermitian matrix $\bm{A}$ is defined by the relation

For a general matrix $\bm{B}$, the spectral norm is defined to be the largest singular value:

The Stable Rank

In several of the applications, we need an analytic measure of the collinearity of the rows and columns of a matrix called the stable rank. For a general matrix $\bm{B}$, the stable rank is defined as

The stable rank is a lower bound for the algebraic rank:

This point follows when we use (21) and (23) to express the two norms in terms of the singular values of $\bm{B}$. In contrast to the algebraic rank, the stable rank is a continuous function of the matrix, so it is more suitable for numerical applications.

Dilations

An extraordinarily fruitful idea from operator theory is to embed matrices within larger block matrices, called dilations. Dilations have an almost magical power. In this work, we will use dilations to extend matrix concentration inequalities from Hermitian matrices to general matrices.

is the map from a general matrix to an Hermitian matrix defined by

It is clear that the Hermitian dilation is a real-linear map. Furthermore, the dilation retains important spectral information. To see why, note that the square of the dilation satisfies

We discover that the squared eigenvalues of $\mathscr{H}(\bm{B})$ coincide with the squared singular values of $\bm{B}$, along with an appropriate number of zeros. As a consequence, $\left\|{\mathscr{H}(\bm{B})}\right\|=\left\|{\bm{B}}\right\|$. Moreover,

We will invoke the identity (28) repeatedly.

One way to justify the first relation in (28) is to introduce the first columns $\bm{u}_{1}$ and $\bm{u}_{2}$ of the unitary matrices $\bm{Q}_{1}$ and $\bm{Q}_{2}$ that appear in the singular value decomposition $\bm{B}=\bm{Q}_{1}\bm{\Sigma}\bm{Q}_{2}^{*}$. Then we may calculate that

Indeed, the spectral norm of $\bm{B}$ equals its largest singular value $\sigma_{1}(\bm{B})$, which coincides with $\bm{u}_{1}^{*}\bm{B}\bm{u}_{2}$ by construction of $\bm{u}_{1}$ and $\bm{u}_{2}$. The second identity relies on a direct calculation. The first inequality follows from the variational representation of the maximum eigenvalue as a Rayleigh quotient; this fact can also be derived as a consequence of (3). The second inequality depends on the definition (22) of the spectral norm of an Hermitian matrix.

Other Matrix Norms

There are a number of other matrix norms that arise sporadically in this work. The Schatten 1-norm of a matrix can be defined as the sum of its singular values:

because of the Cauchy–Schwarz inequality.

Probability with Matrices

We continue with some material from probability, focusing on connections with matrices.

We prefer to avoid abstraction and unnecessary technical detail, so we frame the standing assumption that all random variables are sufficiently regular that we are justified in computing expectations, interchanging limits, and so forth. The manipulations we perform are valid if we assume that all random variables are bounded, but the results hold in broader circumstances if we instate appropriate regularity conditions.

Since the expectation operator is linear, we typically do not use parentheses with it. We instate the convention that powers and products take precedence over the expectation operator. In particular,

This position helps us reduce the clutter of parentheses. We sometimes include extra delimiters when it is helpful for clarity.

Some Scalar Random Variables

We use consistent notation for some of the basic scalar random variables.

We reserve the letter $\gamma$ for a $\textsc{normal}(0,1)$ random variable. That is, $\gamma$ is a real Gaussian with mean zero and variance one.

We reserve the letter $\varrho$ for a random variable that takes the two values $\pm 1$ with equal probability.

A $\textsc{bernoulli}(p)$ random variable takes the value one with probability $p$ and the value zero with probability $1-p$, where $p\in$. We use the letters $\delta$ and $\xi$ for Bernoulli random variables.

Random Matrices

It is more natural to think of the entries of $\bm{Z}$ as complex random variables that may or may not be correlated with each other. We reserve the letters $\bm{X}$ and $\bm{Y}$ for random Hermitian matrices, while the letter $\bm{Z}$ denotes a general random matrix.

A finite sequence $\{\bm{Z}_{k}\}$ of random matrices is independent when

Expectation

The expectation of a random matrix $\bm{Z}=[Z_{jk}]$ is simply the matrix formed by taking the componentwise expectation. That is,

Under mild assumptions, expectation commutes with linear and real-linear maps. Indeed, expectation commutes with multiplication by a fixed matrix:

In particular, the product rule for the expectation of independent random variables extends to matrices:

We use these identities liberally, without any further comment.

Inequalities for Expectation

Markov’s inequality states that a nonnegative (real) random variable $X$ obeys the probability bound

The Markov inequality is a central tool for establishing concentration inequalities.

Jensen’s inequality describes how averaging interacts with convexity. Let $\bm{Z}$ be a random matrix, and let $h$ be a real-valued function on matrices. Then

We use this result many times without direct reference.

The Variance of a Random Hermitian Matrix

The variance of a real random variable $Y$ is defined as the expected squared deviation from the mean:

There are a number of natural extensions of this concept in the matrix setting that play a role in our theory.

Suppose that $\bm{Y}$ is a random Hermitian matrix. We can define a matrix-valued variance:

The matrix $\operatorname{\mathbf{Var}}(\bm{Y})$ is always positive semidefinite. We can interpret the $(j,k)$ entry of this matrix as the covariance between the $j$th and $k$th columns of $\bm{Y}$:

where we have written $\bm{y}_{:j}$ for the $j$th column of $\bm{Y}$.

The matrix-valued variance contains a lot of information about the fluctuations of the random matrix. We can summarize $\operatorname{\mathbf{Var}}(\bm{Y})$ using a single number $v(\bm{Y})$, which we call the matrix variance statistic:

To understand what this quantity means, one may wish to rewrite it as

Roughly speaking, the matrix variance statistic describes the maximum variance of $\bm{Y}\bm{u}$ for any unit vector $\bm{u}$.

The Variance of a Sum of Independent, Random Hermitian Matrices

The matrix-valued variance interacts beautifully with a sum of independent random matrices. Consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random Hermitian matrices with common dimension $d$. Introduce the sum $\bm{Y}=\sum_{k}\bm{X}_{k}$. Then

This identity matches the familiar result for the variance of a sum of independent scalar random variables. It follows that the matrix variance statistic satisfies

The fact that the sum remains inside the norm is very important. Indeed, the best general inequalities between $v(\bm{Y})$ and the matrix variance statistics $v(\bm{X}_{k})$ of the summands are

These relations can be improved in some special cases. For example, when the matrices $\bm{X}_{k}$ are identically distributed, the left-hand inequality becomes an identity.

The Variance of a Rectangular Random Matrix

We will often work with non-Hermitian random matrices. In this case, we need to account for the fact that a general matrix has two different squares. Suppose that $\bm{Z}$ is a random matrix with dimension $d_{1}\times d_{2}$. Define

The matrix $\operatorname{\mathbf{Var}}_{1}(\bm{Z})$ is a positive-semidefinite matrix with dimension $d_{1}\times d_{1}$, and it describes the fluctuation of the rows of $\bm{Z}$. The matrix $\operatorname{\mathbf{Var}}_{2}(\bm{Z})$ is a positive-semidefinite matrix with dimension $d_{2}\times d_{2}$, and it reflects the fluctuation of the columns of $\bm{Z}$. For an Hermitian random matrix $\bm{Y}$,

In other words, the two variances coincide in the Hermitian setting.

As before, it is valuable to reduce these matrix-valued variances to a single scalar parameter. We define the matrix variance statistic of a general random matrix $\bm{Z}$ as

When $\bm{Z}$ is Hermitian, the definition (8) coincides with the original definition (4).

To promote a deeper appreciation for the formula (8), let us explain how it arises from the Hermitian dilation (26). By direct calculation,

The first identity is the definition (3) of the matrix-valued variance. The second line follows from the formula (27) for the square of the dilation. The last identity depends on the definition (7) of the two matrix-valued variances. Therefore, using the definitions (4) and (8) of the matrix variance statistics,

The second identity holds because the spectral norm of a block-diagonal matrix is the maximum norm achieved by one of the diagonal blocks.

The Variance of a Sum of Independent Random Matrices

As in the Hermitian case, the matrix-valued variances interact nicely with an independent sum. Consider a finite sequence $\{\bm{S}_{k}\}$ of independent random matrices with the same dimension. Form the sum $\bm{Z}=\sum_{k}\bm{S}_{k}$. Repeating the calculation leading up to (6), we find that

In summary, the matrix variance statistic of an independent sum satisfies

Notes

Everything in this chapter is firmly established. We have culled the results that are relevant to our discussion. Let us give some additional references for readers who would like more information.

Our treatment of matrix analysis is drawn from Bhatia’s excellent books on matrix analysis [Bha97, Bha07]. The two books [HJ13, HJ94] of Horn & Johnson also serve as good general references. Higham’s work [Hig08] is a generous source of information about matrix functions. Other valuable resources include Carlen’s lecture notes [Car10], the book of Petz [Pet11], and the book of Hiai & Petz [HP14].

Probability with Matrices

The classic introduction to probability is the two-volume treatise [Fel68, Fel71] of Feller. The book [GS01] of Grimmett & Stirzaker offers a good treatment of probability theory and random processes at an intermediate level. For a more theoretical presentation, consider the book [Shi96] of Shiryaev.

There are too many books on random matrix theory for us to include a comprehensive list; here is a selection that the author finds useful. Tao’s book [Tao12] gives a friendly introduction to some of the major aspects of classical and modern random matrix theory. The lecture notes [Kem13] of Kemp are also extremely readable. The survey of Vershynin [Ver12] provides a good summary of techniques from asymptotic convex geometry that are relevant to random matrix theory. The works of Mardia, Kent, & Bibby [MKB79] and Muirhead [Mui82] present classical results on random matrices that are particularly useful in statistics, while Bai & Silverstein [BS10] contains a comprehensive modern treatment. Nica and Speicher [NS06] offer an entrée to the beautiful field of free probability. Mehta’s treatise [Meh04] was the first book on random matrix theory available, and it remains solid.

Chapter 3 The Matrix Laplace Transform Method

This chapter contains the core part of the analysis that ultimately delivers matrix concentration inequalities. Readers who are only interested in the concentration inequalities themselves or the example applications may wish to move on to Chapters 4, 5, and 6.

In the scalar setting, the Laplace transform method provides a simple but powerful way to develop concentration inequalities for a sum of independent random variables. This technique is sometimes referred to as the “Bernstein trick” or “Chernoff bounding.” For a primer, we recommend [BLM13, Chap. 2].

In the matrix setting, there is a very satisfactory extension of this argument that allows us to prove concentration inequalities for a sum of independent random matrices. As in the scalar case, the matrix Laplace transform method is both easy to use and incredibly useful. In contrast to the scalar case, the arguments that lead to matrix concentration are no longer elementary. The purpose of this chapter is to install the framework we need to support these results. Fortunately, in practical applications, all of the technical difficulty remains invisible.

We first define matrix analogs of the moment generating function and the cumulant generating function, which pack up information about the fluctuations of a random Hermitian matrix. Section 2 explains how we can use the matrix mgf to obtain probability inequalities for the maximum eigenvalue of a random Hermitian matrix. The next task is to develop a bound for the mgf of a sum of independent random matrices using information about the summands. In §3, we discuss the challenges that arise; §4 presents the ideas we need to overcome these obstacles. Section 5 establishes that the classical result on additivity of cumulants has a companion in the matrix setting. This result allows us to develop a collection of abstract probability inequalities in §6 that we can specialize to obtain matrix Chernoff bounds, matrix Bernstein bounds, and so forth.

Matrix Moments and Cumulants

At the heart of the Laplace transform method are the moment generating function (mgf) and the cumulant generating function (cgf) of a random variable. We begin by presenting matrix versions of the mgf and cgf.

Let $\bm{X}$ be a random Hermitian matrix. The matrix moment generating function $\bm{M}_{\bm{X}}$ and the matrix cumulant generating function $\bm{\Xi}_{\bm{X}}$ are given by

Note that the expectations may not exist for all values of $\theta$.

The matrix mgf $\bm{M}_{\bm{X}}$ and matrix cgf $\bm{\Xi}_{\bm{X}}$ contain information about how much the random matrix $\bm{X}$ varies. We aim to exploit the data encoded in these functions to control the eigenvalues.

Let us take a moment to expand on Definition 1.1; this discussion is not important for subsequent developments. Observe that the matrix mgf and cgf have formal power series expansions:

The matrix variance was introduced in (3). Higher-order cumulants are harder to write down and interpret.

The Matrix Laplace Transform Method

In the scalar setting, the Laplace transform method allows us to obtain tail bounds for a random variable in terms of its mgf. The starting point for our theory is the observation that a similar result holds in the matrix setting.

In words, we can control the tail probabilities of the extreme eigenvalues of a random matrix by producing a bound for the trace of the matrix mgf. The proof of this fact parallels the classical argument, but there is a twist.

We begin with (1). Fix a positive number $\theta$, and observe that

The first identity depends on the Spectral Mapping Theorem, Proposition 1.3, and the fact that the exponential function is increasing. The inequality follows because the exponential of an Hermitian matrix is positive definite, and (13) shows that the maximum eigenvalue of a positive-definite matrix is dominated by the trace. Combine the latter two displays to reach

This inequality is valid for any positive $\theta$, so we may take an infimum to achieve the tightest possible bound.

To prove (2), we use a similar approach. Fix a negative number $\theta$, and calculate that

In the proof of Proposition 2.1, it may seem crude to bound the maximum eigenvalue by the trace. In fact, our overall approach leads to matrix concentration inequalities that are sharp for specific examples (see the discussion in §§2, 2, and 2), so we must conclude that the loss in this bound is sometimes inevitable. At the same time, this maneuver allows us to exploit some amazing convexity properties of the trace exponential.

We can adapt the proof of Proposition 2.1 to obtain bounds for the expectation of the maximum eigenvalue of a random Hermitian matrix. This argument does not have a perfect analog in the scalar setting.

Let $\bm{Y}$ be a random Hermitian matrix. Then

We establish the bound (4); the proof of (5) is quite similar. Fix a positive number $\theta$, and calculate that

The first identity holds because the maximum eigenvalue is a positive-homogeneous map, as stated in (4). The second relation is Jensen’s inequality. The third follows when we use the Spectral Mapping Theorem, Proposition 1.3, to draw the eigenvalue map through the exponential. The final inequality depends on the fact (13) that the trace of a positive-definite matrix dominates the maximum eigenvalue. ∎

The Failure of the Matrix Mgf

We would like the use the Laplace transform bounds from Section 2 to study a sum of independent random matrices. In the scalar setting, the Laplace transform method is effective for studying an independent sum because the mgf and the cgf decompose. In the matrix case, the situation is more subtle, and the goal of this section is to indicate where things go awry.

Consider an independent sequence $\{X_{k}\}$ of real random variables. The mgf of the sum satisfies a multiplication rule:

The first identity is the definition of an mgf. The second relation holds because the exponential map converts a sum of real scalars to a product, and the third relation requires the independence of the random variables. The last identity, again, is the definition.

At first, we might imagine that a similar relationship holds for the matrix mgf. Consider an independent sequence $\{\bm{X}_{k}\}$ of random Hermitian matrices. Perhaps,

Unfortunately, this hope shatters when we subject it to interrogation.

It is not hard to find the reason that (2) fails. The identity (1) depends on the fact that the scalar exponential converts a sum into a product. In contrast, for Hermitian matrices,

If we introduce the trace, the situation improves somewhat:

The result (3) is known as the Golden–Thompson inequality, a famous theorem from statistical physics. Unfortunately, the analogous bound may fail for three matrices:

It seems that we have reached an impasse.

What if we consider the cgf instead? The cgf of a sum of independent real random variables satisfies an addition rule:

The relation (4) follows when we extract the logarithm of the multiplication rule (1). This result looks like a more promising candidate for generalization because a sum of Hermitian matrices remains Hermitian. We might hope that

As stated, this putative identity also fails. Nevertheless, the addition rule (4) admits a very satisfactory extension to matrices. In contrast with the scalar case, the proof involves much deeper considerations.

A Theorem of Lieb

To find the appropriate generalization of the addition rule for cgfs, we turn to the literature on matrix analysis. Here, we discover a famous result of Elliott Lieb on the convexity properties of the trace exponential function.

Fix an Hermitian matrix $\bm{H}$ with dimension $d$. The function

is a concave map on the convex cone of $d\times d$ positive-definite matrices.

In the scalar case, the analogous function $a\mapsto\exp(h+\log a)$ is linear, so this result describes a new type of phenomenon that emerges when we move to the matrix setting. We present a complete proof of Theorem 4.1 in Chapter 8.

For now, let us focus on the consequences of this remarkable result. Lieb’s Theorem is valuable to us because the Laplace transform bounds from Section 2 involve the trace exponential function. To highlight the connection, we rephrase Theorem 4.1 in probabilistic terms.

Let $\bm{H}$ be a fixed Hermitian matrix, and let $\bm{X}$ be a random Hermitian matrix of the same dimension. Then

The first identity follows from the interpretation (17) of the matrix logarithm as the functional inverse of the matrix exponential. Theorem 4.1 shows that the trace function is concave in $\bm{Y}$, so Jensen’s inequality (2) allows us to draw the expectation inside the function. ∎

Subadditivity of the Matrix Cgf

We are now prepared to generalize the addition rule (4) for scalar cgfs to the matrix setting. The following result is fundamental to our approach to random matrices.

Consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random, Hermitian matrices of the same dimension. Then

The parallel between the additivity rule (4) and the subadditivity rule (2) is striking. With our level of preparation, it is easy to prove this result. We just apply the bound from Corollary 4.2 repeatedly.

This argument is legitimate because $\bm{H}_{m}$ is independent from $\bm{X}_{m}$.

The formulation (2) follows from (1) when we substitute the expression (1) for the matrix cgf and make some algebraic simplifications. ∎

Master Bounds for Sums of Independent Random Matrices

Finally, we can present some general results on the behavior of a sum of independent random matrices. At this stage, we simply combine the Laplace transform bounds with the subadditivity of the matrix cgf to obtain abstract inequalities. Later, we will harness properties of the summands to develop more concrete estimates that apply to specific examples of interest.

Consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random, Hermitian matrices of the same size. Then

Substitute the subadditivity rule for matrix cgfs, Lemma 5.1, into the two matrix Laplace transform results, Proposition 2.1 and Proposition 2.2. ∎

In this chapter, we have focused on probability inequalities for the extreme eigenvalues of a sum of independent random matrices. Nevertheless, these results also give information about the spectral norm of a sum of independent, random, rectangular matrices because we can apply them to the Hermitian dilation (26) of the sum. Instead of presenting a general theorem, we find it more natural to extend individual results to the non-Hermitian case.

Notes

This section includes some historical discussion about the results we have described in this chapter, along with citations for the results that we have established.

The idea of lifting the “Bernstein trick” to the matrix setting is due to two researchers in quantum information theory, Rudolf Ahlswede and Andreas Winter, who were working on a problem concerning transmission of information through a quantum channel [AW02]. Their paper contains a version of the matrix Laplace transform result, Proposition 2.1, along with a substantial number of related foundational ideas. Their work is one of the major inspirations for the tools that are described in these notes.

The statement of Proposition 2.1 and the proof that we present appear in the paper [Oli10b] of Roberto Oliveira. The subsequent result on expectations, Proposition 2.2, first appeared in the paper [CGT12a].

Subadditivity of Cumulants

The major impediment to applying the matrix Laplace transform method is the need to produce a bound for the trace of the matrix moment generating function (the trace mgf). This is where all the technical difficulty in the argument resides.

Ahlswede & Winter [AW02, App.] proposed an approach for bounding the trace mgf of an independent sum, based on a repeated application of the Golden–Thompson inequality (3). Their argument leads to a cumulant bound of the form

when the random Hermitian matrices $\bm{X}_{k}$ have dimension $d$. In other words, Ahlswede & Winter bound the cumulant of a sum in terms of the sum of the maximum eigenvalues of the cumulants. There are cases where the bound (1) is equivalent with Lemma 5.1. For example, the estimates coincide when each matrix $\bm{X}_{k}$ is identically distributed. In general, however, the estimate (1) leads to fundamentally weaker results than our bound from Lemma 5.1. In the worst case, the approach of Ahlswede & Winter may produce an unnecessary factor of the dimension $d$ in the exponent. See [Tro11c, §§3.7, 4.8] for details.

The first major technical advance beyond the original argument of Ahlswede & Winter appeared in a paper [Oli10a] of Oliveira. He developed a more effective way to deploy the Golden–Thompson inequality, and he used this technique to establish a matrix version of Freedman’s inequality [Fre75]. In the scalar setting, Freedman’s inequality extends the Bernstein concentration inequality to martingales; Oliveira obtained the analogous extension of Bernstein’s inequality for matrix-valued martingales. When specialized to independent sums, his result is quite similar with the matrix Bernstein inequality, Theorem 6.2, apart from the precise values of the constants. Oliveira’s method, however, does not seem to deliver the full spectrum of matrix concentration inequalities that we discuss in these notes.

The approach here, based on Lieb’s Theorem, was introduced in the article [Tro11c] by the author of these notes. This paper was apparently the first to recognize that Lieb’s Theorem has probabilistic content, as stated in Corollary 4.2. This idea leads to Lemma 5.1, on the subadditivity of cumulants, along with the master tail bounds from Theorem 6.1. Note that the two articles [Oli10a, Tro11c] are independent works.

For a detailed discussion of the benefits of Lieb’s Theorem over the Golden–Thompson inequality, see [Tro11c, §4]. In summary, to get the sharpest concentration results for random matrices, Lieb’s Theorem appears to be indispensible. The approach of Ahlswede & Winter seems intrinsically weaker. Oliveira’s argument has certain advantages, however, in that it extends from matrices to the fully noncommutative setting [JZ12].

Subsequent research on the underpinnings of the matrix Laplace transform method has led to a martingale version of the subadditivity of cumulants [Tro11a, Tro11b]; these works also depend on Lieb’s Theorem. The technical report [GT14] shows how to use a related result, called the Lieb–Seiringer Theorem [LS05], to obtain upper and lower tail bounds for all eigenvalues of a sum of independent random Hermitian matrices.

Noncommutative Moment Inequalities

There is a closely related, and much older, line of research on noncommutative moment inequalities. These results provide information about the expected trace of a power of a sum of independent random matrices. The matrix Laplace transform method, as encapsulated in Theorem 6.1, gives analogous bounds for the exponential moments.

Research on noncommutative moment inequalities dates to an important paper [LP86] of Françoise Lust-Piquard, which contains an operator extension of the Khintchine inequality. Her result, now called the noncommutative Khintchine inequality, controls the trace moments of a sum of fixed matrices, each modulated by an independent Rademacher random variable; see Section 2 for more details.

In recent years, researchers have generalized many other moment inequalities for a sum of scalar random variables to matrices (and beyond). For instance, the Rosenthal–Pinelis inequality for a sum of independent zero-mean random variables admits a matrix version [JZ13, MJC+14, CGT12a]. We present a variant of the latter result below in (6). See the paper [JX05] for a good overview of some other noncommutative moment inequalities.

Finally, and tangentially, we mention that a different notion of matrix moments and cumulants plays a central role in the theory of free probability [NS06].

Quantum Statistical Mechanics

A curious feature of the theory of matrix concentration inequalities is that the most powerful tools come from the mathematical theory of quantum statistical mechanics. This field studies the bulk statistical properties of interacting quantum systems, and it would seem quite distant from the field of random matrix theory. The connection between these two areas has emerged because of research on quantum information theory, which studies how information can be encoded, operated upon, and transmitted via quantum mechanical systems.

The Golden–Thompson inequality is a major result from quantum statistical mechanics. Bhatia’s book [Bha97, Sec. IX.3] contains a detailed treatment of this result from the perspective of matrix theory. For an account with more physical content, see the book of Thirring [Thi02]. The fact that the Golden–Thompson inequality fails for three matrices can be obtained from simple examples, such as combinations of Pauli spin matrices [Bha97, Exer. IX.8.4].

Lieb’s Theorem [Lie73, Thm. 6] was first established in an important paper of Elliott Lieb on the convexity of trace functions. His main goal was to establish concavity properties for a function that measures the amount of information in a quantum system. See the notes in Chapter 8 for a more detailed discussion.

Chapter 4 Matrix Gaussian Series & Matrix Rademacher Series

In this chapter, we present our first set of matrix concentration inequalities. These results provide spectral information about a sum of fixed matrices, each modulated by an independent scalar random variable. This type of formulation is surprisingly versatile, and it captures a range of interesting examples. Our main goal, however, is to introduce matrix concentration in the simplest setting possible.

To be more precise about our scope, let us introduce the concept of a matrix Gaussian series. Consider a finite sequence $\{\bm{B}_{k}\}$ of fixed matrices with the same dimension, along with a finite sequence $\{\gamma_{k}\}$ of independent standard normal random variables. We will study the spectral norm of the random matrix

This expression looks abstract, but it has concrete modeling power. For example, we can express a Gaussian Wigner matrix, one of the classical random matrices, in this fashion. But the real value of this approach is that we can use matrix Gaussian series to represent many kinds of random matrices built from Gaussian random variables. This technique allows us to attack problems that classical methods do not handle gracefully. For instance, we can easily study a Toeplitz matrix with Gaussian entries.

Similar ideas allow us to treat a matrix Rademacher series, a sum of fixed matrices modulated by random signs. (Recall that a Rademacher random variable takes the values $\pm 1$ with equal probability.) The results in this case are almost identical with the results for matrix Gaussian series, but they allow us to consider new problems. As an example, we can study the expected spectral norm of a fixed real matrix after flipping the signs of the entries at random.

In §1, we begin with an overview of our results for matrix Gaussian series; very similar results also hold for matrix Rademacher series. Afterward, we discuss the accuracy of the theoretical bounds. The subsequent sections, §§2–4, describe what the matrix concentration inequalities tell us about some classical and not-so-classical examples of random matrices. Section 5 includes an overview of a more substantial application in combinatorial optimization. The final part §6 contains detailed proofs of the bounds. We conclude with bibliographical notes.

A Norm Bound for Random Series with Matrix Coefficients

Consider a finite sequence $\{b_{k}\}$ of real numbers and a finite sequence $\{\gamma_{k}\}$ of independent standard normal random variables. Form the random series $Z=\sum_{k}\gamma_{k}b_{k}$. A routine invocation of the scalar Laplace transform method demonstrates that

It turns out that the inequality (1) extends directly to the matrix setting.

Consider a finite sequence $\{\bm{B}_{k}\}$ of fixed complex matrices with dimension $d_{1}\times d_{2}$, and let $\{\gamma_{k}\}$ be a finite sequence of independent standard normal variables. Introduce the matrix Gaussian series

Let $v(\bm{Z})$ be the matrix variance statistic of the sum:

The same bounds hold when we replace $\{\gamma_{k}\}$ by a finite sequence $\{\varrho_{k}\}$ of independent Rademacher random variables.

The proof of Theorem 1.1 appears below in §6.

Let us take a moment to discuss the content of Theorem 1.1. The main message is that the expectation of $\left\|{\bm{Z}}\right\|$ is controlled by the matrix variance statistic $v(\bm{Z})$. Furthermore, $\left\|{\bm{Z}}\right\|$ has a subgaussian tail whose decay rate depends on $v(\bm{Z})$.

The matrix variance statistic $v(\bm{Z})$ defined in (3) specializes the general formulation (8). The second expression (4) follows from the additivity property (11) for the variance of an independent sum. When the summands are Hermitian, observe that the two terms in the maximum coincide. The formulas (3) and (4) are a direct extension of the variance that arises in the scalar bound (1).

As compared with (1), a new feature of the bound (6) is the dimensional factor $d_{1}+d_{2}$. When $d_{1}=d_{2}=1$, the matrix bound reduces to the scalar result (1). In this case, at least, we have lost nothing by lifting the Laplace transform method to matrices. The behavior of the matrix tail bound (6) is more subtle than the behavior of the scalar tail bound (1). See Figure 1 for an illustration.

Optimality of the Bounds for Matrix Gaussian Series

One may wonder whether Theorem 1.1 provides accurate information about the behavior of a matrix Gaussian series. The answer turns out to be complicated. Here is the executive summary: the expectation bound (5) is always quite good, but the tail bound (6) is sometimes quite bad. The rest of this section expands on these claims.

Let $\bm{Z}$ be a matrix Gaussian series of the form (2). We will argue that

In other words, the matrix variance $v(\bm{Z})$ is roughly the correct scale for $\left\|{\bm{Z}}\right\|^{2}$. This pair of estimates is a significant achievement because it is quite challenging to compute the norm of a matrix Gaussian series in general. Indeed, the literature contains very few examples where explicit estimates are available, especially if one desires reasonable constants.

We begin with the lower bound in (7), which is elementary. Indeed, since the spectral norm is convex, Jensen’s inequality ensures that

The first identity follows from (24), and the last is the definition (8) of the matrix variance.

The upper bound in (7) is a consequence of the tail bound (6):

In the first step, rewrite the expectation using integration by parts, and then split the integral at a positive number $E$. In the first term, we bound the probability by one, while the second term results from the tail bound (6). Afterward, we compute the integrals explicitly. Finally, select $E^{2}=2v(\bm{Z})\log(d_{1}+d_{2})$ to complete the proof of (7).

At this point, one may ask whether it is possible to improve either side of the inequality (7). The answer is negative unless we have additional information about the Gaussian series beyond the matrix variance statistic $v(\bm{Z})$.

In Chapter 7, we will describe a technique that allows us to moderate the dimensional factor in (7) for some types of matrix series. But we cannot remove the dimensional factor entirely with current technology.

What about the tail bound (6) for the norm of the Gaussian series? Here, our results are less impressive. It turns out that the large-deviation behavior of the spectral norm of a matrix Gaussian series $\bm{Z}$ is controlled by a statistic $v_{\star}(\bm{Z})$ called the weak variance:

The best general inequalities between the matrix variance statistic and the weak variance are

There are examples of matrix Gaussian series that saturate the lower or the upper inequality.

The classical concentration inequality [BLM13, Thm. 5.6] for a function of independent Gaussian random variables implies that

When studying concentration of random variables, it is quite common that we need to use one method to assess the expected value of the random variable and a separate technique to determine the probability of a large deviation.

The primary value of matrix concentration inequalities inheres in the estimates that they provide for the expectation of the spectral norm (or maximum eigenvalue or minimum eigenvalue) of a random matrix.

In many cases, matrix concentration bounds provide reasonable information about the tail decay, but there are other situations where the tail bounds are feeble. In this event, we recommend applying a scalar concentration inequality to control the tails.

Example: Some Gaussian Matrices

Let us try out our methods on two types of Gaussian matrices that have been studied extensively in the classical literature on random matrix theory. In these cases, precise information about the spectral distribution is available, which provides a benchmark for assessing our results. We find that bounds based on Theorem 1.1 lead to very reasonable estimates, but they are not sharp. The advantage of our approach is that it applies to every example, whereas we are making comparisons with specialized techniques that only illuminate individual cases. Similar conclusions hold for matrices with independent Rademacher entries.

We begin with a family of Gaussian Wigner matrices. A $d\times d$ matrix $\bm{W}_{d}$ from this ensemble is real-symmetric with a zero diagonal; the entries above the diagonal are independent normal variables with mean zero and variance one:

where $\{\gamma_{jk}:1\leq j<k\leq d\}$ is an independent family of standard normal variables. We can represent this matrix compactly as a Gaussian series:

For example, see [BS10, Thm. 5.1]. To make (2) precise, we assume that $\{\bm{W}_{d}\}$ is an independent sequence of Gaussian Wigner matrices, indexed by the dimension $d$.

Theorem 6.1 provides a simple way to bound the norm of a Gaussian Wigner matrix. We just need to compute the matrix variance statistic $v(\bm{W}_{d})$. The formula (4) for $v(\bm{W}_{d})$ asks us to form the sum of the squared coefficients from the representation (1):

Since the terms in (1) are Hermitian, we have only one sum of squares to consider. We have also used the facts that $\mathbf{E}_{jk}\mathbf{E}_{kj}=\mathbf{E}_{jj}$ while $\mathbf{E}_{jk}\mathbf{E}_{jk}=\bm{0}$ because of the condition $j<k$ in the limits of summation. We see that

The bound (5) for the expectation of the norm gives

In conclusion, our techniques overestimate $\left\|{\bm{W}_{d}}\right\|$ by a factor of about $\sqrt{0.5\log d}$. The result (3) is not perfect, but it only takes two lines of work. In contrast, the classical result (2) depends on a long moment calculation that involves challenging combinatorial arguments.

Rectangular Gaussian Matrices

Next, we consider a $d_{1}\times d_{2}$ rectangular matrix with independent standard normal entries:

where $\{\gamma_{jk}\}$ is an independent family of standard normal variables. We can express this matrix efficiently using a Gaussian series:

There is an elegant estimate [DS02, Thm. 2.13] for the norm of this matrix:

Theorem 1.1 yields another bound on the expected norm of the matrix $\bm{G}$. In order to compute the matrix variance statistic $v(\bm{G})$, we calculate the sums of the squared coefficients from the representation (4):

The matrix variance statistic (3) satisfies

The leading term is roughly correct because

The logarithmic factor in (6) does not belong, but it is rather small in comparison with the leading terms. Once again, we have produced a reasonable result with a short argument based on general principles.

Example: Matrices with Randomly Signed Entries

Next, we turn to an example that is superficially similar with the matrix discussed in §2 but is less understood. Consider a fixed $d_{1}\times d_{2}$ matrix $\bm{B}$ with real entries, and let $\{\varrho_{jk}\}$ be an independent family of Rademacher random variables. Consider the $d_{1}\times d_{2}$ random matrix

In other words, we obtain the random matrix $\bm{B}_{\pm}$ by randomly flipping the sign of each entry of $\bm{B}$. The expected norm of this matrix satisfies the bound

where the leading factor $v^{1/2}$ satisfies

Theorem 1.1 leads to a quick proof of a slightly weaker result. We simply need to compute the matrix variance statistic $v(\bm{B}_{\pm})$. To that end, note that

Therefore, using the formula (4), we find that

We see that $v(\bm{B}_{\pm})$ coincides with $v$, the leading term (2) in the established estimate (1)! Now, Theorem 1.1 delivers the bound

Observe that the estimate (3) for the norm matches the correct bound (1) up to the logarithmic factor. Yet again, we obtain a result that is respectably close to the optimal one, even though it is not quite sharp.

The main advantage of using results like Theorem 1.1 to analyze this random matrix is that we can obtain a good result with a minimal amount of arithmetic. The analysis that leads to (1) involves a specialized combinatorial argument.

Example: Gaussian Toeplitz Matrices

Matrix concentration inequalities offer an effective tool for analyzing random matrices whose dependency structures are more complicated than those of the classical ensembles. In this section, we consider Gaussian Toeplitz matrices, which have applications in signal processing.

We construct an (unsymmetric) $d\times d$ Gaussian Toeplitz matrix $\bm{\Gamma}_{d}$ by populating the first row and first column of the matrix with independent standard normal variables; the entries along each diagonal of the matrix take the same value:

where $\{\gamma_{k}\}$ is an independent family of standard normal variables. As usual, we represent the Gaussian Toeplitz matrix as a matrix Gaussian series:

It follows that $\bm{C}^{k}$ shifts a vector up by $k$ places, introducing zeros at the bottom, while $(\bm{C}^{k})^{*}$ shifts a vector down by $k$ places, introducing zeros at the top.

We can analyze this example quickly using Theorem 1.1. First, note that

To obtain the matrix variance statistic (4), we calculate the sum of the squares of the coefficient matrices that appear in (1). In this instance, the two terms in the variance are the same. We find that

In the second line, we (carefully) switch the order of summation and rewrite the identity matrix as a sum of diagonal standard basis matrices. We reach

An application of Theorem 1.1 leads us to conclude that

It turns out that the inequality (3) is correct up to the precise value of the constant, which does not seem to be known. Nevertheless, the limiting value is available for the top eigenvalue of a (scaled) symmetric Toeplitz matrix whose first row contains independent standard normal variables [SV13, Thm. 1]. From this result, we may conclude that

Here, we take $\{\bm{\Gamma}_{d}\}$ to be a sequence of unsymmetric Gaussian Toeplitz matrices, indexed by the ambient dimension $d$. Our simple argument gives the right scaling for this problem, and our estimate for the constant lies within 21% of the optimal value!

Application: Rounding for the MaxQP Relaxation

Our final application involves a more substantial question from combinatorial optimization. One of the methods that has been proposed for solving a certain optimization problem leads to a matrix Rademacher series, and the analysis of this method requires the spectral norm bounds from Theorem 1.1. A detailed treatment would take us too far afield, so we just sketch the context and indicate how the random matrix arises.

There are many types of optimization problems that are computationally difficult to solve exactly. One approach to solving these problems is to enlarge the constraint set in such a way that the problem becomes tractable, a process called “relaxation.” After solving the relaxed problem, we can use a randomized “rounding” procedure to map the solution back to the constraint set for the original problem. If we can perform the rounding step without changing the value of the objective function substantially, then the rounded solution is also a decent solution to the original optimization problem.

One difficult class of optimization problems has a matrix decision variable, and it requires us to maximize a quadratic form in the matrix variable subject to a set of convex quadratic constraints and a spectral norm constraint [Nem07]. This problem is referred to as MaxQP. The desired solution $\bm{B}$ to this problem is a $d_{1}\times d_{2}$ matrix. The solution needs to satisfy several different requirements, but we focus on the condition that $\left\|{\bm{B}}\right\|\leq 1$.

There is a natural relaxation of the MaxQP problem. When we solve the relaxation, we obtain a family $\{\bm{B}_{k}:k=1,2,\dots,n\}$ of $d_{1}\times d_{2}$ matrices that satisfy the constraints

In fact, these two bounds are part of the specification of the relaxed problem. To round the family of matrices back to a solution of the original problem, we form the random matrix

where $\{\varrho_{k}\}$ is an independent family of Rademacher random variables. The scaling factor $\alpha>0$ can be adjusted to guarantee that the norm constraint $\left\|{\bm{Z}}\right\|\leq 1$ holds with high probability.

What is the expected norm of $\bm{Z}$? Theorem 1.1 yields

Here, the matrix variance statistic satisfies

owing to the constraint (1) on the matrices $\bm{B}_{1},\dots,\bm{B}_{n}$. It follows that the scaling parameter $\alpha$ should satisfy

The important fact here is that the scaling parameter $\alpha$ is usually small as compared with the other parameters of the problem ($d_{1},d_{2}$, $n$, and so forth). Therefore, the scaling does not have a massive effect on the value of the objective function. Ultimately, this approach leads to a technique for solving the MaxQP problem that produces a feasible point whose objective value is within a factor of $\sqrt{2\log(d_{1}+d_{2})}$ of the maximum objective value possible.

Analysis of Matrix Gaussian & Rademacher Series

We began this chapter with a concentration inequality, Theorem 1.1, for the norm of a matrix Gaussian series, and we have explored a number of different applications of this result. This section contains a proof of this theorem.

As the development in Chapter 3 suggests, random Hermitian matrices provide the natural setting for establishing matrix concentration inequalities. Therefore, we begin our treatment with a detailed statement of the matrix concentration inequality for a Gaussian series with Hermitian matrix coefficients.

Consider a finite sequence $\{\bm{A}_{k}\}$ of fixed Hermitian matrices with dimension $d$, and let $\{\gamma_{k}\}$ be a finite sequence of independent standard normal variables. Introduce the matrix Gaussian series

Let $v(\bm{Y})$ be the matrix variance statistic of the sum:

The same bounds hold when we replace $\{\gamma_{k}\}$ by a finite sequence of independent Rademacher random variables.

The proof of this result occupies the rest of the section.

Discussion

Before we proceed to the analysis, let us take a moment to compare Theorem 6.1 with the result for general matrix series, Theorem 1.1.

First, we consider the matrix variance statistic $v(\bm{Y})$ defined in (1). Since $\bm{Y}$ has zero mean, this definition coincides with the general formula (4). The second expression, in terms of the coefficient matrices, follows from the additivity property (6) for the variance of a sum of independent, random Hermitian matrices.

Next, bounds for the minimum eigenvalue $\lambda_{\min}(\bm{Y})$ follow from the results for the maximum eigenvalue because $-\bm{Y}$ has the same distribution as $\bm{Y}$. Therefore,

The second identity holds because of the relationship (5) between minimum and maximum eigenvalues. Similar considerations lead to a lower tail bound for the minimum eigenvalue:

This result follows directly from the upper tail bound (3).

This observation points to the most important difference between the Hermitian case and the general case. Indeed, Theorem 6.1 concerns the extreme eigenvalues of the random series $\bm{Y}$ instead of the norm. This change amounts to producing one-sided tail bounds instead of two-sided tail bounds. For Gaussian and Rademacher series, this improvement is not really useful, but there are random Hermitian matrices whose minimum and maximum eigenvalues exhibit different types of behavior. For these problems, it can be extremely valuable to examine the two tails separately. See Chapter 5 and 6 for some results of this type.

Analysis for Hermitian Gaussian Series

We continue with the proof that matrix Gaussian series exhibit the behavior described in Theorem 6.1. Afterward, we show how to adapt the argument to address matrix Rademacher series. Our main tool is Theorem 6.1, the set of master bounds for independent sums. To use this result, we must identify the cgf of a fixed matrix modulated by a Gaussian random variable.

Suppose that $\bm{A}$ is a fixed Hermitian matrix, and let $\gamma$ be a standard normal random variable. Then

We may assume $\theta=1$ by absorbing $\theta$ into the matrix $\bm{A}$. It is well known that the moments of a standard normal variable satisfy

The formula for the odd moments holds because a standard normal variable is symmetric. One way to establish the formula for the even moments is to use integration by parts to obtain a recursion for the $(2q)$th moment in terms of the $(2q-2)$th moment.

The first identity holds because the odd terms vanish from the series representation (15) of the matrix exponential when we take the expectation. To compute the cgf, we extract the logarithm of the mgf and recall (17), which states that the matrix logarithm is the functional inverse of the matrix exponential. ∎

We quickly reach results on the maximum eigenvalue of a matrix Gaussian series with Hermitian coefficients.

Consider a finite sequence $\{\bm{A}_{k}\}$ of Hermitian matrices with dimension $d$, and let $\{\gamma_{k}\}$ be a finite sequence of independent standard normal variables. Define the matrix Gaussian series

The second line follows when we introduce the cgf from Lemma 6.2. To reach the third inequality, we bound the trace by the dimension times the maximum eigenvalue. The fourth line is the Spectral Mapping Theorem, Proposition 1.3. Use the formula (1) to identify the matrix variance statistic $v(\bm{Y})$ in the exponent. The infimum is attained at $\theta=\sqrt{2v(\bm{Y})^{-1}\log d}$. This choice leads to (2).

Next, we turn to the proof of the upper tail bound (3) for $\lambda_{\max}(\bm{Y})$. Invoke the master tail bound (3) from Theorem 6.1, and calculate that

The steps here are the same as in the previous calculation. The infimum is achieved at $\theta=t/v(\bm{Y})$, which yields (3). ∎

Analysis for Hermitian Rademacher Series

The inequalities for matrix Rademacher series involve arguments closely related to the proofs for matrix Gaussian series, but we require one additional piece of reasoning to obtain the simplest results. First, let us compute bounds for the matrix mgf and cgf of a Hermitian matrix modulated by a Rademacher random variable.

Suppose that $\bm{A}$ is a fixed Hermitian matrix, and let $\varrho$ be a Rademacher random variable. Then

First, we establish a scalar inequality. Comparing Taylor series,

The inequality holds because $(2q)!\geq(2q)(2q-2)\cdots(4)(2)=2^{q}q!$.

To compute the matrix mgf, we may assume $\theta=1$. By direct calculation,

The semidefinite bound follows when we apply the Transfer Rule (14) to the inequality (6).

To determine the matrix cgf, observe that

We are prepared to develop some probability inequalities for the maximum eigenvalue of a Rademacher series with Hermitian coefficients.

Consider a finite sequence $\{\bm{A}_{k}\}$ of Hermitian matrices, and let $\{\varrho_{k}\}$ be a finite sequence of independent Rademacher variables. Define the matrix Rademacher series

The bounds for the extreme eigenvalues of $\bm{Y}$ follow from an argument almost identical with the proof in the Gaussian case. The only point that requires justification is the inequality

To obtain this result, we introduce the semidefinite bound, Lemma 6.3, for the Rademacher cgf into the trace exponential. The left-hand side increases after this substitution because of the fact (16) that the trace exponential function is monotone with respect to the semidefinite order. ∎

Analysis of Matrix Series with Rectangular Coefficients

Finally, we consider a series with non-Hermitian matrix coefficients modulated by independent Gaussian or Rademacher random variables. The bounds for the norm of a rectangular series follow instantly from the bounds for the norm of an Hermitian series because of a formal device. We simply apply the Hermitian results to the Hermitian dilation (26) of the series.

Consider a finite sequence $\{\bm{B}_{k}\}$ of $d_{1}\times d_{2}$ complex matrices, and let $\{\zeta_{k}\}$ be a finite sequence of independent random variables, either standard normal or Rademacher.

Recall from Definition 1.5 that the Hermitian dilation is the map

The second expression for $\bm{Y}$ holds because the Hermitian dilation is real-linear. Since we have written $\bm{Y}$ as a matrix series with Hermitian coefficients, we may analyze it using Theorem 6.1. We just need to express the conclusions in terms of the random matrix $\bm{Z}$.

First, we employ the fact (28) that the Hermitian dilation preserves spectral information:

Therefore, bounds on $\lambda_{\max}(\bm{Y})$ deliver bounds on $\left\|{\bm{Z}}\right\|$. In view of the calculation (10) for the variance statistic of a dilation, we have

Recall that the matrix variance statistic $v(\bm{Z})$ defined in (3) coincides with the general definition from (8). Now, invoke Theorem 6.1 to obtain Theorem 1.1. ∎

Notes

We give an overview of research related to matrix Gaussian series, along with references for the specific random matrices that we have analyzed.

The main results, Theorem 1.1 and Theorem 6.1, have an interesting history. In the precise form presented here, these two statements first appeared in [Tro11c], but we can trace them back more than two decades.

In his work [Oli10b, Thm. 1], Oliveira established the mgf bounds presented in Lemma 6.2 and Lemma 6.3. He also developed an ingenious improvement on the arguments of Ahlswede & Winter [AW02, App.], and he obtained a bound similar with Theorem 6.1. The constants in Oliveira’s result are worse, but the dependence on the dimension is better because it depends on the number of summands. We do not believe that the approach Ahlswede & Winter describe in [AW02] can deliver any of these results.

Recently, there have been some minor improvements to the dimensional factor that appears in Theorem 6.1. We discuss these results and give citations in Chapter 7.

The Noncommutative Khintchine Inequality

Our theory about matrix Rademacher and Gaussian series should be compared with a classic result, called the noncommutative Khintchine inequality, that was originally due to Lust-Piquard [LP86]; see also the follow-up work [LPP91]. In its simplest form, this inequality concerns a matrix Rademacher series with Hermitian coefficients:

The noncommutative Khintchine inequality states that

The minimum value of the constant $C_{2q}=(2q)!/(2^{q}\,q!)$ was obtained in the two papers [Buc01, Buc05]. Traditional proofs of the noncommutative Khintchine inequality are quite involved, but there is now an elementary argument available [MJC+14, Cor. 7.3].

Theorem 6.1 is the exponential moment analog of the polynomial moment bound (1). The polynomial moment inequality is somewhat stronger than the exponential moment inequality. Nevertheless, the exponential results are often more useful in practice. For a more thorough exploration of the relationships between Theorem 6.1 and noncommutative moment inequalities, such as (1), see the discussion in [Tro11c, §4].

Application to Random Matrices

It has also been known for a long time that results such as Theorem 6.1 and inequality (1) can be used to study random matrices.

We believe that the geometric functional analysis literature contains the earliest applications of matrix concentration results to analyze random matrices. In a well-known paper [Rud99], Mark Rudelson—acting on a suggestion of Gilles Pisier—showed how to use the noncommutative Khintchine inequality (1) to study covariance estimation. This work led to a significant amount of activity in which researchers used variants of Rudelson’s argument to prove other types of results. See, for example, the paper [RV07]. This approach is powerful, but it tends to require some effort to use.

In parallel, other researchers in noncommutative probability theory also came to recognize the power of noncommutative moment inequalities in random matrix theory. The paper [JX08] contains a specific example. Unfortunately, this literature is technically formidable, which makes it difficult for outsiders to appreciate its achievements.

The work [AW02] of Ahlswede & Winter led to the first “packaged” matrix concentration inequalities of the type that we describe in these lecture notes. For the first few years after this work, most of the applications concerned quantum information theory and random graph theory. The paper [Gro11] introduced the method of Ahlswede & Winter to researchers in mathematical signal processing and statistics, and it served to popularize matrix concentration bounds.

At this point, the available matrix concentration inequalities were still significantly suboptimal. The main advances, in [Oli10a, Tro11c], led to optimal matrix concentration results of the kind that we present in these lecture notes. These results allow researchers to obtain reasonably accurate analyses of a wide variety of random matrices with very little effort.

Wigner and Marčenko–Pastur

Wigner matrices first emerged in the literature on nuclear physics, where they were used to model the Hamiltonians of reactions involving heavy atoms [Meh04, §1.1]. Wigner [Wig55] showed that the limiting spectral distribution of a certain type of Wigner matrix follows the semicircle law. See the book [Tao12, §2.4] of Tao for an overview and the book [BS10, Chap. 2] of Bai & Silverstein for a complete treatment. The Bai–Yin law [BY93] states that, up to scaling, the maximum eigenvalue of a Wigner matrix converges almost surely to two. See [Tao12, §2.3] or [BS10, Chap. 5] for more information. The analysis of the Gaussian Wigner matrix that we present here, using Theorem 6.1, is drawn from [Tro11c, §4].

The first rigorous work on a rectangular Gaussian matrix is due to Marčenko & Pastur [MP67], who established that the limiting distribution of the squared singular values follows a distribution that now bears their names. The Bai–Yin law [BY93] gives an almost-sure limit for the largest singular value of a rectangular Gaussian matrix. The expectation bound (5) appears in a survey article [DS02] by Davidson & Szarek. The latter result is ultimately derived from a comparison theorem for Gaussian processes due to Férnique [Fer75] and amplified by Gordon [Gor85]. Our approach, using Theorem 1.1, is based on [Tro11c, §4].

Randomly Signed Matrices

Matrices with randomly signed entries have not received much attention in the literature. The result (1) is due to Yoav Seginer [Seg00]. There is also a well-known paper [Lat05] by Rafał Latała that provides a bound for the expected norm of a Gaussian matrix whose entries have nonuniform variance. Riemer & Schütt [RS13] have extended the earlier results. The very recent paper [BV14] of Afonso Bandeira and Ramon Van Handel contains an elegant new proof of Seginer’s result based on a general theorem for random matrices with independent entries. The analysis here, using Theorem 1.1, is drawn from [Tro11c, §4].

Gaussian Toeplitz Matrices

Relaxation and Rounding of MaxQP

The idea of using semidefinite relaxation and rounding to solve the MaxQP problem is due to Arkadi Nemirovski [Nem07]. He obtained nontrivial results on the performance of his method using some matrix moment calculations, but he was unable to reach the sharpest possible bound. Anthony So [So09] pointed out that matrix moment inequalities imply an optimal result; he also showed that matrix concentration inequalities have applications to robust optimization. The presentation here, using Theorem 1.1, is essentially equivalent with the approach in [So09], but we have achieved slightly better bounds for the constants.

Chapter 5 A Sum of Random Positive-Semidefinite Matrices

This chapter presents matrix concentration inequalities that are analogous with the classical Chernoff bounds. In the matrix setting, Chernoff-type inequalities allow us to control the extreme eigenvalues of a sum of independent, random, positive-semidefinite matrices.

More formally, we consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random Hermitian matrices that satisfy

Introduce the sum $\bm{Y}=\sum_{k}\bm{X}_{k}$. Our goal is to study the expectation and tail behavior of $\lambda_{\max}(\bm{Y})$ and $\lambda_{\min}(\bm{Y})$. Bounds on the maximum eigenvalue $\lambda_{\max}(\bm{Y})$ give us information about the norm of the matrix $\bm{Y}$, a measure of how much the action of the matrix can dilate a vector. Bounds for the minimum eigenvalue $\lambda_{\min}(\bm{Y})$ tell us when the matrix $\bm{Y}$ is nonsingular; they also provide evidence about the norm of the inverse $\bm{Y}^{-1}$, when it exists.

The matrix Chernoff inequalities are quite powerful, and they have numerous applications. We demonstrate the relevance of this theory by considering two examples. First, we show how to study the norm of a random submatrix drawn from a fixed matrix, and we explain how to check when the random submatrix has full rank. Second, we develop an analysis to determine when a random graph is likely to be connected. These two problems are closely related to basic questions in statistics and in combinatorics.

In contrast, the matrix Bernstein inequalities, appearing in Chapter 6, describe how much a random matrix deviates from its mean value. As such, the matrix Bernstein bounds are more suitable than the matrix Chernoff bounds for problems that concern matrix approximations. Matrix Bernstein inequalities are also more appropriate when the variance $v(\bm{Y})$ is small in comparison with the upper bound $L$ on the summands.

Section 1 presents the main results on the expectations and the tails of the extreme eigenvalues of a sum of independent, random, positive-semidefinite matrices. Section 2 explains how the matrix Chernoff bounds provide spectral information about a random submatrix drawn from a fixed matrix. In §3, we use the matrix Chernoff bounds to study when a random graph is connected. Afterward, in §4 we explain how to prove the main results.

The Matrix Chernoff Inequalities

In the scalar setting, the Chernoff inequalities describe the behavior of a sum of independent, nonnegative random variables that are subject to a uniform upper bound. These results are often applied to study the number $Y$ of successes in a sequence of independent—but not necessarily identical—Bernoulli trials with small probabilities of success. In this case, the Chernoff bounds show that $Y$ behaves like a Poisson random variable. The random variable $Y$ concentrates near the expected number of successes. Its lower tail has Gaussian decay, while its upper tail drops off faster than that of an exponential random variable. See [BLM13, §2.2] for more background.

In the matrix setting, we encounter similar phenomena when we consider a sum of independent, random, positive-semidefinite matrices whose eigenvalues meet a uniform upper bound. This behavior emerges from the next theorem, which closely parallels the scalar Chernoff theorem.

Consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random, Hermitian matrices with common dimension $d$. Assume that

The proof of Theorem 1.1 appears below in §4.

Let us consider some facets of Theorem 1.1.

In many situations, it is easier to work with streamlined versions of the expectation bounds:

We obtain these results by selecting $\theta=1$ in both (3) and (4) and evaluating the numerical constants.

We can also weaken the tail bounds (5) and (6) to reach

The first bound shows that the lower tail of $\lambda_{\min}(\bm{Y})$ decays at a subgaussian rate with variance $L/\mu_{\min}$. The second bound manifests that the upper tail of $\lambda_{\max}(\bm{Y})$ decays faster than that of an exponential random variable with mean $L/\mu_{\max}$. This is the same type of prediction we receive from the scalar Chernoff inequalities.

As with other matrix concentration results, the tail bounds (5) and (6) can overestimate the actual tail probabilities for the extreme eigenvalues of $\bm{Y}$, especially at large deviations from the mean. The value of the matrix Chernoff theorem derives from the estimates (3) and (4) for the expectation of the minimum and maximum eigenvalue of $\bm{Y}$. Scalar concentration inequalities may provide better estimates for tail probabilities.

We can moderate the dimensional factor $d$ in the bounds for $\lambda_{\max}(\bm{Y})$ from Theorem 1.1 when the random matrix $\bm{Y}$ has limited spectral content in most directions. We take up this analysis in Chapter 7.

Next, let us present an important refinement [CGT12a, Thm. A.1] of the bound (8) that can be very useful in practice:

This estimate may be regarded as a matrix version of Rosenthal’s inequality [Ros70]. Observe that the uniform bound $L$ appearing in (8) always exceeds the large parenthesis on the right-hand side of (9). Therefore, the estimate (9) is valuable when the summands are unbounded and, especially, when they have heavy tails. See the notes at the end of the chapter for more information.

Optimality of the Matrix Chernoff Bounds

In this section, we explore how well bounds such as Theorem 1.1 and inequality (9) describe the behavior of a random matrix $\bm{Y}$ formed as a sum of independent, random positive-semidefinite matrices.

We will demonstrate that both terms in the matrix Rosenthal inequality (9) are necessary. More precisely,

The appearance of $\mu_{\max}$ on the left-hand side of (10) is a consequence of Jensen’s inequality. Indeed, the maximum eigenvalue is convex, so

To justify the other term, apply the fact that the summands $\bm{X}_{k}$ are positive semidefinite to conclude that

We have used the fact that $\lambda_{\max}(\bm{A}+\bm{H})\geq\lambda_{\max}(\bm{A})$ whenever $\bm{H}$ is positive semidefinite. Average the last two displays to develop the left-hand side of (10). The right-hand side of (10) is obviously just (9).

A simple example suffices to show that the logarithm cannot always be removed from the second term in (8) or from (9). For each natural number $n$, consider the $d\times d$ random matrix

where $\big\{\delta_{ik}^{(n)}\big\}$ is an independent family of $\textsc{bernoulli}(n^{-1})$ random variables and $\mathbf{E}_{kk}$ is the $d\times d$ matrix with a one in the $(k,k)$ entry an zeros elsewhere. An easy application of (8) delivers

Using the Poisson limit of a binomial random variable and the Skorokhod representation, we can construct an independent family $\{Q_{k}\}$ of $\textsc{poisson}(1)$ random variables for which

Therefore, the logarithm on the second term in (8) cannot be reduced by a factor larger than the iterated logarithm $\log\log d$. This modest loss comes from approximations we make when developing the estimate for the mean. The tail bound (6) accurately predicts the order of $\lambda_{\max}(\bm{Y}_{n})$ in this example.

The latter example depends on the commutativity of the summands and the infinite divisibility of the Poisson distribution, so it may seem rather special. Nevertheless, the logarithm really does belong in many (but not all!) examples that arise in practice. In particular, it is necessary in the application to random submatrices in §2.

The upper expectation bound (4) is quite satisfactory, but the situation is murkier for the lower expectation bound (3). The mean term appears naturally in the lower bound:

This estimate is a consequence of Jensen’s inequality and the concavity of the minimum eigenvalue. On the other hand, it is not clear what the correct form of the second term in (3) should be for a general sum of random positive-semidefinite matrices.

Nevertheless, a simple example demonstrates that the lower Chernoff bound (3) is numerically sharp in some situations. Let $\bm{X}$ be a $d\times d$ random positive-semidefinite matrix that satisfies

The lower Chernoff bound (3) implies that

On the other hand, $\lambda_{\max}(\bm{Y}_{n})>0$ if and only if each diagonal matrix $d\,\mathbf{E}_{ii}$ appears at least once among the summands $\bm{X}_{1},\dots,\bm{X}_{n}$. To determine the probability that this event occurs, notice that this question is an instance of the coupon collector problem [MR95, §3.6]. The probability of collecting all $d$ coupons within $n$ draws undergoes a phase transition from about zero to about one at $n=d\log d$. By refining this argument [Tro11d], we can verify that both lower Chernoff bounds (3) and (5) provide a numerically sharp lower bound for the value of $n$ where the phase transition occurs. In other words, the lower matrix Chernoff bounds are themselves sharp.

Example: A Random Submatrix of a Fixed Matrix

The matrix Chernoff inequality can be used to bound the extreme singular values of a random submatrix drawn from a fixed matrix. Theorem 1.1 might not seem suitable for this purpose because it deals with eigenvalues, but we can connect the method with the problem via a simple transformation. The results in this section have found applications in randomized linear algebra, sparse approximation, machine learning, and other fields. See the notes at the end of the chapter for some additional discussion and references.

Let $\bm{B}$ be a fixed $d\times n$ matrix, and let $\bm{b}_{:k}$ denote the $k$th column of this matrix. The matrix can be expressed as the sum of its columns:

The symbol $\mathbf{e}_{k}$ refers to the standard basis (column) vector with a one in the $k$th component and zeros elsewhere; the length of the vector $\mathbf{e}_{k}$ is determined by context.

We consider a simple model for a random column submatrix. Let $\{\delta_{k}\}$ be an independent sequence of $\textsc{bernoulli}(p/n)$ random variables. Define the random matrix

That is, we include each column independently with probability $p/n$, which means that there are typically about $p$ nonzero columns in the matrix. We do not remove the other columns; we just zero them out.

In this section, we will obtain bounds on the expectation of the extreme singular values $\sigma_{1}(\bm{Z})$ and $\sigma_{d}(\bm{Z})$ of the $d\times n$ random matrix $\bm{Z}$. More precisely,

That is, the random submatrix $\bm{Z}$ gets its “fair share” of the squared singular values of the original matrix $\bm{B}$. There is a fluctuation term that depends on largest norm of a column of $\bm{B}$ and the logarithm of the number $d$ of rows in $\bm{B}$. This result is very useful because a positive lower bound on $\sigma_{d}(\bm{Z})$ ensures that the rows of the random submatrix $\bm{Z}$ are linearly independent.

To study the singular values of $\bm{Z}$, it is convenient to define a $d\times d$ random, positive-semidefinite matrix

Note that $\delta_{k}^{2}=\delta_{k}$ because $\delta_{k}$ only takes the values zero and one. The eigenvalues of $\bm{Y}$ determine the singular values of $\bm{Z}$, and vice versa. In particular,

where we arrange the singular values of $\bm{Z}$ in weakly decreasing order $\sigma_{1}(\bm{Z})\geq\dots\geq\sigma_{d}(\bm{Z})$.

The matrix Chernoff inequality provides bounds for the expectations of the eigenvalues of $\bm{Y}$. To apply the result, first calculate

Define $L=\max_{k}\left\|{\bm{b}_{:k}}\right\|^{2}$, and observe that $\left\|{\smash{\delta_{k}\bm{b}_{:k}\bm{b}_{:k}^{*}}}\right\|\leq L$ for each index $k$. The simplified matrix Chernoff bounds (7) and (8) now deliver the result (1).

A Random Row and Column Submatrix

Next, we consider a model for a random set of rows and columns drawn from a fixed $d\times n$ matrix $\bm{B}$. In this case, it is helpful to use matrix notation to represent the extraction of a submatrix. Define independent random projectors

where $\{\delta_{k}\}$ is an independent family of $\textsc{bernoulli}(p/d)$ random variables and $\{\xi_{k}\}$ is an independent family of $\textsc{bernoulli}(r/n)$ random variables. Then

is a random submatrix of $\bm{B}$ with about $p$ nonzero rows and $r$ nonzero columns.

The notations $\bm{b}_{j:}$ and $\bm{b}_{:k}$ refer to the $j$th row and $k$th column of the matrix $\bm{B}$, while $b_{jk}$ is the $(j,k)$ entry of the matrix. In other words, the random submatrix $\bm{Z}$ gets its share of the total squared norm of the matrix $\bm{B}$. The fluctuation terms reflect the maximum row norm and the maximum column norm of $\bm{B}$, as well as the size of the largest entry. There is also a weak dependence on the ambient dimensions $d$ and $n$.

The argument has much in common with the calculations for a random column submatrix, but we need to do some extra work to handle the interaction between the random row sampling and the random column sampling.

To begin, we express the squared norm $\left\|{\bm{Z}}\right\|^{2}$ in terms of the maximum eigenvalue of a random positive-semidefinite matrix:

We have used the fact that $\bm{RR}^{*}=\bm{R}$, and the notation $(\bm{PB})_{:k}$ refers to the $k$th column of the matrix $\bm{PB}$. Observe that the random positive-semidefinite matrix on the right-hand side has dimension $d$. Invoking the matrix Chernoff inequality (8), conditional on the choice of $\bm{P}$, we obtain

The required calculation is analogous with the one in the Section 1, so we omit the details. To reach a deterministic bound, we still have two more expectations to control.

Next, we examine the term in (3) that involves the maximum eigenvalue:

The first identity holds because $\lambda_{\max}(\bm{CC}^{*})=\lambda_{\max}(\bm{C}^{*}\bm{C})$ for any matrix $\bm{C}$, and $\bm{PP}^{*}=\bm{P}$. Observe that the random positive-semidefinite matrix on the right-hand side has dimension $n$, and apply the matrix Chernoff inequality (8) again to reach

Recall that $\lambda_{\max}(\bm{B}^{*}\bm{B})=\left\|{\bm{B}}\right\|^{2}$ to simplify this expression slightly.

Last, we develop a bound on the maximum column norm in (3). This result also follows from the matrix Chernoff inequality, but we need to do a little work to see why. There are more direct proofs, but this approach is closer in spirit to the rest of our proof.

We are going to treat the maximum column norm as the maximum eigenvalue of a sum of independent, random diagonal matrices. Observe that

To activate the matrix Chernoff bound, we need to compute the two parameters that appear in (8). First, the uniform upper bound $L$ satisfies

Second, to compute $\mu_{\max}$, note that

Take the maximum eigenvalue of this expression to reach

Therefore, the matrix Chernoff inequality implies

On average, the maximum squared column norm of a random submatrix $\bm{PB}$ with approximately $p$ nonzero rows gets its share $p/d$ of the maximum squared column norm of $\bm{B}$, plus a fluctuation term that depends on the magnitude of the largest entry of $\bm{B}$ and the logarithm of the number $n$ of columns.

Combine the three bounds (3), (4), and (5) to reach the result (2). We have simplified numerical constants to make the expression more compact.

Application: When is an Erdős–Rényi Graph Connected?

Random graph theory concerns probabilistic models for the interactions between pairs of objects. One basic question about a random graph is to ask whether there is a path connecting every pair of vertices or whether there are vertices segregated in different parts of the graph. It is possible to address this problem by studying the eigenvalues of random matrices, a challenge that we take up in this section.

Recall that an undirected graph is a pair $G=(V,E)$. The elements of the set $V$ are called vertices. The set $E$ is a collection of unordered pairs $\{u,v\}$ of distinct vertices, called edges. We say that the graph has an edge between vertices $u$ and $v$ in $V$ if the pair $\{u,v\}$ appears in $E$. For simplicity, we assume that the vertex set $V=\{1,\dots,n\}$. The degree $\textrm{deg}(k)$ of the vertex $k$ is the number of edges in $E$ that include the vertex $k$.

There are several natural matrices associated with an undirected graph. The adjacency matrix of the graph $G$ is an $n\times n$ symmetric matrix $\bm{A}$ whose entries indicate which edges are present:

We have assumed that edges connect distinct vertices, so the diagonal entries of the matrix $\bm{A}$ equal zero. Next, define a diagonal matrix $\bm{D}=\operatorname{diag}(\textrm{deg}(1),\dots,\textrm{deg}(n))$ whose entries list the degrees of the vertices. The Laplacian $\bm{\Delta}$ and normalized Laplacian $\bm{H}$ of the graph are the matrices

These matrices and their spectral properties play a dominant role in modern graph theory. For example, the graph $G$ is connected if and only if the second-smallest eigenvalue of $\bm{\Delta}$ is strictly positive. The second smallest eigenvalue of $\bm{H}$ controls the rate at which a random walk on the graph $G$ converges to the stationary distribution (under appropriate assumptions). See the book [GR01] for more information about these connections.

The Model of Erdős & Rényi

The simplest possible example of a random graph is the independent model $G(n,p)$ of Erdős and Rényi [ER60]. The number $n$ is the number of vertices in the graph, and $p\in(0,1)$ is the probability that two vertices are connected. More precisely, here is how to construct a random graph in $G(n,p)$. Between each pair of distinct vertices, we place an edge independently at random with probability $p$. In other words, the adjacency matrix takes the form

The family $\{\xi_{jk}:1\leq j<k\leq n\}$ consists of mutually independent $\textsc{bernoulli}(p)$ random variables. Figure 1 shows one realization of the adjacency matrix of an Erdős–Rényi graph.

Let us explain how to represent the adjacency matrix and Laplacian matrix of an Erdős–Rényi graph as a sum of independent random matrices. The adjacency matrix $\bm{A}$ of a random graph in $G(n,p)$ can be written as

This expression is a straightforward translation of the definition (1) into matrix form. Similarly, the Laplacian matrix $\bm{\Delta}$ of the random graph can be expressed as

To verify the formula (3), observe that the presence of an edge between the vertices $j$ and $k$ increases the degree of $j$ and $k$ by one. Therefore, when $\xi_{jk}=1$, we augment the $(j,j)$ and $(k,k)$ entries of $\bm{\Delta}$ to reflect the change in degree, and we mark the $(j,k)$ and $(k,j)$ entries with $-1$ to reflect the presence of the edge between $j$ and $k$.

Connectivity of an Erdős–Rényi Graph

We will obtain a near-optimal bound for the range of parameters where an Erdős–Rényi graph $G(n,p)$ is likely to be connected. We can accomplish this goal by showing that the second smallest eigenvalue of the $n\times n$ random Laplacian matrix $\bm{\Delta}=\bm{D}-\bm{A}$ is strictly positive. We will solve the problem by using the matrix Chernoff inequality to study the second-smallest eigenvalue of the random Laplacian $\bm{\Delta}$.

We need to form a random matrix $\bm{Y}$ that consists of independent positive-semidefinite terms and whose minimum eigenvalue coincides with the second-smallest eigenvalue of $\bm{\Delta}$. Our approach is to compress the matrix $\bm{Y}$ to the orthogonal complement of the vector $\mathbf{e}$ of ones. To that end, we introduce an $(n-1)\times n$ partial isometry $\bm{R}$ that satisfies

Now, consider the $(n-1)\times(n-1)$ random matrix

Recall that $\{\xi_{jk}\}$ is an independent family of $\textsc{bernoulli}(p)$ random variables, so the summands are mutually independent. The Conjugation Rule (12) ensures that each summand remains positive semidefinite. Furthermore, the Courant–Fischer theorem implies that the minimum eigenvalue of $\bm{Y}$ coincides with the second-smallest eigenvalue of $\bm{\Delta}$ because the smallest eigenvalue of $\bm{\Delta}$ has eigenvector $\mathbf{e}$.

To apply the matrix Chernoff inequality, we show that $L=2$ is an upper bound for the eigenvalues of each summand in (4). We have

The first bound follows from the submultiplicativity of the spectral norm. To obtain the second bound, note that $\xi_{jk}$ takes 0–1 values. The matrix $\bm{R}$ is a partial isometry so its norm equals one. Finally, a direct calculation shows that $\bm{T}=\mathbf{E}_{jj}+\mathbf{E}_{kk}-\mathbf{E}_{jk}-\mathbf{E}_{kj}$ satisfies the polynomial $\bm{T}^{2}=2\bm{T}$, so each eigenvalue of $\bm{T}$ must equal zero or two.

Next, we compute the expectation of the matrix $\bm{Y}$.

The first identity follows when we apply linearity of expectation to (5) and then use linearity of matrix multiplication to draw the sum inside the conjugation by $\bm{R}$. The term $(n-1)\,\mathbf{I}_{n}$ emerges when we sum the diagonal matrices. The term $\mathbf{ee}^{*}-\mathbf{I}_{n}$ comes from the off-diagonal matrix units, once we note that the matrix $\mathbf{ee}^{*}$ has one in each component. The last identity holds because of the properties of $\bm{R}$ displayed in (4). We conclude that

To arrive at a probability inequality for the second-smallest eigenvalue $\lambda_{2}^{\uparrow}(\bm{\Delta})$ of the matrix $\bm{\Delta}$, we apply the tail bound (5) to the matrix $\bm{Y}$. We obtain, for $t\in(0,1)$,

for an Erdős–Rényi graph $G(n,p)$ to be connected with high probability as $n\to\infty$. This bound is quite close to the optimal result, which lacks the factor two on the right-hand side. It is possible to make this reasoning more precise, but it does not seem worth the fuss.

Proof of the Matrix Chernoff Inequalities

The first step toward the matrix Chernoff inequalities is to develop an appropriate semidefinite bound for the mgf and cgf of a random positive-semidefinite matrix. The method for establishing this result mimics the proof in the scalar case: we simply bound the exponential with a linear function.

Suppose that $\bm{X}$ is a random matrix that satisfies $0\leq\lambda_{\min}(\bm{X})$ and $\lambda_{\max}(\bm{X})\leq L$. Then

By assumption, each eigenvalue of $\bm{X}$ lies in the interval $[0,L]$. Thus, the Transfer Rule (14) implies that

Expectation respects the semidefinite order, so

To obtain the semidefinite bound for the cgf, we simply take the logarithm of the semidefinite bound for the mgf. This operation preserves the semidefinite order because of the property (18) that the logarithm is operator monotone. ∎

We break the proof of the matrix inequality into two pieces. First, we establish the bounds on the maximum eigenvalue, which are slightly easier. Afterward, we develop the bounds on the minimum eigenvalue.

Consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random Hermitian matrices with common dimension $d$. Assume that

Next, we turn to the upper bound (6) for the upper tail of the maximum eigenvalue. Substitute the cgf bounds (1) into the master inequality (3) to reach

The steps here are identical with the previous argument. To complete the proof, make the change of variables $t\mapsto(1+\varepsilon)\mu_{\max}$. Then the infimum is achieved at $\theta=L^{-1}\log(1+\varepsilon)$, which leads to the tail bound (6). ∎

The lower bounds follow from a related argument that is slightly more delicate.

Once again, consider a finite sequence $\{\bm{X}_{k}\}$ of independent, random Hermitian matrices with dimension $d$. Assume that

Note that $g(\theta)<0$ for $\theta<0$, which alters a number of the steps in the argument.

Most of the steps are the same as in the proof of the upper bound (4), so we focus on the differences. Since the factor $\theta^{-1}$ in the first and second lines is negative, upper bounds on the trace reduce the value of the expression. We move to the fourth line by invoking the property $\lambda_{\max}(\alpha\bm{A})=\alpha\lambda_{\min}(\bm{A})$ for $\alpha<0$, which follows from (4) and (5). This piece of algebra depends on the fact that $g(\theta)<0$ when $\theta<0$. To obtain the result (3), we change variables: $\theta\mapsto-\theta/L$.

Finally, we establish the bound (5) for the lower tail of the minimum eigenvalue. Introduce the cgf bounds (2) into the master inequality (4) to reach

The justifications here match those in with the previous argument. Finally, we make the change of variables $t\mapsto(1-\varepsilon)\mu_{\min}$. The infimum is attained at $\theta=L^{-1}\log(1-\varepsilon)$, which yields the tail bound (5). ∎

Notes

As usual, we continue with an overview of background references and related work.

Scalar Chernoff inequalities date to the paper [Che52, Thm. 1] by Herman Chernoff. The original result provides probability bounds for the number of successes in a sequence of independent but non-identical Bernoulli trials. Chernoff’s proof combines the scalar Laplace transform method with refined bounds on the mgf of a Bernoulli random variable. It is very common to encounter simplified versions of Chernoff’s result, such as [Lug09, Exer. 8] or [MR95, §4.1].

In their paper [AW02], Ahlswede & Winter developed a matrix version of the Chernoff inequality. The matrix mgf bound, Lemma 4.1, essentially appears in their work. Ahlswede & Winter focus on the case of independent and identically distributed random matrices, in which case their results are roughly equivalent with Theorem 1.1. For the general case, their approach leads to matrix expectation statistics of the form

It is clear that their $\mu_{\min}^{\textrm{AW}}$ may be substantially smaller than the quantity $\mu_{\min}$ we defined in Theorem 1.1. Similarly, their $\mu_{\max}^{\textrm{AW}}$ may be substantially larger than the quantity $\mu_{\max}$ that drives the upper Chernoff bounds.

The tail bounds from Theorem 1.1 are drawn from [Tro11c, §5], but the expectation bounds we present are new. The technical report [GT14] extends the matrix Chernoff inequality to provide upper and lower tail bounds for all eigenvalues of a sum of random, positive-semidefinite matrices. Chapter 7 contains a slight improvement of the bounds for the maximum eigenvalue in Theorem 1.1.

Let us mention a few other results that are related to the matrix Chernoff inequality. First, Theorem 1.1 has a lovely information-theoretic formulation where the tail bounds are stated in terms of an information divergence. To establish this result, we must restructure the proof and eliminate some of the approximations. See [AW02, Thm. 19] or [Tro11c, Thm. 5.1].

Second, the problem of bounding the minimum eigenvalue of a sum of random, positive-semidefinite matrices has a special character. The reason, roughly, is that a sum of independent, nonnegative random variables cannot easily take the value zero. A closely related phenomenon holds in the matrix setting, and it is possible to develop estimates that exploit this observation. See [Oli13, Thm. 3.1] and [KM13, Thm. 1.3] for two wildly different approaches.

The Matrix Rosenthal Inequality

The matrix Rosenthal inequality (9) is one of the earliest matrix concentration bounds. In his paper [Rud99], Rudelson used the noncommutative Khintchine inequality (1) to establish a specialization of (9) to rank-one summands. A refinement appears in [RV07], and explicit constants were first derived in [Tro08c]. We believe that the paper [CGT12a] contains the first complete statement of the moment bound (9) for general positive-semidefinite summands; see also the work [MZ11]. The constants in [CGT12a, Thm. A.1], and hence in (9), can be improved slightly by using the sharp version of the noncommutative Khintchine inequality from [Buc01, Buc05]. Let us stress that all of these results follow from easy variations of Rudelson’s argument.

The work [MJC+14, Cor. 7.4] provides a self-contained and completely elementary proof of a matrix Rosenthal inequality that is closely related to (9). This result depends on different principles from the works mentioned in the last paragraph.

Random Submatrices

The problem of studying a random submatrix drawn from a fixed matrix has a long history. An early example is the paving problem from operator theory, which asks for a maximal well-conditioned set of columns (or a well-conditioned submatrix) inside a fixed matrix. Random selection provides a natural way to approach this question. The papers of Bourgain & Tzafriri [BT87, BT91] and Kashin & Tzafriri [KT94] study random paving using sophisticated tools from functional analysis. See the paper [NT14] for a summary of research on randomized methods for constructing pavings. Very recently, Adam Marcus, Dan Spielman, & Nikhil Srivastava [MSS14] have solved the paving problem completely.

The article [Tro11d] contains the observation that the matrix Chernoff inequality is an ideal tool for studying random submatrices. It applies this technique to study a random matrix that arises in numerical linear algebra [HMT11], and it achieves an optimal estimate for the minimum singular value of the random matrix that arises in this setting. Our analysis of a random column submatrix is based on this work. The analysis of a random row and column submatrix is new. The paper [CD12], by Chrétien and Darses, uses matrix Chernoff bounds in a more sophisticated way to develop tail bounds for the norm of a random row and column submatrix.

Random Graphs

The analysis of random graphs and random hypergraphs appeared as one of the earliest applications of matrix concentration inequalities [AW02]. Christofides and Markström developed a matrix Hoeffding inequality to aid in this purpose [CM08]. Later, Oliveira wrote two papers [Oli10a, Oli11] on random graph theory based on matrix concentration. We recommend these works for further information.

To analyze the random graph Laplacian, we compressed the Laplacian to a subspace so that the minimum eigenvalue of the compression coincides with the second-smallest eigenvalue of the original Laplacian. This device can be extended to obtain tail bounds for all the eigenvalues of a sum of independent random matrices. See the technical report [GT14] for a development of this idea.

Chapter 6 A Sum of Bounded Random Matrices

In this chapter, we describe matrix concentration inequalities that generalize the classical Bernstein bound. The matrix Bernstein inequalities concern a random matrix formed as a sum of independent, random matrices that are bounded in spectral norm. The results allow us to study how much this type of random matrix deviates from its mean value in the spectral norm.

Formally, we consider an finite sequence $\{\bm{S}_{k}\}$ of random matrices of the same dimension. Assume that the matrices satisfy the conditions

Form the sum $\bm{Z}=\sum_{k}\bm{S}_{k}$. The matrix Bernstein inequality controls the expectation and tail behavior of $\left\|{\bm{Z}}\right\|$ in terms of the matrix variance statistic $v(\bm{Z})$ and the uniform bound $L$.

The matrix Bernstein inequality is a powerful tool with a huge number of applications. In these pages, we can only give a coarse indication of how researchers have used this result, so we have chosen to focus on problems that use random sampling to approximate a specified matrix. This model applies to the sample covariance matrix in the introduction. In this chapter, we outline several additional examples. First, we consider the technique of randomized sparsification, in which we replace a dense matrix with a sparse proxy that has similar spectral behavior. Second, we explain how to develop a randomized algorithm for approximate matrix multiplication, and we establish an error bound for this method. Third, we develop an analysis of random features, a method for approximating kernel matrices that has become popular in contemporary machine learning.

As these examples suggest, the matrix Bernstein inequality is very effective for studying randomized approximations of a given matrix. Nevertheless, when the matrix Chernoff inequality, Theorem 1.1, happens to apply to a problem, it often delivers better results.

Section 1 describes the matrix Bernstein inequality. Section 2 explains how to use the Bernstein inequality to study randomized methods for matrix approximation. In §§3, 4, and 5, we apply the latter result to three matrix approximation problems. We conclude with the proof of the matrix Bernstein inequality in §6.

A Sum of Bounded Random Matrices

In the scalar setting, the label “Bernstein inequality” applies to a very large number of concentration results. Most of these bounds have extensions to matrices. For simplicity, we focus on the most famous of the scalar results, a tail bound for the sum $Z$ of independent, zero-mean random variables that are subject to a uniform bound. In this case, the Bernstein inequality shows that $Z$ concentrates around zero. The tails of $Z$ make a transition from subgaussian decay at moderate deviations to subexponential decay at large deviations. See [BLM13, §2.7] for more information about Bernstein’s inequality.

In analogy, the simplest matrix Bernstein inequality concerns a sum of independent, zero-mean random matrices whose norms are bounded above. The theorem demonstrates that the norm of the sum acts much like the scalar random variable $Z$ that we discussed in the last paragraph.

Consider a finite sequence $\{\bm{S}_{k}\}$ of independent, random matrices with common dimension $d_{1}\times d_{2}$. Assume that

Let $v(\bm{Z})$ be the matrix variance statistic of the sum:

Let us spend a few moments to discuss the matrix Bernstein inequality, Theorem 1.1, its consequences, and some of the improvements that are available.

First, observe that the matrix variance statistic $v(\bm{Z})$ appearing in (1) coincides with the general definition (8) because $\bm{Z}$ has zero mean. To reach (2), we have used the additivity law (11) for an independent sum to express the matrix variance statistic in terms of the summands. Observe that, when the summands $\bm{S}_{k}$ are Hermitian, the two terms in the maximum coincide.

Next, let us explain how to interpret the tail bound (4). The main difference between this result and the scalar Bernstein bound is the appearance of the dimensional factor $d_{1}+d_{2}$, which reduces the range of $t$ where the inequality is informative. To get a better idea of what this result means, it is helpful to make a further estimate:

In other words, for moderate values of $t$, the tail probability decays as fast as the tail of a Gaussian random variable whose variance is comparable with $v(\bm{Z})$. For larger values of $t$, the tail probability decays at least as fast as that of an exponential random variable whose mean is comparable with $L$. As usual, we insert a warning that the tail behavior reported by the matrix Bernstein inequality can overestimate the actual tail behavior.

Last, it is helpful to remember that the matrix Bernstein inequality extends to a sum of uncentered random matrices. In this case, the result describes the spectral-norm deviation of the random sum from its mean value. For reference, we include the statement here.

Consider a finite sequence $\{\bm{S}_{k}\}$ of independent random matrices with common dimension $d_{1}\times d_{2}$. Assume that each matrix has uniformly bounded deviation from its mean:

and let $v(\bm{Z})$ denote the matrix variance statistic of the sum:

This result follows as an immediate corollary of Theorem 1.1.

The bounds in Theorem 1.1 are stated in terms of the ambient dimensions $d_{1}$ and $d_{2}$ of the random matrix $\bm{Z}$. The dependence on the ambient dimension is not completely natural. For example, consider embedding the random matrix $\bm{Z}$ into the top corner of a much larger matrix which is zero everywhere else. It turns out that we can achieve results that reflect only the “intrinsic dimension” of $\bm{Z}$. We turn to this analysis in Chapter 7.

In addition, there are many circumstances where the uniform upper bound $L$ that appears in (3) does not accurately reflect the tail behavior of the random matrix. For instance, the summands themselves may have very heavy tails. In such emergencies, the following expectation bound [CGT12a, Thm. A.1] can be a lifesaver.

This result is a matrix formulation of the Rosenthal–Pinelis inequality [Pin94, Thm. 4.1].

Finally, let us reiterate that there are other types of matrix Bernstein inequalities. For example, we can sharpen the tail bound (4) to obtain a matrix Bennett inequality. We can also relax the boundedness assumption to a weaker hypothesis on the growth of the moments of each summand $\bm{S}_{k}$. In the Hermitian setting, the result can also discriminate the behavior of the upper and lower tails, which is a consequence of Theorem 6.1 below. See the notes at the end of this chapter and the annotated bibliography for more information.

Optimality of the Matrix Bernstein Inequality

To use the matrix Bernstein inequality, Theorem 1.1, and its relatives with intelligence, one must appreciate their strengths and weaknesses. We will focus on the matrix Rosenthal–Pinelis inequality (6). Nevertheless, similar insights are relevant to the estimate (3).

Let us present lower bounds to demonstrate that the matrix Rosenthal–Pinelis inequality (6) requires both terms that appear. First, the quantity $v(\bm{Z})$ cannot be omitted because Jensen’s inequality implies that

Under a natural hypothesis, the second term on the right-hand side of (6) also is essential. Suppose that each summand $\bm{S}_{k}$ is a symmetric random variable; that is, $\bm{S}_{k}$ and $-\bm{S}_{k}$ have the same distribution. In this case, an involved argument [LT91, Prop. 6.10] leads to the bound

There are examples where the right-hand side of (7) is comparable with the uniform upper bound $L$ on the summands, but this is not always so.

In summary, when the summands $\bm{S}_{k}$ are symmetric, we have matching estimates

We see that the bound (6) must include some version of each term that appears, but the logarithms are not always necessary.

First, let us show that the variance term in (6) must contain a logarithm. For each natural number $n$, consider the $d\times d$ random matrix $\bm{Z}$ of the form

where $\{\varrho_{ik}\}$ is an independent family of Rademacher random variables. An easy application of the bound (6) implies that

Using the central limit theorem and the Skorokhod representation, we can construct an independent family $\{\gamma_{k}\}$ of standard normal random variables for which

Therefore, we cannot remove the logarithm from the variance term in (6).

Next, let us justify the logarithm on the norm of the summands in (6). For each natural number $n$, consider a $d\times d$ random matrix $\bm{Z}$ of the form

where $\big\{\delta_{ik}^{(n)}\big\}$ is an independent family of $\textsc{bernoulli}(n^{-1})$ random variables. The matrix Rosenthal–Pinelis inequality (6) ensures that

Using the Poisson limit of a binomial random variable and the Skorohod representation, we can construct an independent family $\{Q_{k}\}$ of $\textsc{poisson}(1)$ random variables for which

In short, the bound we derived from (6) requires the logarithm on the second term, but it is suboptimal by a $\log\log$ factor. The upper matrix Chernoff inequality (6) correctly predicts the appearance of the iterated logarithm in this example, as does the matrix Bennett inequality.

The last two examples rely heavily on the commutativity of the summands as well as the infinite divisibility of the normal and Poisson distributions. As a consequence, it may appear that the logarithms only appear in very special contexts. In fact, many (but not all!) examples that arise in practice do require the logarithms that appear in the matrix Bernstein inequality. It is a subject of ongoing research to obtain a simple criterion for deciding when the logarithms belong.

Example: Matrix Approximation by Random Sampling

In applied mathematics, we often need to approximate a complicated target object by a more structured object. In some situations, we can solve this problem using a beautiful probabilistic approach called empirical approximation. The basic idea is to construct a “simple” random object whose expectation equals the target. We obtain the approximation by averaging several independent copies of the simple random object. As the number of terms in this average increases, the approximation becomes more complex, but it represents the target more faithfully. The challenge is to quantify this tradeoff.

In particular, we often encounter problems where we need to approximate a matrix by a more structured matrix. For example, we may wish to find a sparse matrix that is close to a given matrix, or we may need to construct a low-rank matrix that is close to a given matrix. Empirical approximation provides a mechanism for obtaining these approximations. The matrix Bernstein inequality offers a natural tool for assessing the quality of the randomized approximation.

This section develops a general framework for empirical approximation of matrices. Subsequent sections explain how this technique applies to specific examples from the fields of randomized linear algebra and machine learning.

Let $\bm{B}$ be a target matrix that we hope to approximate by a more structured matrix. To that end, let us represent the target as a sum of “simple” matrices:

The idea is to identify summands with desirable properties that we want our approximation to inherit. The examples in this chapter depend on decompositions of the form (1).

Along with the decomposition (1), we need a set of sampling probabilities:

We want to ascribe larger probabilities to “more important” summands. Quantifying what “important” means is the most difficult aspect of randomized matrix approximation. Choosing the right sampling distribution for a specific problem requires insight and ingenuity.

Given the data (1) and (2), we may construct a “simple” random matrix $\bm{R}$ by sampling:

Our goal is to quantify the approximation error as a function of the complexity $n$ of the approximation:

Error Estimate for Matrix Sampling Estimators

We can obtain an error estimate for the approximation scheme described in Section 1 as an immediate corollary of the matrix Bernstein inequality, Theorem 1.1.

Let $\bm{B}$ be a fixed $d_{1}\times d_{2}$ matrix. Construct a $d_{1}\times d_{2}$ random matrix $\bm{R}$ that satisfies

Since $\bm{R}$ is an unbiased estimator of the target matrix $\bm{B}$, we can write

Now, each of the summands is subject to an upper bound:

The first relation is the triangle inequality; the second is Jensen’s inequality. The last estimate follows from our assumption that $\left\|{\bm{R}}\right\|\leq L$.

To control the matrix variance statistic $v(\bm{Z})$, first note that

The first identity follows from the expression (2) for the matrix variance statistic, and the second holds because the summands $\bm{S}_{k}$ are identically distributed. We may calculate that

The last line follows from the definition (4) of $m_{2}(\bm{R})$.

We are prepared to apply the matrix Bernstein inequality, Theorem 1.1, to the random matrix $\bm{Z}=\sum_{k}\bm{S}_{k}$. This operation results in the statement of the corollary. ∎

Discussion

One of the most common applications of the matrix Bernstein inequality is to analyze empirical matrix approximations. As a consequence, Corollary 2.1 is one of the most useful forms of the matrix Bernstein inequality. Let us discuss some of the important aspects of this result.

First, let us examine how many samples $n$ suffice to bring the approximation error bound in Corollary 2.1 below a specified positive tolerance $\varepsilon$. Examining inequality (5), we find that

Roughly, the number $n$ of samples should be on the scale of the per-sample second moment $m_{2}(\bm{R})$ and the uniform upper bound $L$.

The bound (7) also reveals an unfortunate aspect of empirical matrix approximation. To make the tolerance $\varepsilon$ small, the number $n$ of samples must increase proportional with $\varepsilon^{-2}$. In other words, it takes many samples to achieve a highly accurate approximation. We cannot avoid this phenomenon, which ultimately is a consequence of the central limit theorem.

On a more positive note, it is quite valuable that the error bounds (5) and (6) involve the spectral norm. This type of estimate simultaneously controls the error in every linear function of the approximation:

The Schatten $1$-norm $\left\|{\cdot}\right\|_{S_{1}}$ is defined in (29). These bounds also control the error in each singular value $\sigma_{j}(\bar{\bm{R}}_{n})$ of the approximation:

When there is a gap between two singular values of $\bm{B}$, we can also obtain bounds for the discrepancy between the associated singular vectors of $\bar{\bm{R}}_{n}$ and $\bm{B}$ using perturbation theory.

To construct a good sampling estimator $\bm{R}$, we ought to control both $m_{2}(\bm{R})$ and $L$. In practice, this demands considerable creativity. This observation hints at the possibility of achieving a bias–variance tradeoff when approximating $\bm{B}$. To do so, we can drop all of the “unimportant” terms in the representation (1), i.e., those whose sampling probabilities are small. Then we construct a random approximation $\bm{R}$ only for the “important” terms that remain. Properly executed, this process may decrease both the per-sample second moment $m_{2}(\bm{R})$ and the upper bound $L$. The idea is analogous with shrinkage in statistical estimation.

Corollary 2.1 extends beyond the sampling model based on the finite expansion (1). Indeed, we can consider a more general decomposition of the target matrix $\bm{B}$:

where $\mu$ is a probability measure on a sample space $\Omega$. As before, the idea is to represent the target matrix $\bm{B}$ as an average of “simple” matrices $\bm{B}(\omega)$. The main difference is that the family of simple matrices may now be infinite. In this setting, we construct the random approximation $\bm{R}$ so that

As we will discuss, this abstraction is important for applications in machine learning.

Another fundamental point about sampling estimators is that they are usually suboptimal. In other words, the matrix sampling estimator may incur an error substantially worse than the error in the best structured approximation of the target matrix.

To see why, let us consider a simple form of low-rank approximation by random sampling. The method here does not have practical value, but it highlights the reason that sampling estimators usually do not achieve ideal results. Suppose that $\bm{B}$ has singular value decomposition

Given the SVD, we can construct a random rank-one approximation $\bm{R}$ of the form

Per Corollary 2.1, the error in the associated sampling estimator $\bar{\bm{R}}_{n}$ of $\bm{B}$ satisfies

On the other hand, a best rank-$n$ approximation of $\bm{B}$ takes the form $\bm{B}_{n}=\sum_{j=1}^{n}\sigma_{j}\bm{u}_{j}\bm{v}_{j}^{*}$, and it incurs error

The second relation is Markov’s inequality, which provides an accurate estimate only when the singular values $\sigma_{1},\dots,\sigma_{n+1}$ are comparable. In that case, the sampling estimator arrives within a logarithmic factor of the optimal error. But there are many matrices whose singular values decay quickly, so that $\sigma_{n+1}\ll(n+1)^{-1}$. In the latter situation, the error in the sampling estimator is much worse than the optimal error.

We often encounter papers that develop Frobenius-norm error bounds for matrix approximations, perhaps because the analysis is more elementary. But one must recognize that Frobenius-norm error bounds are not acceptable in most cases of practical interest:

Frobenius-norm error bounds are typically vacuous.

In particular, this phenomenon occurs in data analysis whenever we try to approximate a matrix that contains white or pink noise.

To illustrate this point, let us consider the ubiquitous problem of approximating a low-rank matrix corrupted by additive white Gaussian noise:

Now, the spectral-norm error in the desired approximation satisfies

On the other hand, the Frobenius-norm error in the desired approximation satisfies

We see that the Frobenius-norm error can be quite large, even when we find the required approximation.

Here is another way to look at the same fact. Suppose we construct an approximation $\widehat{\bm{B}}$ of the matrix $\bm{B}$ from (8) whose Frobenius-norm error is comparable with the optimal error:

Application: Randomized Sparsification of a Matrix

Many tasks in data analysis involve large, dense matrices that contain a lot of redundant information. For example, an experiment that tabulates many variables about a large number of subjects typically results in a low-rank data matrix because subjects are often similar with each other. Many questions that we pose about these data matrices can be addressed by spectral computations. In particular, factor analysis involves a singular value decomposition.

When the data matrix is approximately low rank, it has fewer degrees of freedom than its ambient dimension. Therefore, we can construct a simpler approximation that still captures most of the information in the matrix. One method for finding this approximation is to replace the dense target matrix by a sparse matrix that is close in spectral-norm distance. An elegant way to identify this sparse proxy is to randomly select a small number of entries from the original matrix to retain. This is a type of empirical approximation.

Sparsification has several potential advantages. First, it is considerably less expensive to store a sparse matrix than a dense matrix. Second, many algorithms for spectral computation operate more efficiently on sparse matrices.

In this section, we examine a very recent approach to randomized sparsification due to Kundu & Drineas [KD14]. The analysis is an immediate consequence of Corollary 2.1. See the notes at the end of the chapter for history and references.

Let $\bm{B}$ be a fixed $d_{1}\times d_{2}$ complex matrix. The sparsification problem requires us to find a sparse matrix $\widehat{B}$ that has small distance from $\bm{B}$ with respect to the spectral norm. We can achieve this goal using an empirical approximation strategy.

First, let us express the target matrix as a sum of its entries:

Now, we introduce a $d_{1}\times d_{2}$ random matrix $\bm{R}$ that has exactly one nonzero entry:

We use the convention that $0/0=0$ so that we do not need to treat zero entries separately. It is immediate that

Therefore, $\bm{R}$ is an unbiased estimate of $\bm{B}$.

Although the expectation of $\bm{R}$ is correct, its variance is quite high. Indeed, $\bm{R}$ has only one nonzero entry, while $\bm{B}$ typically has many nonzero entries. To reduce the variance, we combine several independent copies of the simple estimator:

Performance of Randomized Sparsification

The randomized sparsification method is clearly a type of empirical approximation, so we can use Corollary 2.1 to perform the analysis. We will establish the following error bound.

The short proof of (2) appears below in Section 3.

Placing the error (2) on a relative scale, we see that

The stable rank $\operatorname{srank}(\bm{B})$, defined in (25), emerges naturally as a quantity of interest.

Now, suppose that the sparsity level $n$ satisfies

where the tolerance $\varepsilon\in(0,1]$. We determine that

Since the stable rank always exceeds one and we have assumed that $\varepsilon\leq 1$, this estimate implies that

We discover that it is possible to replace the matrix $\bm{B}$ by a matrix with at most $n$ nonzero entries while achieving a small relative error in the spectral norm. When $\operatorname{srank}(\bm{B})\ll\min\{d_{1},\ d_{2}\}$, we can achieve a dramatic reduction in the number of nonzero entries needed to carry the spectral information in the matrix $\bm{B}$.

Analysis of Randomized Sparsification

Let us proceed with the analysis of randomized sparsification. To apply Corollary 2.1, we need to obtain bounds for the per-sample variance $m_{2}(\bm{R})$ and the uniform upper bound $L$. The key to both calculations is to obtain appropriate lower bounds on the sampling probabilities $p_{ij}$. Indeed,

Each estimate follows by neglecting one term in (3).

First, we turn to the uniform bound on the random matrix $\bm{R}$. We have

Second, we turn to the computation of the per-sample second moment $m_{2}(\bm{R})$. We have

The semidefinite inequality holds because each matrix ${\left|{\smash{b_{ij}}}\right|}^{2}\mathbf{E}_{ii}$ is positive semidefinite and because of the second bound in (3). Similarly,

This is the required estimate for the per-sample second moment.

Application: Randomized Matrix Multiplication

Numerical linear algebra (NLA) is a well-established and important part of computer science. Some of the basic problems in this area include multiplying matrices, solving linear systems, computing eigenvalues and eigenvectors, and solving linear least-squares problems. Historically, the NLA community has focused on developing highly accurate deterministic methods that require as few floating-point operations as possible. Unfortunately, contemporary applications can strain standard NLA methods because problems have continued to become larger. Furthermore, on modern computer architectures, computational costs depend heavily on communication and other resources that the standard algorithms do not manage very well.

In response to these challenges, researchers have started to develop randomized algorithms for core problems in NLA. In contrast to the classical algorithms, these new methods make random choices during execution to achieve computational efficiencies. These randomized algorithms can also be useful for large problems or for modern computer architectures. On the other hand, randomized methods can fail with some probability, and in some cases they are less accurate than their classical competitors.

Matrix concentration inequalities are one of the key tools used to design and analyze randomized algorithms for NLA problems. In this section, we will describe a randomized method for matrix multiplication developed by Magen & Zouzias [MZ11, Zou13]. We will analyze this algorithm using Corollary 2.1. Turn to the notes at the end of the chapter for more information about the history.

One of the basic tasks in numerical linear algebra is to multiply two matrices with compatible dimensions. Suppose that $\bm{B}$ is a $d_{1}\times N$ complex matrix and that $\bm{C}$ is an $N\times d_{2}$ complex matrix, and we wish to compute the product $\bm{BC}$. The straightforward algorithm forms the product entry by entry:

This approach takes $\operatorname{O}(N\cdot d_{1}d_{2})$ arithmetic operations. There are algorithms, such as Strassen’s divide-and-conquer method, that can reduce the cost, but these approaches are not considered practical for most applications.

Suppose that the inner dimension $N$ is substantially larger than the outer dimensions $d_{1}$ and $d_{2}$. In this setting, both matrices $\bm{B}$ and $\bm{C}$ are rank-deficient, so the columns of $\bm{B}$ contain a lot of linear dependencies, as do the rows of $\bm{C}$. As a consequence, a random sample of columns from $\bm{B}$ (or rows from $\bm{C}$) can be used as a proxy for the full matrix. Formally, the key to this approach is to view the matrix product as a sum of outer products:

As usual, $\bm{b}_{:j}$ denotes the $j$th column of $\bm{B}$, while $\bm{c}_{j:}$ denotes the $j$th row of $\bm{C}$. We can approximate this sum using the empirical method.