How close are the eigenvectors and eigenvalues of the sample and actual covariance matrices?

TL;DR

This paper establishes non-asymptotic bounds on how closely the eigenvectors and eigenvalues of sample covariance matrices approximate those of the true covariance matrix, depending on the number of samples and distribution properties.

Contribution

It provides new non-asymptotic concentration bounds for eigenvectors, eigenspaces, and eigenvalues for a broad class of distributions, including finite second moment and bounded support.

Findings

Inner product between sample and true eigenvectors decreases with eigenvalue distance

Conditions identified for distinguishing principal components with limited samples

Bounds applicable to distributions with finite second moments and bounded support

Abstract

How many samples are sufficient to guarantee that the eigenvectors and eigenvalues of the sample covariance matrix are close to those of the actual covariance matrix? For a wide family of distributions, including distributions with finite second moment and distributions supported in a centered Euclidean ball, we prove that the inner product between eigenvectors of the sample and actual covariance matrices decreases proportionally to the respective eigenvalue distance. Our findings imply non-asymptotic concentration bounds for eigenvectors, eigenspaces, and eigenvalues. They also provide conditions for distinguishing principal components based on a constant number of samples.

Figures

Equations

$$\displaystyle m=O\Bigg{(}\frac{k_{j}^{2}}{(\lambda_{i}-\lambda_{j})^{2}}\Bigg{)}\quad\text{and}\quad m=O\Bigg{(}\frac{k_{i}^{2}}{\lambda_{i}^{2}}\Bigg{)}$$

$$\tilde{C}=\sum _{p=1}^m\frac{(x_p-\bar{x})(x_p-\bar{x}{)}^{\ast }}{m},$$

$$\left|⟨{\tilde{u}}_j,u_j⟩\right|\le \frac{2\Vert \delta C{u}_j{\Vert }_2}{\left|{\lambda }_i-{\lambda }_j\right|},$$

$$\displaystyle{\mathbf{P}\left(|\langle\widetilde{u}_{i},u_{j}\rangle|\geq t\right)}\leq\frac{1}{m}\,\Big{(}\frac{2k_{j}}{t\,|\lambda_{i}-\lambda_{j}|}\Big{)}^{2},$$

$$\mathbf{P}\left(\frac{|{\tilde{\lambda }}_i-{\lambda }_i|}{{\lambda }_i}\ge t\right)\le \frac{1}{m}{\left(\frac{k_i}{{\lambda }_{i}t}\right)}^2.$$

$$\mathbf{P}\left(|⟨{\tilde{u}}_i,u_j⟩|\ge t\right)\le \text{exp}\left(1-\frac{cm{\Phi }_{ij}(t{)}^2}{{\lambda }_j{\Vert x\Vert }_{{\Psi }_2}^2}\right),$$

$$(C+\delta C)(u_i+\delta u_i)=({\lambda }_i+\delta {\lambda }_i)(u_i+\delta u_i).$$

$$C\delta u_i+\delta C{u}_i+\delta C\delta u_i={\lambda }_i\delta u_i+\delta {\lambda }_i{u}_i+\delta {\lambda }_i\delta u_i,$$

$$\sum _{\ell =1}^n{\beta }_{ij}{u}_j^{\ast }C{u}_{\ell }+u_j^{\ast }\delta C{u}_i+\sum _{\ell =1}^n{\beta }_{ij}{u}_j^{\ast }\delta C{u}_{\ell }={\lambda }_i\sum _{\ell =1}^n{\beta }_{ij}{u}_j^{\ast }{u}_{\ell }+\delta {\lambda }_i{u}_j^{\ast }{u}_i+\delta {\lambda }_i\sum _{\ell =1}^n{\beta }_{ij}{u}_j^{\ast }{u}_{\ell }$$

$$\delta {\lambda }_i{u}_j^{\ast }{u}_i+({\lambda }_i+\delta {\lambda }_i-{\lambda }_j){\beta }_{ik}=u_j^{\ast }\delta C{u}_i+\sum _{\ell =1}^n{\beta }_{ij}{u}_j^{\ast }\delta C{u}_{\ell }.$$

$$\delta {\lambda }_i{u}_j^{\ast }{u}_i+({\tilde{\lambda }}_i-{\lambda }_j)({\tilde{u}}_i^{\ast }{u}_j-u_i^{\ast }{u}_j)=u_j^{\ast }\delta C{u}_i+\sum _{\ell =1}^n({\tilde{u}}_i^{\ast }{u}_{\ell })u_j^{\ast }\delta C{u}_{\ell }-u_j^{\ast }\delta C{u}_i.$$

$$({\tilde{\lambda }}_i-{\lambda }_j{)}^2({\tilde{u}}_i^{\ast }{u}_j{)}^2={\left(\sum _{\ell =1}^n({\tilde{u}}_i^{\ast }{u}_{\ell })(u_j^{\ast }\delta C{u}_{\ell })\right)}^2.$$

$$({\tilde{\lambda }}_i-{\lambda }_j{)}^2({\tilde{u}}_i^{\ast }{u}_j{)}^2$$

$$\sum _{\ell =1}^n(u_j^{\ast }\delta C{u}_{\ell }{)}^2$$

