Tail bounds for all eigenvalues of a sum of random matrices

Introduction

The field of nonasymptotic random matrix theory has traditionally focused on the problem of bounding the extreme eigenvalues of a random matrix. In some circumstances, however, we may also be interested in studying the behavior of the interior eigenvalues. In this case, classical tools do not readily apply. Indeed, the interior eigenvalues are determined by the min-max of a random process, which is very challenging to control.

This paper demonstrates that it is possible to combine the matrix Laplace transform method detailed in [Tro11c] with the Courant–Fischer characterization of eigenvalues to obtain nontrivial bounds on the interior eigenvalues of a sum of random self-adjoint matrices. This approach expands the scope of the matrix probability inequalities from [Tro11c] so that they provide interesting information about the bulk spectrum.

As one application of our approach, we investigate estimates for the covariance matrix of a centered stationary random process. We show that the eigenvalues of the sample covariance matrix provide relative-error approximations to the eigenvalues of the covariance matrix. We focus on Gaussian processes, but our arguments can be extended to other distributions. The following theorem distills the results in section 7.

We believe that this paper contains the first general-purpose tools for studying the full spectrum of a finite-dimensional random matrix. The literature on random matrix theory (RMT) contains some complementary results, but they do not seem to apply with the same generality. Methods from RMT fall into two rough categories: asymptotic methods and nonasymptotic methods. We discuss the relevant results from each in turn.

The modern asymptotic theory began in the 1950s when physicists observed that, on certain scales, the behavior of a quantum system is described by the spectrum of a random matrix [Meh04]. They further observed the phenomenon of universality: as the dimension increases, the spectral statistics become independent of the distribution of the random matrix; instead, they are determined by the symmetries of the distribution [Dei07]. Since these initial observations, physicists, statisticians, engineers, and mathematicians have found manifold applications of the asymptotic theory in high-dimensional statistics [Joh01, Joh07, El 08], physics [GMGW98, Meh04], wireless communication [TV04, ST06], and pure mathematics [RS96, BK99], to mention only a few areas.

Asymptotic random matrix theory has developed primarily through the examination of specific classes of random matrices. We mention two well-studied classes. Sample covariance matrices take the form $n^{-1}\bm{B}_{n}\bm{B}_{n}^{*}$, where the columns of $\bm{B}_{n}$ comprise $n$ independent observations. Wigner matrices are Hermitian matrices whose superdiagonal entries are independent, zero-mean, and have unit variance and whose diagonal entries are i.i.d., real, and have finite variance.

The fundamental object of study in asymptotic random matrix theory is the empirical spectral distribution function (ESD). Given a random Hermitian matrix $\bm{A}$ of order $n$, its ESD

is a random distribution function which encodes the statistics of the spectrum of $\bm{A}.$ Wigner’s theorem [Wig55], the seminal result of the asymptotic theory, establishes that if $\{\bm{A}_{n}\}$ is a sequence of independent, symmetric $n\times n$ matrices with i.i.d. $\mathcal{N}(0,1)$ entries on and above the diagonal, then the expected ESD of $n^{-1/2}\bm{A}_{n}$ converges weakly in probability, as $n$ approaches infinity, to the semicircular law given by

Thus, at least in the limiting sense, the spectra of these random matrices are well characterized. Development of the classical asymptotic theory has been driven by the natural question raised by Wigner’s result: to what extent is the semicircular law, and more generally, the existence of a limiting spectral distribution (LSD) universal?

The literature on the existence and universality of LSDs is massive; we mention only the highlights. It is now known that the semicircular law is universal for Wigner matrices. Suppose that $\{\bm{A}_{n}\}$ is a sequence of independent $n\times n$ Wigner matrices. Grenander established that if all the moments are finite, then the ESD of $n^{-1/2}\bm{A}_{n}$ converges weakly to the semicircular law in probability [Gre63]. Arnold showed that, assuming a finite fourth moment, the ESD almost surely converges weakly to the semicircular law [Arn71]. Around the same time, Marc̆enko and Pastur determined the form of the limiting spectral distribution of sample covariance matrices [MP67].

More recently, Tao and Vu confirmed the long-conjectured circular law hypothesis. Let $\{\bm{C}_{n}\}$ be a sequence of independent $n\times n$ matrices whose entries are i.i.d. and have unit variance. Then the ESD of $n^{-1/2}\bm{C}_{n}$ converges weakly to the uniform measure on the unit disk, both in probability and almost surely [TV10b].