$$\left|⟨{\tilde{u}}_j,u_j⟩\right|\le \frac{2\Vert \delta C{u}_j{\Vert }_2}{\left|{\lambda }_i-{\lambda }_j\right|},$$

$$(\delta {\lambda }_i+{\lambda }_i-{\lambda }_j)({\tilde{u}}_i^{\ast }{u}_j)=\sum _{\ell =1}^n({\tilde{u}}_i^{\ast }{u}_{\ell })(u_j^{\ast }\delta C{u}_{\ell }).$$

$$\left|{\tilde{u}}_i^{\ast }{u}_j\right|$$

$$\left|{\tilde{u}}_i^{\ast }{u}_j\right|\le \frac{\Vert \delta C{u}_j{\Vert }_2}{|{\tilde{\lambda }}_i-{\lambda }_j|}\le \frac{\Vert \delta C{u}_j{\Vert }_2}{|{\lambda }_i-{\lambda }_j|}<\frac{2\Vert \delta C{u}_j{\Vert }_2}{\left|{\lambda }_i-{\lambda }_j\right|}.$$

$$\displaystyle{\mathbf{P}\left(|\langle\widetilde{u}_{i},u_{j}\rangle|\geq t\right)}\leq\frac{1}{m}\,\Big{(}\frac{2\,k_{j}}{t\,|\lambda_{i}-\lambda_{j}|}\Big{)}^{2}$$

$$\mathbf{P}\left(|X-\alpha |\ge t\right)\le \frac{\mathbf{Var}\left[X\right]+(\mathbf{E}\left[X\right]-\alpha {)}^2}{t^2}.$$

$$\mathbf{P}\left(|⟨{\tilde{u}}_i,u_j⟩|\ge t\right)$$

$$\mathbf{E}\left[\Vert \delta C{u}_j{\Vert }_2^2\right]$$

$$=\mathbf{E}\left[u_j^{\ast }(\tilde{C}-{\lambda }_j)(\tilde{C}-{\lambda }_j)u_j\right]$$

$$=\mathbf{E}\left[u_j^{\ast }{\tilde{C}}^2{u}_j\right]+{\lambda }_j^2-2{\lambda }_j{u}_j^{\ast }\mathbf{E}\left[\tilde{C}\right]u_j$$

$$=\mathbf{E}\left[u_j^{\ast }{\tilde{C}}^2{u}_j\right]-{\lambda }_j^2.$$

$$\mathbf{E}\left[u_j^{\ast }{\tilde{C}}^2{u}_j\right]$$

$$=\sum _{p\ne q}{u}_j^{\ast }\frac{\mathbf{E}\left[x_p{x}_p^{\ast }\right]\mathbf{E}\left[x_q{x}_q^{\ast }\right]}{m^2}{u}_j+\sum _{p=1}^m{u}_j^{\ast }\frac{\mathbf{E}\left[x_p{x}_p^{\ast }{x}_p{x}_p^{\ast }\right]}{m^2}{u}_j$$

$$=\frac{m(m-1)}{m^2}{\lambda }_j^2+\frac{1}{m}{u}_j^{\ast }\mathbf{E}\left[x{x}^{\ast }x{x}^{\ast }\right]u_j$$

$$=(1-\frac{1}{m}){\lambda }_j^2+\frac{1}{m}{u}_j^{\ast }\mathbf{E}\left[x{x}^{\ast }x{x}^{\ast }\right]u_j$$

$$\mathbf{E}\left[\Vert \delta C{u}_j{\Vert }_2^2\right]$$

$$=\frac{u_j^{\ast }\mathbf{E}\left[x{x}^{\ast }x{x}^{\ast }\right]u_j-{\lambda }_j^2}{m}=\frac{\mathbf{E}\left[\Vert x{x}^{\ast }{u}_j{\Vert }_2^2\right]-{\lambda }_j^2}{m}.$$

$$\mathbf{P}\left(\sum _{i\ne j}{w}_{ij}⟨{\tilde{u}}_i,u_j{⟩}^2>t\right)\le \sum _{i\ne j}\frac{4w_{ij}{k}_j^2}{mt({\lambda }_i-{\lambda }_j{)}^2},$$

Full text

1 Introduction

The covariance matrix $C$ of an $n$-dimensional distribution is an integral part of data analysis, with numerous occurrences in machine learning and signal processing. It is therefore crucial to understand how close it is to the sample covariance, i.e., the matrix $\widetilde{C}$ estimated from a finite number of samples $m$. Following developments in the tools for the concentration of measure, Vershynin showed that a sample size of $m=O(n)$ is up to iterated logarithmic factors sufficient for all distributions with finite fourth moment supported in a centered Euclidean ball of radius $O(\sqrt{n})$ [26]. Similar results hold also for sub-exponential distributions [1] and distributions with finite second moment [20].

We take an alternative standpoint and ask if we can do better when only a subset of the spectrum is of interest. Concretely, our objective is to characterize how many samples are sufficient to guarantee that an eigenvector and/or eigenvalue of the sample and actual covariance matrices are, respectively, sufficiently close. Our approach is motivated by the observation that methods that utilize the covariance commonly prioritize the estimation of principal eigenspaces. For instance, in (local) principal component analysis we are usually interested in the first few eigenvectors [9, 18], whereas when reducing the dimension of a distribution one commonly projects it to the span of the first few eigenvectors [16, 12].