Although the convergence rate of the ESD has considerable practical interest, it was not until 1993 that theoretical results became available when Bai showed that for Wigner matrices [Bai93a] and sample covariance matrices [Bai93b] the expected ESDs of $n^{-1/2}\bm{A}_{n}$ and $n^{-1}\bm{B}_{n}\bm{B}_{n}^{*},$ respectively, both converge pointwise at a rate of O($n^{-1/4}$). Later, Bai and coauthors established the pointwise convergence in probability of the ESD of the normalized Wigner matrix $n^{-1/2}\bm{A}_{n}$ [BMT97] and greatly improved the convergence rates [BMT99, BMT02, BMY03]. The strongest result to date is due to Bai et al., who have shown that, if the entries of the Wigner matrix possess finite sixth moments, then pointwise convergence in probability of the ESD of $n^{-1/2}\bm{A}_{n}$ occurs at the rate of O($n^{-1/2}$) [BHPZ11].

Classically, individual eigenvalues have been studied through the limiting behavior of the extremal eigenvalues and the asymptotic joint distribution of several eigenvalues. Much is known about the limiting distribution of the largest eigenvalues of Wigner and covariance matrices. Geman showed that if the columns of $\bm{B}_{n}$ are drawn from a sufficiently regular distribution, then the largest eigenvalue of the sample covariance matrix $n^{-1}\bm{B}_{n}\bm{B}_{n}^{*}$ converges almost surely to a limit [Gem80]. Bai, Yin, and coauthors showed that the existence of a fourth moment is both necessary and sufficient for the existence of such a limit [YBK88, BSY88]. They also identified necessary and sufficient conditions for the existence of limits for the smallest and largest eigenvalues of a normalized Wigner matrix $n^{-1/2}\bm{A}_{n}$ [BY88b]. El Karoui has recently described the limiting behavior of the leading eigenvalues of a large class of sample covariance matrices [El 07].

Less is known about the rate of convergence of the eigenvalues, but some results are available. Write the eigenvalues of a self-adjoint matrix $\bm{A}$ in nonincreasing order $\lambda_{1}\geq\ldots\geq\lambda_{n}.$ For $1\leq j\leq n,$ the classical location $\gamma_{j}$ of the $j$th eigenvalue of the normalized Wigner matrix $n^{-1/2}\bm{A}_{n}$ is defined via the relation

where $\rho_{sc}$ is the density associated with the semicircular law. Intuitively, the facts that $F^{\frac{1}{\sqrt{n}}\bm{A}_{n}}\rightarrow F^{sc}$ and $F^{\frac{1}{\sqrt{n}}\bm{A}_{n}}(\lambda_{j})=j/n$ suggest that $\frac{1}{\sqrt{n}}\lambda_{j}\rightarrow\gamma_{j}.$ Indeed, it follows from [BY88a, BY88b] that

asymptotically almost surely. Under the assumption that the entries exhibit uniform subgaussian decay, Erdös, Yau, and Yin have strengthened this result by showing that, up to log factors, the eigenvalues of $n^{-1/2}\bm{A}_{n}$ are within O($n^{-2/3}$) of their classical position with high probability [EYY10]. More generally, Tao and Vu have established the universality of a result due to Gustavsson [Gus05] in the complex Gaussian Wigner case: $(\log n)^{-1/2}(\sqrt{n}\lambda_{j}-n\gamma_{j})$ is asymptotically normally distributed [TV11]. Further, they have shown that eigenvalues in the bulk of the spectrum ($j=\Omega(n)$) of a Wigner matrix satisfy

for some universal constant $c>0$ [TV10a].

In contrast to the asymptotic theory, which remains to a large extent driven by the study of particular classes of random matrices, the nonasymptotic theory has developed as a collection of techniques for addressing the behavior of a broad range of random matrices. The nonasymptotic theory has its roots in geometric functional analysis in the 1970s, where random matrices were used to investigate the local properties of Banach spaces [LM93, SD01, Ver10]. Since then, the nonasymptotic theory has found applications in areas including theoretical computer science [Ach03, Vem04, SS08], machine learning [DM05], optimization [Nem07, So09], and numerical linear algebra [DM10, HMT11, Mah11].