Our finding is that the “spectral leaking” occurring in the eigenvector estimation is strongly localized w.r.t. the eigenvalue axis. In other words, the eigenvector $\widetilde{u}_{i}$ of the sample covariance is less far likely to lie in the span of an eigenvector $u_{j}$ of the actual covariance when the eigenvalue distance $|\lambda_{i}-\lambda_{j}|$ is large and the concentration of the distribution in the direction of $u_{j}$ is small. This phenomenon agrees with the intuition that principal components of high variance are easier to estimate, exactly because they are more likely to appear in the samples of the distribution. In addition, it suggests that it might be possible to obtain good estimates of well separated principal eigenspaces from fewer than $n$ samples.

We provide a mathematical argument confirming this phenomenon. Under fairly general conditions, we prove that

[TABLE]

samples are asymptotically almost surely (a.a.s). sufficient to guarantee that $|\langle\widetilde{u}_{i},u_{j}\rangle|$ and $|\delta\lambda_{i}|/\lambda_{i}$, respectively, is small for all distributions with finite second moment. Here, $k_{j}^{2}$ is a measure of the kurtosis of the distribution in the direction of $u_{j}$. We also attain a high probability bound for distributions supported in a centered and scaled Euclidean ball, and show how our results can be used to characterize the sensitivity of principal component analysis to a limited sample set.

To the best of our knowledge, these are the first non-asymptotic results concerning the eigenvectors of the sample covariance of distributions with finite second moment. Previous studies have intensively investigated the limiting distribution of the eigenvalues of a sample covariance matrix [24, 4], such as the smallest and largest eigenvalues [6] and the eigenvalue support [7]. Eigenvectors and eigenprojections have attracted less attention; the main research thrust entails using tools from the theory of large-dimensional matrices to characterize limiting distributions [3, 13, 22, 5] and it has limited applicability in the non-asymptotic setting where the sample size $m$ is small and $n$ cannot be arbitrary large.

Differently, our arguments follow from a combination of techniques from perturbation analysis and non-asymptotic concentration of measure. Moreover, in contrast to standard perturbation bounds [11, 27] commonly used to reason about eigenspaces [14, 15], they can be used to characterize weighted linear combinations of $\langle\widetilde{u}_{i},u_{j}\rangle^{2}$ over $i$ and $j$, and they do not depend on the minimal eigenvalue gap separating two eigenspaces but rather on all eigenvalue differences. The latter renders them particularly amendable to situations where the eigenvalue gap is not significant but the eigenvalue magnitudes decrease sufficiently fast.

Our work also connects to subspace methods, where the signal and noise spaces are separated by an appropriate eigenspace projection. In their recent work, Shaghaghi and Vorobyov characterized the first two moments of the projection error, a result which implies sample estimates [23]. Their results are particularly tight, but are restricted to specific projectors and Normal distributions. Finally, we remark that there exist alternative estimators for the spectrum of the covariance with better asymptotic properties [2, 19]. Instead, we here focus on the standard estimates, i.e., the eigenvalues and eigenvectors of the sample covariance.

2 Problem statement and main results

Let $x\in\mathbb{C}^{n}$ be a sample of a multivariate distribution and denote by $x_{1},x_{2},\ldots,x_{m}$ the $m$ independent samples used to form the sample covariance, defined as

[TABLE]

where $\bar{x}$ is the sample mean. Denote by $u_{i}$ the eigenvector of $C$ associated with eigenvalue $\lambda_{i}$, and correspondingly for the eigenvectors $\widetilde{u}_{i}$ and eigenvalues $\widetilde{\lambda}_{i}$ of $\widetilde{C}$, such that $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$. We ask:

Problem 1**.**

How many samples are sufficient to guarantee that the inner product $|\langle\widetilde{u}_{i},u_{j}\rangle|=|\widetilde{u}_{i}^{*}u_{j}|$ and the eigenvalue gap $|\delta\lambda_{i}|=|\widetilde{\lambda}_{i}-\lambda_{i}|$ is smaller than some constant $t$ with probability larger than $\epsilon$?

Clearly, when asking that all eigenvectors and eigenvalues of the sample and actual covariance matrices are close, we will require at least as many samples as needed to ensure that $\|\widetilde{C}-C\|_{2}\leq t$ [26]. However, we might do better when only a subset of the spectrum is of interest. The reason is that inner products $|\langle\widetilde{u}_{i},u_{j}\rangle|$ possess strong localized structure along the eigenvalue axis. To illustrate this phenomenon, let us consider the distribution constructed by the $n=784$ pixel values of digit ‘1’ in the MNIST database. Figure 1, compares the eigenvectors $u_{j}$ of the covariance computed from all 6742 images, to the eigenvectors $\widetilde{u}_{i}$ of the sample covariance matrices $\widetilde{C}$ computed from a random subset of $m=10$, 100, 500, and $1000$ samples. For each $i=1,4,20,100$, we depict at $\lambda_{j}$ the average of $|\langle\widetilde{u}_{i},u_{j}\rangle|$ over 100 sampling draws. We observe that: (i) The magnitude of $\langle\widetilde{u}_{i},u_{j}\rangle$ is inversely proportional to their eigenvalue gap $|\lambda_{i}-\lambda_{j}|$. (ii) Eigenvector $\widetilde{u}_{j}$ mostly lies in the span of eigenvectors $u_{j}$ over which the distribution is concentrated.