As is the case in the asymptotic theory, the sharpest and most comprehensive results available in the nonasymptotic theory concern the behavior of Gaussian matrices. The amenability of the Gaussian distribution makes it possible to obtain results such as Szarek’s nonasymptotic analog of the Wigner semicircle theorem for Gaussian matrices [Sza90] and Chen and Dongarra’s bounds on the condition number of Gaussian matrices [CD05]. The properties of less well-behaved random matrices can sometimes be related back to those of Gaussian matrices using probabilistic tools, such as symmetrization; see, e.g., the derivation of Latała’s bound on the norms of zero-mean random matrices [Lat05].

More generally, bounds on extremal eigenvalues can be obtained from knowledge of the moments of the entries. For example, the smallest singular value of a square matrix with i.i.d. zero-mean subgaussian entries with unit variance is O($n^{-1/2}$) with high probability [RV08]. Concentration of measure results, such as Talagrand’s concentration inequality for product spaces [Tal95], have also contributed greatly to the nonasymptotic theory. We mention in particular the work of Achlioptas and McSherry on randomized sparsification of matrices [AM01, AM07], that of Meckes on the norms of random matrices [Mec04], and that of Alon, Krivelevich and Vu [AKV02] on the concentration of the largest eigenvalues of random symmetric matrices, all of which are applications of Talagrand’s inequality. In cases where geometric information on the distribution of the random matrices is available, the tools of empirical process theory—such as the generic chaining, also due to Talagrand [Tal05]—can be used to convert this geometric information into information on the spectra. One natural example of such a case consists of matrices whose rows are independently drawn from a log-concave distribution [MP06, ALPTJ11].

The noncommutative Khintchine inequality (NCKI), which bounds the moments of the norm of a sum of fixed matrices modulated by random signs [LP86, LPP91], is a widely used tool in the nonasymptotic theory. Despite its power, the NCKI is unwieldy. To use it, one must reduce the problem to a suitable form by applying symmetrization and decoupling arguments and exploiting the equivalence between moments and tail bounds. It is often more convenient to apply the NCKI in the guise of a lemma, due to Rudelson [Rud99], that provides an analog of the law of large numbers for sums of rank-one matrices. This result has found many applications, including column-subset selection [RV07] and the fast approximate solution of least-squares problems [DMMS11]. The NCKI and its corollaries do not always yield sharp results because parasitic logarithmic factors arise in many settings.

The current paper is ultimately based on the influential work of Ahlswede and Winter [AW02]. This line of research leads to explicit tail bounds for the maximum eigenvalue of a sum of random matrices. These probability inequalities parallel the classical scalar tail bounds due to Bernstein and others. Matrix probability inequalities allow us to obtain valuable information about the maximum eigenvalue of a random matrix with very little effort. Furthermore, they apply to a wide variety of random matrices. We note, however, that matrix probability inequalities can lead to parasitic logarithmic factors similar to those that emerge from the NCKI.

Major contributions to the literature on matrix probability inequalities include the papers [CM08, Rec09, Gro11]. We emphasize two works of Oliveira [Oli09, Oli10] that go well beyond earlier research. The sharpest current results appear in the works of Tropp [Tro11c, Tro11b, Tro11a]. Recently, Hsu, Kakade, and Zhang [HKZ11] have modified Tropp’s approach to establish matrix probability inequalities that depend on an intrinsic dimension parameter, rather than the ambient dimension.

2. Outline

In section 2, we introduce the notation used in this paper and state a convenient version of the Courant–Fischer theorem. In section 3, we use the Courant–Fischer theorem to extend the Laplace transform technique from [Tro11c] to apply to all the eigenvalues of self-adjoint matrices, thereby obtaining the minimax Laplace transform. We apply this technique in sections 4 and 5 to develop eigenvalue analogs of the classical Chernoff and Bernstein bounds. The final two sections illustrate, using two familiar problems, that the minimax Laplace technique gives us significantly more information on the spectra of random matrices than current approaches. In section 6, we use the Chernoff bounds to quantify the effects of column sparsification on all the singular values of matrices with orthogonal rows. In section 7, we consider the question of how fast, in relative error, the eigenvalues of empirical covariance matrices converge.

Background and Notation

We establish the notation used in the sequel and state a convenient version of the Courant–Fischer theorem.

One of our central tools is the variational characterization of the eigenvalues of a self-adjoint matrix given by the Courant–Fischer theorem. For integers $d$ and $n$ satisfying $1\leq d\leq n$, the complex Stiefel manifold

Let $\bm{A}$ be a self-adjoint matrix with dimension $n$. Then

Tail Bounds For Interior Eigenvalues

In this section we develop a generic bound on the tail probabilities of eigenvalues of sums of independent, random, self-adjoint matrices. We establish this bound by supplementing the matrix Laplace transform methodology of [Tro11c] with Proposition 2.1 and a new result, due to Lieb and Seiringer [LS05], on the concavity of a certain trace function on the cone of positive-definite matrices.

First we observe that the Courant–Fischer theorem allows us relate the behavior of the $k$th eigenvalue of a matrix to the behavior of the largest eigenvalue of an appropriate compression of the matrix.

Let $\theta$ be a fixed positive number. Then

The first identity follows from the positive homogeneity of eigenvalue maps and the second from the monotonicity of the scalar exponential function. The final two relations are Markov’s inequality and (2.1).

To continue, we need to bound the expectation. Interchange the order of the exponential and the minimum; then apply the spectral mapping theorem to see that

The first inequality is Jensen’s. The second inequality follows because the exponential of a self-adjoint matrix is positive definite, so its largest eigenvalue is smaller than its trace.

Combine these observations and take the infimum over all positive $\theta$ to complete the argument. ∎

We are interested in the case where the matrix $\bm{X}$ in Theorem 3.1 can be expressed as a sum of independent random matrices. In this case, we use the following result to develop the right-hand side of the Laplace transform bound (3.1).

Consider a finite sequence $\{\bm{X}_{j}\}$ of independent, random, self-adjoint matrices with dimension $n$ and a sequence $\{\bm{A}_{j}\}$ of fixed self-adjoint matrices with dimension $n$ that satisfy the relations

Theorem 3.2 is an extension of Lemma 3.4 of [Tro11c], which establishes the special case (3.4). The proof depends upon a recent result due to Lieb and Seiringer [LS05, Thm. 3] that extends Lieb’s earlier result [Lie73, Thm. 6].

First, note that (3.2) and the operator monotonicity of the matrix logarithm yield the following inequality for each $k$:

The first inequality follows from Proposition 3.1 and Jensen’s inequality, and the second depends on (3.5) and the monotonicity of the trace exponential. Iterate this argument to complete the proof. ∎

Our main result follows from combining Theorem 3.1 and Theorem 3.2.

Consider a finite sequence $\{\bm{X}_{j}\}$ of independent, random, self-adjoint matrices with dimension $n$, and let $k\leq n$ be an integer.

Let $\{\bm{A}_{j}\}$ be a sequence of self-adjoint matrices that satisfy the semidefinite relations

The first bound in Theorem 3.3 requires less detailed information on how compression affects the summands but correspondingly does not give as sharp results as the second.

In the following two sections, we use the minimax Laplace transform method to derive Chernoff and Bernstein inequalities for the interior eigenvalues of a sum of independent random matrices. Tail bounds for the eigenvalues of matrix Rademacher and Gaussian series, eigenvalue Hoeffding, and matrix martingale eigenvalue tail bounds can all be derived in a similar manner; see [Tro11c] for relevant details.

Chernoff bounds

Classical Chernoff bounds establish that the tails of a sum of independent nonnegative random variables decay subexponentially. [Tro11c] develops Chernoff bounds for the maximum and minimum eigenvalues of a sum of independent positive-semidefinite matrices. We extend this analysis to study the interior eigenvalues.

The sequence $\{\bm{X}_{j}\}$ associated with $\Psi$ will always be clear from context. We have the following result.

Consider a finite sequence $\{\bm{X}_{j}\}$ of independent, random, positive-semidefinite matrices with dimension $n.$ Given an integer $k\leq n$, define

where $\Psi$ is a function that satisfies (4.1).

If it is difficult to estimate $\Psi(\bm{V}_{+})$ or $\Psi(\bm{V}_{-}),$ one can resort to the weaker estimates

Theorem 4.1 follows from Theorem 3.3 using an appropriate bound on the matrix moment generating functions. The following lemma is due to Ahlswede and Winter [AW02]; see also [Tro11c, Lem. 5.8].

We consider the case where $\Psi(\bm{V}_{+})=1;$ the general case follows by homogeneity. Define

Bound the trace by the maximum eigenvalue, taking into account the reduced dimension of the summands:

The equality follows from the spectral mapping theorem. Identify the quantity $\mu_{k}$; then combine the last two inequalities to obtain

The right-hand side is minimized when $\theta=\log(1+\delta),$ which gives the desired upper tail bound. ∎

As before, we consider the case where $\Psi(\bm{V}_{-})=1.$ Clearly,