We formalize these statements in two steps.

2.1 Perturbation arguments

First, we work in the setting of Hermitian matrices and notice the following inequality:

Theorem 3.2.

For Hermitian matrices $C$ and $\widetilde{C}=\delta C+C$, with eigenvectors $u_{j}$ and $\widetilde{u}_{i}$ respectively, the inequality

[TABLE]

holds for $\text{sgn}(\lambda_{i}>\lambda_{j})\,2\widetilde{\lambda}_{i}>\text{sgn}(\lambda_{i}>\lambda_{j})(\lambda_{i}+\lambda_{j})$ and $\lambda_{i}\neq\lambda_{j}$.

The above stands out from standard eigenspace perturbation results, such as the sin($\Theta$) Theorem [11] and its variants [14, 15, 27] for three main reasons: First, Theorem 3.2 characterizes the angle between any pair of eigenvectors allowing us to jointly bound any linear combination of inner-products. Though this often proves handy (c.f. Section 5), it was not possible using sin($\Theta$)-type arguments. Second, classical bounds are not appropriate for a probabilistic analysis as they feature ratios of dependent random variables (corresponding to perturbation terms). In the analysis of spectral clustering, this complication was dealt with by assuming that $|{\lambda}_{i}-\lambda_{j}|\leq|\widetilde{\lambda}_{i}-\lambda_{j}|$ [15]. We weaken this condition at a cost of a multiplicative factor of 2. In contrast to previous work, we also prove that the condition is met a.a.s. Third, previous bounds are expressed in terms of the minimal eigenvalue gap between eigenvectors lying in the interior and exterior of the subspace of interest. This is a limiting factor in practice as it renders the results only amenable to situations where there is a very large eigenvalue gap separating the subspaces. The proposed result improves upon this by considering every eigenvalue difference.

2.2 Spectral concentration

The second part of our analysis focuses on the covariance and has a statistical flavor. In particular, we give an answer to Problem 1 for various families of distributions.

In the context of distributions with finite second moment, we prove in Section 4.1 that:

Theorem 4.1.

For any two eigenvectors $\widetilde{u}_{i}$ and $u_{j}$ of the sample and actual covariance respectively, and for any real number $t>0$:

[TABLE]

s.t. the same conditions as Theorem 3.2.

For eigenvalues, we provide the following guarantee:

Theorem 4.2.

For any eigenvalues $\lambda_{i}$ and $\widetilde{\lambda}_{i}$ of $C$ and $\widetilde{C}$, respectively, and for any $t>0$, we have

[TABLE]

Term $k_{j}=(\mathbf{E}\left[\|xx^{*}u_{j}\|_{2}^{2}\right]-\lambda_{j}^{2})^{1/2}$ captures the tendency of the distribution to fall in the span of $u_{j}$: the smaller the tail in the direction of $u_{j}$ the less likely we are going to confuse $\widetilde{u}_{i}$ with $u_{j}$.

For normal distributions, we have that $k_{j}^{2}=\lambda_{j}^{2}+\lambda_{j}\text{tr}(C)$ and the number of samples needed for $|\langle\widetilde{u}_{i},u_{j}\rangle|$ to be small is $m=O({\text{tr}(C)}/{\lambda_{i}^{2}})$ when $\lambda_{j}=O(1)$ and $m=O(\lambda_{i}^{-2})$ when $\quad\lambda_{j}=O(\text{tr}(C)^{-1})$. Thus for normal distributions, principal components $u_{i}$ and $u_{j}$ with $\text{min}\{\lambda_{i}/\lambda_{j},\lambda_{i}\}=\Omega(\text{tr}(C)^{1/2})$ can be distinguished given a constant number of samples. On the other hand, estimating $\lambda_{i}$ with small relative error requires $m=O({\text{tr}(C)}/{\lambda_{i}})$ samples and can thus be achieved from very few samples when $\lambda_{i}$ is large111Though the same cannot be stated about the absolute error $|\delta\lambda_{i}|$, that is smaller for small $\lambda_{i}$..

In Section 4.2, we also give a sharp bound for the family of distributions supported within a ball (i.e., $\|x\|\leq r$ a.s.).

Theorem 4.3.

For sub-gaussian distributions supported within a centered Euclidean ball of radius $r$, there exists an absolute constant $c$, independent of the sample size, such that for any real number $t>0$,

[TABLE]

where $\Phi_{ij}(t)=\frac{|\lambda_{i}-\lambda_{j}|\,t-2\lambda_{j}}{2\,(r^{2}/\lambda_{j}-1)^{1/2}}-2\left\lVert x\right\rVert_{\Psi_{2}}$ and s.t. the same conditions as Theorem 3.2.

Above, $\left\lVert x\right\rVert_{\Psi_{2}}$ is the sub-gaussian norm, for which we usually have $\left\lVert x\right\rVert_{\Psi_{2}}=O(1)$ [25]. As such, the theorem implies that, whenever $\lambda_{i}\gg\lambda_{j}=O(1)$, the sample requirement is with high probability $m=O(r^{2}/{\lambda_{i}^{2}})$.

These theorems solidify our experimental findings shown in Figure 1. Moreover, combined with the analysis of principal component estimation given in Section 5, they provide a concrete characterization of the relation between the spectrum (eigenvectors and eigenvalues) of the sample and actual covariance matrix as a function of the number of samples, the eigenvalue gap, and the distribution properties.

3 Perturbation arguments

Before focusing on the sample covariance matrix, it helps to study $\langle\widetilde{u}_{i},u_{j}\rangle$ in the setting of Hermitian matrices. The presentation of the results is split in three parts. Section 3.1 starts by studying some basic properties of inner products of the form $\langle\widetilde{u}_{i},u_{j}\rangle$, for any $i$ and $j$. The results are used in Section 3.2 to provide a first bound on the angle between two eigenvectors, and refined in Section 3.3.

3.1 Basic observations

We start by noticing an exact relation between the angle of a perturbed eigenvector and the actual eigenvectors of $C$.

Lemma 3.1**.**

For every $i$ and $j$ in $1,2,\ldots,n$, the relation $(\widetilde{\lambda}_{i}-\lambda_{j})\,(\widetilde{u}_{i}^{*}u_{j})=\sum_{\ell=1}^{n}(\widetilde{u}_{i}^{*}u_{\ell})\,(u_{j}^{*}\delta Cu_{\ell})$ holds .

Proof.

The proof follows from a modification of a standard argument in perturbation theory. We start from the definition $\widetilde{C}\,\widetilde{u}_{i}=\widetilde{\lambda}_{i}\,\widetilde{u}_{i}$ and write

[TABLE]

Expanded, the above expression becomes

[TABLE]

where we used the fact that $Cu_{i}=\lambda_{i}u_{i}$ to eliminate two terms. To proceed, we substitute $\delta u_{i}=\sum_{j=1}^{n}\beta_{ij}u_{\ell}$, where $\beta_{ij}=\delta u_{i}^{*}u_{\ell}$, into  (6) and multiply from the left by $u_{j}^{*}$, resulting to:

[TABLE]

Cancelling the unnecessary terms and rearranging, we have

[TABLE]

At this point, we note that $(\lambda_{i}+\delta\lambda_{i}-\lambda_{j})=\widetilde{\lambda}_{i}-\lambda_{j}$ and furthermore that $\beta_{ik}=\widetilde{u}_{i}^{*}u_{j}-u_{i}^{*}u_{j}$. With this in place, equation (8) becomes

[TABLE]

The proof completes by noticing that, in the left hand side, all terms other than $(\widetilde{\lambda}_{i}-\lambda_{j})\,\widetilde{u}_{i}^{*}u_{j}$ fall-off, either due to $u_{i}^{*}u_{j}=0$, when $i\neq k$, or because $\delta\lambda_{i}=\widetilde{\lambda}_{i}-\lambda_{j}$, o.w. ∎

As the expression reveals, $\langle\widetilde{u}_{i},u_{j}\rangle$ depends on the orientation of $\widetilde{u}_{i}$ with respect to all other $u_{\ell}$. Moreover, the angles between eigenvectors depend not only on the minimal gap between the subspace of interest and its complement space (as in the sin($\Theta$) theorem), but on every difference $\widetilde{\lambda}_{i}-\lambda_{j}$. This is a crucial ingredient to a tight bound, that will be retained throughout our analysis.

3.2 Bounding arbitrary angles

We proceed to decouple the inner products.

Theorem 3.1**.**

For any Hermitian matrices $C$ and $\widetilde{C}=\delta C+C$, with eigenvectors $u_{j}$ and $\widetilde{u}_{i}$ respectively, we have that $|\widetilde{\lambda}_{i}-\lambda_{j}|\left\lvert\langle\widetilde{u}_{i},u_{j}\rangle\right\rvert\leq\|\delta C\,u_{j}\|_{2}$.

Proof.

We rewrite Lemma 3.1 as

[TABLE]

We now use the Cauchy-Schwartz inequality

[TABLE]

where in the last step we exploited Lemma 3.2. The proof concludes by taking a square root at both sides of the inequality. ∎

Lemma 3.2**.**

$\sum\limits_{\ell=1}^{n}(u_{j}^{*}\delta Cu_{\ell})^{2}=\|\delta C\,u_{j}\|_{2}^{2}$.

Proof.

We first notice that $u_{j}^{*}\delta Cu_{\ell}$ is a scalar and equal to its transpose. Moreover, $\delta C$ is Hermitian as the difference of two Hermitian matrices. We therefore have that

[TABLE]

matching our claim. ∎

3.3 Refinement

As a last step, we move all perturbation terms to the numerator, at the expense of a multiplicative constant factor.

Theorem 3.2**.**

For Hermitian matrices $C$ and $\widetilde{C}=\delta C+C$, with eigenvectors $u_{j}$ and $\widetilde{u}_{i}$ respectively, the inequality

[TABLE]

holds for $\text{sgn}(\lambda_{i}>\lambda_{j})\,2\widetilde{\lambda}_{i}>\text{sgn}(\lambda_{i}>\lambda_{j})(\lambda_{i}+\lambda_{j})$ and $\lambda_{i}\neq\lambda_{j}$.

Proof.

Adding and subtracting $\lambda_{i}$ from the left side of the expression in Lemma 3.1 gives