Apply Lemma 4.2 to see that, for $\theta>0,$

Using reasoning analogous to that in the proof of the upper bound, we justify the first of the following inequalities:

The remaining equalities follow from the fact that $-g(\theta)<0$ and the definition of $\mu_{k}.$

The right-hand side is minimized when $\theta=-\log(1-\delta),$ which gives the desired lower tail bound. ∎

Bennett and Bernstein inequalities

The classical Bennett and Bernstein inequalities use the variance or knowledge of the moments of the summands to control the probability that a sum of independent random variables deviates from its mean. In [Tro11c], matrix Bennett and Bernstein inequalities are developed for the extreme eigenvalues of self-adjoint random matrix sums. We establish that the interior eigenvalues satisfy analogous inequalities.

The sequence $\{\bm{X}_{j}\}$ associated with $\Psi$ will always be clear from context.

Consider a finite sequence $\{\bm{X}_{j}\}$ of independent, random, self-adjoint matrices with dimension $n$, all of which have zero mean. Given an integer $k\leq n$, define

where the function $h(u)=(1+u)\log(1+u)-u$ for $u\geq 0.$ The function $\Psi$ satisfies (5.1) above.

Results (i) and (ii) are, respectively, matrix analogs of the classical Bennett and Bernstein inequalities. As in the scalar case, the Bennett inequality reflects a Poisson-type decay in the tails of the eigenvalues. The Bernstein inequality states that small deviations from the eigenvalues of the expected matrix are roughly normally distributed while larger deviations are subexponential. The split Bernstein inequalities (iii) make explicit the division between these two regimes.

As stated, Theorem 5.1 estimates the probability that the eigenvalues of a sum are large. Using the identity

Theorem 5.1 can be applied to estimate the probability that eigenvalues of a sum are small.

To prove Theorem 5.1, we use the following lemma (Lemma 6.7 in [Tro11c]) to control the moment generating function of a random matrix with bounded maximum eigenvalue.

The maximum eigenvalue in this expression equals $\sigma_{k}^{2}$, thus

The Bennett inequality (i) follows by substituting $\theta=\log(1+t/\sigma_{k}^{2})$ into the right-hand side and simplifying.

The Bernstein inequality (ii) is a consequence of (i) and the fact that

which can be established by comparing derivatives.

The subgaussian and subexponential portions of the split Bernstein inequalities (iii) are verified through algebraic comparisons on the relevant intervals. ∎

Occasionally, as in the application in section 7 to the problem of covariance matrix estimation, one desires a Bernstein-type tail bound that applies to summands that do not have bounded maximum eigenvalues. In this case, if the moments of the summands satisfy sufficiently strong growth restrictions, one can extend classical scalar arguments to obtain results such as the following Bernstein bound for subexponential matrices.

Consider a finite sequence $\{\bm{X}_{j}\}$ of independent, random, self-adjoint matrices with dimension $n$, all of which satisfy the subexponential moment growth condition

where $B$ is a positive constant and $\bm{\Sigma}_{j}^{2}$ are positive-semidefinite matrices. Given an integer $k\leq n$, set

This result is an extension of [Tro11c, Theorem 6.2], which, in turn, generalizes a classical scalar argument [DG98].

As with the other matrix inequalities, Theorem 5.3 follows from an application of Theorem 3.3 and appropriate semidefinite bounds on the moment generating functions of the summands. Thus, the key to the proof lies in exploiting the moment growth conditions of the summands to majorize their moment generating functions. The following lemma, a trivial extension of Lemma 6.8 in [Tro11c], provides what we need.

Let $\bm{X}$ be a random self-adjoint matrix satisfying the subexponential moment growth conditions

We note that $\bm{X}_{j}$ satisfies the growth condition

if and only if the scaled matrix $\bm{X}_{j}/B$ satisfies

Thus, by rescaling, it suffices to consider the case $B=1.$ We now do so.

By Lemma 5.4, the moment generating functions of the summands satisfy

where $g(\theta)=\theta^{2}/(2-2\theta).$ Now we apply Theorem 3.3(i):

To achieve the final simplification, we identified $\mu_{k}$ and $\sigma_{k}^{2}.$ Now, select $\theta=t/(t+\sigma_{k}^{2}).$ Then simplication gives the Bernstein inequality (i).

Algebraic comparisons on the relevant intervals yield the split Bernstein inequalities (ii). ∎

An application to column subsampling