[TABLE]

For $\lambda_{i}\neq\lambda_{j}$, the above expression can be re-written as

[TABLE]

Let us examine the right-hand side inequality carefully. Obviously, when the condition $\left\lvert\lambda_{i}-\lambda_{j}\right\rvert\leq 2\left\lvert\delta\lambda_{i}\right\rvert$ is not met, the right clause of (13) is irrelevant. Therefore, for $\left\lvert\delta\lambda_{i}\right\rvert<\left\lvert\lambda_{i}-\lambda_{j}\right\rvert/2$ the bound simplifies to

[TABLE]

Similar to the proof of Theorem 3.1, applying the Cauchy-Schwartz inequality we have that

[TABLE]

where in the last step we used Lemma 3.2. To finish the proof we notice that, due to Theorem 3.2, whenever $|\lambda_{i}-\lambda_{j}|\leq|\widetilde{\lambda}_{i}-\lambda_{j}|$, one has

[TABLE]

Our bound therefore holds for the union of intervals $\left\lvert\delta\lambda_{i}\right\rvert<\left\lvert\lambda_{i}-\lambda_{j}\right\rvert/2$ and $|\lambda_{i}-\lambda_{j}|\leq|\widetilde{\lambda}_{i}-\lambda_{j}|$, i.e., for $\widetilde{\lambda}_{i}>(\lambda_{i}+\lambda_{j})/2$ when $\lambda_{i}>\lambda_{j}$ and for $\widetilde{\lambda}_{i}<(\lambda_{i}+\lambda_{j})/2$ when $\lambda_{i}<\lambda_{j}$. ∎

4 Spectral concentration

This section builds on the perturbation results of Section 3 to characterize how far any inner product $\langle\widetilde{u}_{i},u_{j}\rangle$ and eigenvalue $\widetilde{\lambda}_{i}$ are from the ideal estimates.

Before proceeding, we remark on some simplifications employed in the following. W.l.o.g., we will assume that the mean $\mathbf{E}\left[x\right]$ is zero and the covariance full rank. Though the case of rank-deficient $C$ is easily handled by substituting the inverse with the Moore-Penrose pseudoinverse, we opt to make the exposition in the simpler setting. In addition, we will assume the perspective of Theorem 3.2, for which the inequality $\text{sgn}(\lambda_{i}>\lambda_{j})\,2\widetilde{\lambda}_{i}>\text{sgn}(\lambda_{i}>\lambda_{j})(\lambda_{i}+\lambda_{j})$ holds. This event occurs a.a.s. when the gap and the sample size are sufficiently large (see Section 4.1.2), but it is convenient to assume that it happens almost surely. Removing this assumption is possible, but is not pursued here as it leads to less elegant and sharp estimates.

4.1 Distributions with finite second moment

Our first flavor of results is based on a variant of the Tchebichef inequality and holds for any distribution with finite second moment.

4.1.1 Concentration of eigenvector angles

We start with the concentration of inner-products $|\langle\widetilde{u}_{i},u_{j}\rangle|$.

Theorem 4.1**.**

For any two eigenvectors $\widetilde{u}_{i}$ and $u_{j}$ of the sample and actual covariance respectively, with $\lambda_{i}\neq\lambda_{j}$, and for any real number $t>0$, we have

[TABLE]

for $\text{sgn}(\lambda_{i}>\lambda_{j})\,2\widetilde{\lambda}_{i}>\text{sgn}(\lambda_{i}>\lambda_{j})(\lambda_{i}+\lambda_{j})$ and $k_{j}=\left(\mathbf{E}\left[\|xx^{*}u_{j}\|_{2}^{2}\right]-\lambda_{j}^{2}\right)^{1/2}$.

Proof.

According to a variant of Tchebichef’s inequality [21], for any random variable $X$ and for any real numbers $t>0$ and $\alpha$:

[TABLE]

Setting $X=\langle\widetilde{u}_{i},u_{j}\rangle$ and $\alpha=0$, we have

[TABLE]

where the last inequality follows from Theorem 3.2. We continue by expanding $\delta C$ using the definition of the eigenvalue decomposition and substituting the expectation.

[TABLE]

In addition,

[TABLE]

and therefore

[TABLE]

Putting everything together, the claim follows. ∎

The following corollary will be very useful when applying our results.

Corollary 4.1**.**

For any weights $w_{ij}$ and real $t>0$:

[TABLE]

where $k_{j}=\left(\mathbf{E}\left[\|xx^{*}u_{j}\|_{2}^{2}\right]-\lambda_{j}^{2}\right)^{1/2}$ and $w_{ij}\neq 0$ when $\lambda_{i}\neq\lambda_{j}$ and $\text{sgn}(\lambda_{i}>\lambda_{j})\,2\widetilde{\lambda}_{i}>\text{sgn}(\lambda_{i}>\lambda_{j})(\lambda_{i}+\lambda_{j})$.

Proof.

We proceed as in the proof of Theorem 4.1:

[TABLE]

The claim follows by computing $\mathbf{E}\left[\|\delta Cu_{j}\|_{2}^{2}\right]$ (as before) and squaring both terms within the probability. ∎

4.1.2 Eigenvalue concentration

A slight modification of the same argument can be used to characterize the eigenvalue relative difference.

Theorem 4.2**.**

For any eigenvalues $\lambda_{i}$ and $\widetilde{\lambda}_{i}$ of $C$ and $\widetilde{C}$, respectively, and for any $t>0$, we have

[TABLE]

where $k_{i}=(\mathbf{E}\left[\|xx^{*}u_{i}\|_{2}^{2}\right]-\lambda_{i})^{1/2}$.

Proof.

Directly from the Bauker-Fike theorem [8] one sees that

[TABLE]

The proof is then identical to the one of Theorem 4.1. ∎

As such, the probability the main condition of our Theorems holds is at least

[TABLE]

4.1.3 The influence of the distribution

As seen by the straightforward inequality $\mathbf{E}\left[\|xx^{*}u_{j}\|_{2}^{2}\right]\leq\mathbf{E}\left[\|x\|_{2}^{4}\right]$, $k_{j}$ connects to the kurtosis of the distribution. However, it also captures the tendency of the distribution to fall in the span of $u_{j}$.

To see this, we will work with the uncorrelated random vectors $\varepsilon=C^{-1/2}x$, which have zero mean and identity covariance.

[TABLE]

If we further set $\hat{\varepsilon}=U^{*}\varepsilon$, we have

[TABLE]

It is therefore easier to untangle the spaces spanned by $\widetilde{u}_{i}$ and $u_{j}$ when the variance of the distribution along the latter space is small (the expression is trivially minimized when $\lambda_{j}\rightarrow 0$) or when the variance is entirely contained along that space (the expression is also small when $\lambda_{i}=0$ for all $i\neq j$). In addition, it can be seen that distributions with fast decaying tails allow for better principal component identification ($\mathbf{E}\left[\hat{\varepsilon}(j)^{4}\right]$ is a measure of kurtosis over the direction of $u_{j}$).

For the particular case of a Normal distribution, we provide a closed-form expression.

Corollary 4.2**.**

For a Normal distribution, we have $k_{j}^{2}=\lambda_{j}\left(\lambda_{j}+\text{tr}(C)\right)$.

Proof.

For a centered and normal distribution with identity covariance, the choice of basis is arbitrary and the vector $\hat{\varepsilon}=U^{*}\varepsilon$ is also zero mean with identity covariance. Moreover, for every $\ell\neq j$ we can write $\mathbf{E}\left[\hat{\varepsilon}(\ell)^{2}\hat{\varepsilon}(j)^{2}\right]=\mathbf{E}\left[\hat{\varepsilon}(\ell)^{2}\right]\mathbf{E}\left[\hat{\varepsilon}(j)^{2}\right]=1$. This implies that

[TABLE]

and, accordingly, $k_{j}^{2}=\lambda_{j}\left(\lambda_{j}+\text{tr}(C)\right)$. ∎

4.2 Distributions supported in a Euclidean ball

Our last result provides a sharper probability estimate for the family of sub-gaussian distributions supported in a centered Euclidean ball of radius $r$, with their $\Psi_{2}$-norm

[TABLE]

where $\mathcal{S}^{n-1}$ is the unit sphere and with the $\psi_{2}$-norm of a random variable $X$ defined as

[TABLE]

Our setting is therefore similar to the one used to study covariance estimation [26]. Due to space constraints, we refer the reader to the excellent review article [25] for an introduction to sub-gaussian distributions as a tool for non-asymptotic analysis of random matrices.

Theorem 4.3**.**

For sub-gaussian distributions supported within a centered Euclidean ball of radius $r$, there exists an absolute constant $c$, independent of the sample size, such that for any real number $t>0$,

[TABLE]

where $\Phi_{ij}(t)=\frac{|\lambda_{i}-\lambda_{j}|\,t-2\lambda_{j}}{2\,(r^{2}/\lambda_{j}-1)^{1/2}}-2\left\lVert x\right\rVert_{\Psi_{2}}$, $\lambda_{i}\neq\lambda_{j}$ and $\text{sgn}(\lambda_{i}>\lambda_{j})\,2\widetilde{\lambda}_{i}>\text{sgn}(\lambda_{i}>\lambda_{j})(\lambda_{i}+\lambda_{j})$.

Proof.

We start from the simple observation that, for every upper bound $B$ of $|\langle\widetilde{u}_{i},u_{j}\rangle|$ the relation ${\mathbf{P}\left(|\langle\widetilde{u}_{i},u_{j}\rangle|>t\right)}\leq{\mathbf{P}\left(B>t\right)}$ holds. To proceed therefore we will construct a bound with a known tail. As we saw in Sections 3.3 and 4.1,

[TABLE]

Assuming further that $\|x\|_{2}\leq r$, and since the numerator is minimized when $\|x_{p}\|_{2}^{2}$ approaches $\lambda_{j}$, we can write for every sample $x=C^{1/2}\varepsilon$:

[TABLE]

which is a shifted and scaled version of the random variable $|\hat{\varepsilon}(j)|=|u_{j}^{*}\varepsilon|$. Setting $a=(\lambda_{j}r^{2}-\lambda_{j}^{2})^{1/2}$, we have

[TABLE]