As an application of our Chernoff bounds, we examine how sampling columns from a matrix with orthonormal rows affects the spectrum. This question has applications in numerical linear algebra and compressed sensing. The special cases of the maximum and minimum eigenvalues have been studied in the literature [Tro08, RV07]. The limiting spectral distributions of matrices formed by sampling columns from similarly structured matrices have also been studied: the results of [GH10] apply to matrices formed by sampling columns from any fixed orthogonal matrix, and [Far10] studies matrices formed by sampling columns and rows from the discrete Fourier transform matrix. We mention in particular [Rud99], the main result of which provides a uniform bound on the tails of all singular values of the sampled matrix. The theorem proven in this section provides bounds which reflect the differences in the tails of the individual singular values, and thus can be viewed as an elaboration of the result in [Rud99].

Let $\bm{U}$ be an $n\times r$ matrix with orthonormal rows. We model the sampling operation using a random diagonal matrix $\bm{D}$ whose entries are independent $\text{Bern}(p)$ random variables. Then the random matrix

can be interpreted as a random column submatrix of $\bm{U}$ with an average of $pr$ nonzero columns. Our goal is to study the behavior of the spectrum of $\widehat{\bm{U}}.$

Recall that the $j$th column of $\bm{U}$ is written $\bm{u}_{j}.$ Consider the following coherence-like quantity associated with $\bm{U}:$

There does not seem to be a simple expression for $\tau_{k}.$ However, by choosing $\bm{V}^{*}$ to be the restriction to an appropriate $k$-dimensional coordinate subspace, we see that $\tau_{k}$ always satisfies

The following theorem shows that the behavior of $s_{k}(\widehat{\bm{U}}),$ the $k$th singular value of $\widehat{\bm{U}},$ can be explained in terms of $\tau_{k}.$

Let $\bm{U}$ be an $n\times r$ matrix with orthonormal rows, and let $p$ be a sampling probability. Define the sampled matrix $\widehat{\bm{U}}$ according to (6.1), and the numbers $\{\tau_{k}\}$ according to (6.2). Then, for each $k=1,\ldots,n,$

where $\bm{u}_{j}$ is the $j$th column of $\bm{U}$ and $d_{j}\sim\text{Bern}(p).$ Compute

Covariance Estimation

We conclude with an extended example that illustrates how this circle of ideas allows one to answer interesting statistical questions. Specifically, we investigate the convergence of the individual eigenvalues of sample covariance matrices, with errors measured in relative precision.

An important challenge is to determine how many samples are needed to ensure that the empirical covariance estimator has a fixed relative accuracy in the spectral norm. That is, given a fixed $\varepsilon,$ how large must $n$ be so that

This estimation problem has been studied extensively. It is now known that for distributions with a finite second moment, $\Omega(p\log p)$ samples suffice [Rud99], and for log-concave distributions, $\Omega(p)$ samples suffice [ALPTJ11]. More broadly, Vershynin [Ver11] conjectures that, for distributions with finite fourth moment, $\Omega(p)$ samples suffice; he establishes this result to within iterated log factors. In [SV11], Srivastava and Vershynin establish that $\Omega(p)$ samples suffice for distributions which have finite $2+\varepsilon$ moments, for some $\varepsilon>0,$ and satisfy an additional regularity condition.

In this section, we derive a relative approximation bound for each eigenvalue of $\bm{C}$ that allows us to confirm this intuition. For simplicity we assume the samples are drawn from a $\mathcal{N}(\bm{0},\bm{C})$ distribution where $\bm{C}$ is full-rank, but the arguments can be extended to cover other distributions.

Write $\lambda_{k}$ for the $k$th eigenvalue of $\bm{C}$, and write $\hat{\lambda}_{k}$ for the $k$th eigenvalue of $\widehat{\bm{C}}_{n}.$ Then for $k=1,\ldots,p,$

The following corollary provides an answer to our question about relative error estimates.

Let $\lambda_{k}$ and $\hat{\lambda}_{k}$ be as in Theorem 7.1. Then

The first bound in Corollary 7.2 tells us how many samples are needed to ensure that $\hat{\lambda}_{k}$ does not overestimate $\lambda_{k}.$ Likewise, the second bound tells us how many samples ensure that $\hat{\lambda}_{k}$ does not underestimate $\lambda_{k}.$

Corollary 7.2 suggests that the relationship of $\hat{\lambda}_{k}$ to $\lambda_{k}$ is determined by the spectrum of $\bm{C}$ in the following manner. When the eigenvalues below $\lambda_{k}$ are small compared with $\lambda_{k}$, the quantity

is small, and so $\hat{\lambda}_{k}$ is not likely to overestimate $\lambda_{k}$. Similarly, when the eigenvalues above $\lambda_{k}$ are comparable with $\lambda_{k}$, the quantity

is small, and so $\hat{\lambda}_{k}$ is not likely to underestimate $\lambda_{k}$.

We now have everything needed to establish Theorem 1.1.

Assuming the stated decay condition, that

The results in Theorem 7.1 and Corollary 7.2 also apply when $\bm{C}$ is rank-deficient: simply replace each occurence of the dimension $p$ in the bounds with $\operatorname{rank}(\bm{C}).$

We now prove Theorem 7.1. This result requires supporting lemmas; we defer their proofs until after a discussion of extensions to Theorem 7.1.

We study the error $|\lambda_{k}(\widehat{\bm{C}}_{n})-\lambda_{k}(\bm{C})|.$ To apply the methods developed in this paper, we pass to a question about the eigenvalues of a difference of two matrices. The first lemma accomplishes this goal by compressing both the population covariance matrix and the sample covariance matrix to a fixed invariant subspace of the population covariance matrix.

We apply this result with $\bm{A}=\bm{C}$ and $\bm{X}=\widehat{\bm{C}}_{n}.$ Because $\widehat{\bm{C}}_{n}$ is unbounded, we apply Theorem 5.3 to handle the estimates in (7.2) and (7.3). To use this theorem, we need the following moment growth estimate for rank-one Wishart matrices.

Let $\bm{\xi}\sim\mathcal{N}(\bm{0},\bm{G}).$ Then for any integer $m\geq 2,$

With these preliminaries addressed, we prove Theorem 7.1.

The factor $n$ comes from the normalization of the sample covariance matrix.

The covariance matrix of $\bm{\eta}_{j}$ is $\bm{C},$ so that of $\bm{W}_{+}^{*}\bm{\eta}_{j}$ is $\bm{W}_{+}^{*}\bm{C}\bm{W}_{+}.$ Apply Lemma 7.4 to verify that $\bm{W}_{+}^{*}\bm{\eta}_{j}\bm{\eta}_{j}\bm{W}_{+}$ satisfies the subexponential moment growth bound required by Theorem 5.3 with

In fact, $\bm{W}_{+}^{*}\bm{C}\bm{W}_{+}$ is the compression of $\bm{C}$ to the invariant subspace corresponding with its bottom $p-k+1$ eigenvalues, so

We are concerned with the maximum eigenvalue of the sum in (7.4), so we take $\bm{V}_{+}=\mathbf{I}$ in the statement of Theorem 5.3 to find that

It follows from the subgaussian branch of the split Bernstein inequality of Theorem 5.3 that

when $t\leq 4n\lambda_{k}(\bm{C}).$ This provides the desired bound on the probability that $\lambda_{k}(\widehat{\bm{C}}_{n})$ overestimates $\lambda_{k}(\bm{C}).$ ∎

The factor $n$ comes from the normalization of the sample covariance matrix.

The covariance matrix of $\bm{\eta}_{j}$ is $\bm{C},$ so that of $\bm{W}_{-}^{*}\bm{\eta}_{j}$ is $\bm{W}_{-}^{*}\bm{C}\bm{W}_{-}.$ Apply Lemma 7.4 to verify that for any integer $m\geq 2,$

Thus, $\bm{W}_{-}^{*}(-\bm{\eta}_{j}\bm{\eta}_{j}^{*})\bm{W}_{-}$ satisfies the subexponential moment growth bound required by Theorem 5.3 with

In fact, $\bm{W}_{-}^{*}\bm{C}\bm{W}_{-}$ is the compression of $\bm{C}$ to the invariant subspace corresponding with its top $k$ eigenvalues, so

We are concerned with the maximum eigenvalue of the sum in (7.5), so we take $\bm{V}_{+}=\mathbf{I}$ in the statement of Theorem 5.3 to find that

It follows from the subgaussian branch of the split Bernstein inequality of Theorem 5.3 that

when $t\leq 4n\lambda_{1}(\bm{C}).$ This provides the desired bound on the probability that $\lambda_{k}(\widehat{\bm{C}}_{n})$ underestimates $\lambda_{k}(\bm{C}).$ ∎

2. Extensions of Theorem 7.1