By Lemma 4.1 however, the left hand side is a sum of independent sub-gaussian variables. Since the summands are not centered, we expand each $|\hat{\varepsilon}_{p}(j)|=z_{p}+\mathbf{E}\left[|\hat{\varepsilon}_{p}(j)|\right]$ in terms of a centered sub-gaussian $z_{p}$ with the same $\psi_{2}$-norm. Furthermore, by Jensen’s inequality and Lemma 4.1

[TABLE]

Therefore, if we set $\Phi_{ij}(t)=\frac{(0.5\,|\lambda_{i}-\lambda_{j}|\,t-\lambda_{j})}{(r^{2}/\lambda_{j}-1)^{1/2}}-2\left\lVert x\right\rVert_{\Psi_{2}}$

[TABLE]

Moreover, by the rotation invariance principle, the left hand side of the last inequality is a sub-gaussian with $\psi_{2}$-norm smaller than $(C\sum_{p=1}^{m}\left\lVert z_{p}\right\rVert_{\psi_{2}}^{2})^{1/2}=(c_{1}m)^{1/2}\left\lVert z\right\rVert_{\psi_{2}}\leq(Cm/\lambda_{j})^{1/2}\left\lVert x\right\rVert_{\Psi_{2}}$, for some absolute constant $c_{1}$. As a consequence, there exists an absolute constant $c_{2}$, such that for each $\theta>0$:

[TABLE]

Substituting $\theta={m\,\Phi_{ij}(t)}/{\lambda_{j}}$, we have

[TABLE]

which is the desired bound. ∎

Lemma 4.1**.**

If $x$ is a sub-gaussian random vector and $\varepsilon=C^{-1/2}x$, then for every $i$, the random variable $\hat{\varepsilon}(i)=u_{i}^{*}\varepsilon$ is also sub-gaussian, with $\left\lVert\hat{\varepsilon}(i)\right\rVert_{\psi_{2}}\leq\left\lVert x\right\rVert_{\Psi_{2}}/\sqrt{\lambda_{i}}$.

Proof.

The fact that $\hat{\varepsilon}(i)$ is sub-gaussian follows easily by the definition of a sub-gaussian random vector, according to which for every $y\in\mathbb{R}^{n}$ the random variable $\langle x,y\rangle$ is sub-gaussian. Setting $y=u_{i}$, the first part is proven. For the bound of the norm, notice that

[TABLE]

where, for the last inequality, we set $y=u_{i}$. ∎

5 Principal component estimation

The last step in our exposition entails moving from statements about eigenvector angles to statements about the eigenvectors themselves. In particular, we are going to focus on the eigenvectors of the sample covariance associated to the largest eigenvalues, commonly referred to as principal components.

The phenomenon we will characterize is portrayed in Figure 2. The experiment in question concerns the $n=784$ dimensional distribution constructed by the 6131 images featuring digit ‘3’ found in the MNIST database. The top sub-figure shows the estimation accuracy for a set of five principal components averaged over 200 sampling draws. The trends suggest that the first three eigenvectors can be estimated up to a satisfactory accuracy by much less than $n$ samples. The sample requirements decrease further if we are satisfied by a tight localization of eigenvectors222This is relevant for instance when we wish to construct projectors over specific principal eigenspaces and we have to ensure that the projection space estimated from the sample covariance closely approximates an eigenspace of the actual covariance., such that $\widetilde{u}_{i}$ is almost entirely contained within the span of the $2\ell$ surrounding eigenvectors of $u_{i}$. As suggested by the middle and bottom sub-figures, by slightly increasing $\ell$, we reduce the number of samples needed to estimate higher principal components.

In the following, we show that our results verify these trends. We are interested to find out how many samples are sufficient to ensure that $\widetilde{u}_{i}$ is almost entirely contained within $\text{span}(u_{i-\ell},\ldots,u_{i},\ldots,u_{i+\ell})$ for some small non-negative integer $\ell$. Setting $\mathcal{S}=\{(j<i-\ell)\cup(i+\ell<j\leq r)\}$, where $r$ is the rank of $C$, we have as a consequence of Corollary 4.1 that for distributions with finite second moment:

[TABLE]

In accordance with the experimental results, equation (38) reveals that it is much easier to estimate the principal components of larger variance, and that, by introducing a small slack in terms of $\ell$ one mitigates the requirement for eigenvalue separation.

It might be also interesting to observe that, for covariance matrices that are (approximately) low-rank, we obtain estimates reminiscent of compressed setting [10], in the sense that the sample requirement becomes a function of the non-zero eigenvalues. Though intuitive, this dependency of the estimation accuracy on the rank was not transparent in known results for covariance estimation [20, 1, 26].

6 Conclusions

The main contribution of this paper was the derivation of non-asymptotic bounds for the concentration of inner-products $|\langle\widetilde{u}_{i},u_{j}\rangle|$ involving eigenvectors of the sample and actual covariance matrices. We also showed how these results can be extended to reason about eigenvectors, eigenspaces, and eigenvalues.

We have identified two interesting directions for further research. The first has to do with obtaining tighter concentration estimates. Especially with regards to our perturbation arguments, we believe that our current bounds on inner products could be sharpened by at least a constant multiplicative factor. We also suspect that a joint analysis of angles could also lead to a significant improvement over Corollary 4.1. The second direction involves using our results for the analysis of methods that utilize the eigenvectors of the covariance for dimensionality reduction. Examples include (fast) principal component projection [12] and regression [17].