Results analogous to Theorem 7.1 can be established for other distributions. If the distribution is bounded, the possibility that $\hat{\lambda}_{k}$ deviates above or below $\lambda_{k}$ can be controlled using the Bernstein inequality of Theorem 5.1. If the distribution is unbounded but has matrix moments that satisfy a sufficiently nice growth condition, the probability that $\hat{\lambda}_{k}$ deviates below $\lambda_{k}$ as well as the probability that it deviates above $\lambda_{k}$ can be bounded using a Bernstein inequality analogous to that in Theorem 5.3.

Theorem 7.1 controls the error in the $k$th sample eigenvalue in terms of all the eigenvalues of the covariance matrix, so it is most useful when the eigenvalues of the covariance matrix satisfy decay conditions such as those given in the statement of Theorem 1.1. If such conditions are not satisfied, the results of [ALPTJ11] on the convergence of empirical covariance matrices of isotropic log-concave random vectors lead to tighter bounds on the probabilities that $\hat{\lambda}_{k}$ overestimates or underestimates $\lambda_{k}.$

To see the relevance of the results in [ALPTJ11], first observe the following consequence of the subadditivity of the maximum eigenvalue mapping:

In conjunction with (7.2), this gives us the following control on the probability that $\lambda_{k}(\bm{X})$ overestimates $\lambda_{k}(\bm{A}):$

In our application, $\bm{X}$ is the empirical covariance matrix and $\bm{A}$ is the actual covariance matrix. The spectral norm dominates the maximum eigenvalue, so

where $\bm{S}$ is the square root of $\bm{W}_{+}^{*}\bm{C}\bm{W}_{+}.$ Now factor out $\bm{S}^{2}$ and identify $\lambda_{k}(\bm{C})=\|\bm{S}^{2}\|$ to obtain

Note that if $\bm{\eta}$ is drawn from a $\mathcal{N}(\bm{0},\bm{C})$ distribution, then the covariance matrix of the transformed sample $\bm{S}^{-1}\bm{W}_{+}^{*}\bm{\eta}$ is the identity:

Similarly, for more general distributions, the bounds on the probability of $\hat{\lambda}_{k}$ overestimating or underestimating $\lambda_{k}$ can be tightened beyond those suggested in Theorem 7.1 by using the results in [ALPTJ11] or [Ver11]. Note, however, that one cannot use knowledge of spectral decay to sharpen the results obtained from [ALPTJ11] and [Ver11] into estimates like those given in Theorem 1.1.

Finally, we note that the techniques developed in the proof of Theorem 7.1 can be used to investigate the spectrum of the error matrices $\widehat{\bm{C}}_{n}-\bm{C}.$

3. Proofs of the supporting lemmas

We now establish the lemmas used in the proof of Theorem 7.1.

The probability that $\lambda_{k}(\bm{X})$ overestimates $\lambda_{k}(\bm{A})$ is controlled with the sequence of inequalities

We use a related approach to study the probability that $\lambda_{k}(\bm{X})$ underestimates $\lambda_{k}(\bm{A}).$ Our choice of $\bm{W}_{-}$ implies that

This establishes the bounds on the probabilities of $\lambda_{k}(\bm{X})$ deviating above or below $\lambda_{k}(\bm{A}).$ ∎

Factor the covariance matrix of $\bm{\xi}$ as $\bm{G}=\bm{U\Lambda U}^{*}$ where $\bm{U}$ is orthogonal and $\bm{\Lambda}=\text{diag}(\lambda_{1},\ldots,\lambda_{p})$ is the matrix of eigenvalues of $\bm{G}$. Let $\bm{\gamma}$ be a $\mathcal{N}(\bm{0},\mathbf{I}_{p})$ random variable. Then $\bm{\xi}$ and $\bm{U\Lambda}^{1/2}\bm{\gamma}$ are identically distributed, so

Consider the $(i,j)$ entry of the bracketed matrix in (7.6):

From this expression, and the independence of the Gaussian variables $\{\gamma_{i}\},$ we see that this matrix is diagonal.

To bound the diagonal entries, use a multinomial expansion to further develop the sum in (7.7) for the $(i,i)$ entry:

Denote the $L_{r}$ norm of a random variable $X$ by

The second inequality is the triangle inequality for $L_{r}$ norms. Now we reverse the multinomial expansion to see that the diagonal terms satisfy the inequality

Combine this result with (7.8) to see that

Complete the proof by using this estimate in (7.6).

References