Lecture 9: Covariance Estimation
Two-sided singular value bounds and the sample covariance
Lecture 8 showed that a tall random matrix is an approximate isometry: for \(m\geq C_0n\), every vector is distorted by a factor \(1\pm\eta\). The proof was a net and a union bound, and it left the dependence of the distortion on \(m\) and \(n\) unresolved. This lecture sharpens the argument to the form \[ \sqrt m-C\sqrt n\ \leq\ s_{\min}(A)\ \leq\ s_{\max}(A)\ \leq\ \sqrt m+C\sqrt n, \] promised at the time, and relaxes the hypothesis from independent entries to independent rows. The relaxation is what makes the result statistical. If the rows are samples from a distribution, then \(\frac1mA^{\mathsf T}A\) is the sample covariance, and the theorem says how many samples determine a covariance matrix.
1 The covariance estimation problem
Let \(X\) be a mean-zero random vector in \(\mathbb R^n\) with \(\mathbb E\|X\|_2^2<\infty\), and let \[ \Sigma:=\mathbb EXX^{\mathsf T} \] be its covariance matrix, an \(n\times n\) positive semidefinite matrix with \(\langle\Sigma x,x\rangle=\mathbb E\langle X,x\rangle^2\). Given independent samples \(X_1,\ldots,X_m\) distributed as \(X\), the natural estimate of \(\Sigma\) is the sample covariance \[ \Sigma_m:=\frac1m\sum_{i=1}^mX_iX_i^{\mathsf T}. \] It is unbiased, \(\mathbb E\Sigma_m=\Sigma\), and the law of large numbers applied to each of the \(n^2\) entries gives \(\Sigma_m\to\Sigma\) almost surely as \(m\to\infty\) with \(n\) fixed. The high-dimensional question is quantitative: how large must \(m\) be, as a function of \(n\), so that \(\Sigma_m\approx\Sigma\) in operator norm with high probability?
Some dependence on \(n\) is unavoidable. The matrix \(\Sigma_m\) is a sum of \(m\) rank-one matrices, so it has rank at most \(m\). If \(\Sigma\) is invertible and \(m<n\), then \(\Sigma_m\) has a nontrivial kernel while \(\Sigma\) does not, and \[ \|\Sigma_m-\Sigma\|_{\mathrm{op}}\geq\lambda_{\min}(\Sigma)>0 . \] So at least \(m\geq n\) samples are needed. The theorem of this lecture says that \(m\) proportional to \(n\) suffices.
The operator norm is the right way to measure the error because it controls what the covariance is used for. The quadratic form \(x\mapsto\langle\Sigma x,x\rangle\) is the variance of the data in direction \(x\), and \(\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq\varepsilon\) says that every directional variance is estimated to within \(\varepsilon\). Principal component analysis reads off the leading eigenvectors of \(\Sigma\), and by the Davis–Kahan theorem the leading eigenvectors of \(\Sigma_m\) are close to those of \(\Sigma\) whenever the operator norm error is small compared to the gap between consecutive eigenvalues of \(\Sigma\). An entrywise error bound would not give this: \(n^2\) small entrywise errors can add up to an operator norm error of size \(n\) times larger.
2 Sub-Gaussian random vectors
The hypothesis on the data is the vector version of the sub-Gaussian assumption of Lecture 4. Recall that a mean-zero random variable is \(V\)-subgaussian if \(\mathbb Ee^{\lambda\xi}\leq e^{\lambda^2V/2}\) for every real \(\lambda\).
Definition 1 (Sub-Gaussian and isotropic random vectors) A mean-zero random vector \(X\) in \(\mathbb R^n\) is \(K\)-subgaussian if for every \(x\in\mathbb R^n\) the random variable \(\langle X,x\rangle\) is \(K\|x\|_2^2\)-subgaussian. It is isotropic if \(\mathbb EXX^{\mathsf T}=I_n\), that is, if \(\mathbb E\langle X,x\rangle^2=\|x\|_2^2\) for every \(x\).
A \(K\)-subgaussian vector has every one-dimensional marginal \(\langle X,x\rangle\) sub-Gaussian at the same scale, in proportion to the length of \(x\). Isotropy fixes that scale: the variance of a mean-zero \(V\)-subgaussian variable is at most \(V\) (compare the Taylor expansions at \(\lambda=0\) of the two sides of the defining inequality), so an isotropic \(K\)-subgaussian vector always has \(K\geq1\). We use this normalization throughout.
Gaussian vectors. If \(X\sim N(0,I_n)\), then \(\langle X,x\rangle\sim N(0,\|x\|_2^2)\), which is \(\|x\|_2^2\)-subgaussian. So \(X\) is isotropic and \(1\)-subgaussian.
Independent coordinates. If the coordinates of \(X\) are independent, mean-zero, of variance \(1\), and \(K\)-subgaussian, then \(X\) is isotropic, and by the Pythagorean theorem of Lecture 5, \(\langle X,x\rangle=\sum_ix_iX_i\) is \(\big(\sum_ix_i^2K\big)\)-subgaussian. So \(X\) is \(K\)-subgaussian. Random sign vectors are the case \(K=1\). This is exactly the row structure assumed in Lecture 8.
Random directions. If \(X\) is uniform on the sphere of radius \(\sqrt n\), then \(X\) is isotropic by symmetry, and it is \(C\)-subgaussian for an absolute constant \(C\). We do not prove this now; it follows from the concentration of measure on the sphere, which is the subject of Lecture 11.
Anisotropic data. If \(Z\) is isotropic and \(K\)-subgaussian and \(\Sigma\) is positive semidefinite, then \(X=\Sigma^{1/2}Z\) has covariance \(\Sigma\) and \(\langle X,x\rangle=\langle Z,\Sigma^{1/2}x\rangle\) is \(K\|\Sigma^{1/2}x\|_2^2=K\langle\Sigma x,x\rangle\)-subgaussian. So the variance proxy of each marginal is \(K\) times its true variance. This is the form in which the covariance theorem is stated.
The examples suggest the natural hypothesis for covariance estimation: the variance proxy of each marginal should be comparable to its true variance. The next lemma shows that under this hypothesis, the data can be standardized.
Lemma 1 (Reduction to the isotropic case) Let \(X\) be a mean-zero random vector with invertible covariance \(\Sigma\), and suppose that for every \(x\in\mathbb R^n\) the marginal \(\langle X,x\rangle\) is \(K\langle\Sigma x,x\rangle\)-subgaussian. Then \(Z:=\Sigma^{-1/2}X\) is isotropic and \(K\)-subgaussian, and \(X=\Sigma^{1/2}Z\).
Here \(\Sigma^{1/2}\) is the positive definite square root of \(\Sigma\) and \(\Sigma^{-1/2}\) its inverse. For every \(x\), \[ \langle Z,x\rangle=\langle X,\Sigma^{-1/2}x\rangle, \qquad \mathbb E\langle Z,x\rangle^2=\langle\Sigma\Sigma^{-1/2}x,\Sigma^{-1/2}x\rangle=\|x\|_2^2 , \] so \(Z\) is isotropic. By hypothesis \(\langle Z,x\rangle\) is \(K\langle\Sigma\Sigma^{-1/2}x,\Sigma^{-1/2}x\rangle=K\|x\|_2^2\)-subgaussian.
If \(\Sigma\) is singular, \(X\) lies almost surely in the range of \(\Sigma\), a subspace of dimension \(r=\operatorname{rank}\Sigma<n\), and everything below goes through on that subspace with \(n\) replaced by \(r\). We assume \(\Sigma\) invertible to keep the notation light.
3 Nets and quadratic forms
The two facts about nets from Lecture 8 are that \(S^{n-1}\) has an \(\varepsilon\)-net of cardinality at most \((3/\varepsilon)^n\), and that a bound on \(\|Ax\|_2\) over a net transfers to the sphere. For the sharper result we need the transfer for quadratic forms. This is Problem 3 of Exercise Sheet 4, and we record only the statement.
For a symmetric matrix \(B\), the spectral theorem gives \[ \|B\|_{\mathrm{op}}=\sup_{\|x\|_2=1}|\langle Bx,x\rangle| . \]
Lemma 2 (Quadratic forms on a net) Let \(B\) be a symmetric \(n\times n\) matrix, let \(0<\varepsilon<\frac12\), and let \(N\) be an \(\varepsilon\)-net of \(S^{n-1}\). Then \[ \|B\|_{\mathrm{op}}\leq\frac1{1-2\varepsilon}\max_{y\in N}|\langle By,y\rangle| . \] In particular, for \(\varepsilon=\frac14\), \(\|B\|_{\mathrm{op}}\leq2\max_{y\in N}|\langle By,y\rangle|\).
The second fact converts a bound on \(\frac1mA^{\mathsf T}A-I\) into a bound on the singular values of \(A\). It is the same computation as in the proof of the one-vector lemma of Lecture 8, which passed from \(\|Ax\|_2^2\) to \(\|Ax\|_2\).
Lemma 3 (From the Gram matrix to the singular values) Let \(A\) be an \(m\times n\) matrix and \(\delta\geq0\). If \[ \Big\|\frac1mA^{\mathsf T}A-I_n\Big\|_{\mathrm{op}}\leq\max\{\delta,\delta^2\}, \] then \[ (1-\delta)\sqrt m\leq s_{\min}(A)\leq s_{\max}(A)\leq(1+\delta)\sqrt m . \]
For a unit vector \(x\), the hypothesis says \(|\|Ax\|_2^2/m-1|\leq\max\{\delta,\delta^2\}\). Write \(a=\|Ax\|_2/\sqrt m\geq0\) and suppose \(|a-1|>\delta\). If \(a>1+\delta\) then \(a^2-1>2\delta+\delta^2\geq\max\{\delta,\delta^2\}\). If \(a<1-\delta\), which forces \(\delta<1\), then \(1-a^2=(1-a)(1+a)>\delta=\max\{\delta,\delta^2\}\). Both contradict the hypothesis, so \(|a-1|\leq\delta\) for every unit \(x\).
4 Two-sided bounds for matrices with independent rows
Theorem 1 (Two-sided bound, independent sub-Gaussian rows) Let \(A\) be an \(m\times n\) matrix whose rows \(A_1,\ldots,A_m\) are independent, mean-zero, isotropic, \(K\)-subgaussian random vectors in \(\mathbb R^n\). There is an absolute constant \(C\) such that, for every \(t\geq0\), with probability at least \(1-2e^{-t^2}\), \[ \Big\|\frac1mA^{\mathsf T}A-I_n\Big\|_{\mathrm{op}}\leq K\max\{\delta,\delta^2\}, \qquad \delta:=C\Big(\sqrt{\frac nm}+\frac t{\sqrt m}\Big). \] Consequently, on the same event, \[ \sqrt m-CK(\sqrt n+t)\ \leq\ s_{\min}(A)\ \leq\ s_{\max}(A)\ \leq\ \sqrt m+CK(\sqrt n+t). \]
The rows need not have independent coordinates, and they need not be identically distributed. What is used is that each row is an isotropic sub-Gaussian vector and that the rows are independent of one another. When the rows are samples from a distribution, this is exactly the situation of the covariance estimation problem.
Write \(B:=\frac1mA^{\mathsf T}A-I_n\), a symmetric matrix, and note that for every \(x\in\mathbb R^n\), \[ \langle Bx,x\rangle=\frac1m\|Ax\|_2^2-\|x\|_2^2 =\frac1m\sum_{i=1}^m\big(\langle A_i,x\rangle^2-\|x\|_2^2\big). \] Set \(\varepsilon:=K\max\{\delta,\delta^2\}\).
Step 1: the net. Let \(N\) be a \(\frac14\)-net of \(S^{n-1}\) with \(|N|\leq12^n\), which exists by the net lemma of Lecture 8. By Lemma 2, \(\|B\|_{\mathrm{op}}\leq2\max_{y\in N}|\langle By,y\rangle|\), so it suffices to show that \[ \max_{y\in N}|\langle By,y\rangle|\leq\frac\varepsilon2 \] with probability at least \(1-2e^{-t^2}\).
Step 2: one vector. Fix \(y\in S^{n-1}\) and set \(\xi_i:=\langle A_i,y\rangle\). The \(\xi_i\) are independent, mean-zero, with \(\mathbb E\xi_i^2=1\), and \(K\)-subgaussian. By the sub-Gaussian–\(\psi_2\) comparison of Lecture 5, \(\|\xi_i\|_{\psi_2}\leq2\sqrt K\), and by the Orlicz–Hölder identity of the same lecture, \(\|\xi_i^2\|_{\psi_1}=\|\xi_i\|_{\psi_2}^2\leq4K\). The constant \(1\) has \(\|1\|_{\psi_1}=1/\log2<2\), so by the triangle inequality and \(K\geq1\), \[ \|\xi_i^2-1\|_{\psi_1}\leq4K+2\leq6K . \] Now \(\langle By,y\rangle=\frac1m\sum_i(\xi_i^2-1)\) is an average of \(m\) independent, mean-zero, sub-exponential variables. Bernstein’s inequality (Lecture 6, \(\psi_1\) version) gives \[ \mathbb P\Big\{|\langle By,y\rangle|>\frac\varepsilon2\Big\} =\mathbb P\Big\{\Big|\sum_{i=1}^m(\xi_i^2-1)\Big|>\frac{m\varepsilon}2\Big\} \leq2\exp\Big(-c\,m\min\Big\{\frac\varepsilon{12K},\frac{\varepsilon^2}{144K^2}\Big\}\Big). \] Since \(\varepsilon/K=\max\{\delta,\delta^2\}\geq\delta^2\) and \(\varepsilon^2/K^2=\max\{\delta,\delta^2\}^2\geq\delta^2\), the minimum is at least \(\delta^2/144\), and \[ m\delta^2=C^2(\sqrt n+t)^2\geq C^2(n+t^2). \] So the one-vector failure probability is at most \(2\exp\big(-c_1C^2(n+t^2)\big)\) with \(c_1=c/144\).
Step 3: the union bound. Over the net, \[ \mathbb P\Big\{\max_{y\in N}|\langle By,y\rangle|>\frac\varepsilon2\Big\} \leq12^n\cdot2\exp\big(-c_1C^2(n+t^2)\big) =2\exp\big(n\log12-c_1C^2n-c_1C^2t^2\big) \leq2e^{-t^2}, \] provided \(C\) is chosen so that \(c_1C^2\geq\log12\), which also gives \(c_1C^2\geq1\). This proves the first claim.
Step 4: singular values. Since \(K\geq1\), \(K\max\{\delta,\delta^2\}\leq\max\{K\delta,(K\delta)^2\}\), so Lemma 3 applies with \(K\delta\) in place of \(\delta\) and gives \((1-K\delta)\sqrt m\leq s_{\min}(A)\leq s_{\max}(A)\leq(1+K\delta)\sqrt m\). Finally \(K\delta\sqrt m=CK(\sqrt n+t)\).
Compare this with Lecture 8. There, entries were independent and the conclusion was a distortion \(\eta\) that did not improve as \(m/n\) grew. Here the distortion is \(CK\sqrt{n/m}\), which tends to zero. The price is that the theorem says nothing useful about \(s_{\min}(A)\) once \(CK\sqrt{n/m}\geq1\), that is, once \(m\) is only a constant multiple of \(n\). In that regime the union bound is again too weak, and the small-ball estimates of Lecture 8 are still the only tool we have. For Gaussian matrices the constant \(CK\) can be replaced by \(1\), and Lecture 12 proves this by an entirely different method.
The theorem also yields a bound in expectation, which is the form used most often.
Corollary 1 (Expected deviation of the Gram matrix) In the setting of Theorem 1, \[ \mathbb E\Big\|\frac1mA^{\mathsf T}A-I_n\Big\|_{\mathrm{op}}\leq CK\Big(\sqrt{\frac nm}+\frac nm\Big). \]
Let \(Y=\|\frac1mA^{\mathsf T}A-I_n\|_{\mathrm{op}}\) and, for \(t\geq0\), let \[ \delta(t)=C\Big(\sqrt{\frac nm}+\frac t{\sqrt m}\Big), \qquad s(t)=K\big(\delta(t)+\delta(t)^2\big)\geq K\max\{\delta(t),\delta(t)^2\}. \] The function \(s\) is increasing, and Theorem 1 gives \(\mathbb P\{Y>s(t)\}\leq2e^{-t^2}\). Writing the expectation as a tail integral and substituting \(u=s(t)\), \[ \mathbb EY=\int_0^\infty\mathbb P\{Y>u\}\,du \leq s(0)+\int_0^\infty\mathbb P\{Y>s(t)\}\,s'(t)\,dt \leq s(0)+2\int_0^\infty e^{-t^2}s'(t)\,dt . \] Here \(s'(t)=K\delta'(t)(1+2\delta(t))\) with \(\delta'(t)=C/\sqrt m\), so \[ \int_0^\infty e^{-t^2}s'(t)\,dt \leq\frac{CK}{\sqrt m}\int_0^\infty e^{-t^2}\Big(1+2C\sqrt{\frac nm}+\frac{2Ct}{\sqrt m}\Big)dt \leq\frac{C'K}{\sqrt m}\Big(1+\sqrt{\frac nm}+\frac1{\sqrt m}\Big). \] Since \(n\geq1\), each of \(1/\sqrt m\), \(\sqrt n/m\), and \(1/m\) is at most \(\sqrt{n/m}+n/m\). Similarly \(s(0)=K(\delta(0)+\delta(0)^2)=K\big(C\sqrt{n/m}+C^2n/m\big)\). Adding up gives the claim with a new absolute constant.
5 Covariance estimation
Theorem 2 (Covariance estimation) Let \(X\) be a mean-zero random vector in \(\mathbb R^n\) with invertible covariance \(\Sigma\), and suppose that for every \(x\in\mathbb R^n\) the marginal \(\langle X,x\rangle\) is \(K\langle\Sigma x,x\rangle\)-subgaussian. Let \(X_1,\ldots,X_m\) be independent copies of \(X\) and \(\Sigma_m\) their sample covariance. Then \[ \mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq CK\Big(\sqrt{\frac nm}+\frac nm\Big)\|\Sigma\|_{\mathrm{op}}, \] and for every \(u\geq0\), with probability at least \(1-2e^{-u}\), \[ \|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq CK\Big(\sqrt{\frac{n+u}m}+\frac{n+u}m\Big)\|\Sigma\|_{\mathrm{op}} . \]
By Lemma 1, \(Z_i:=\Sigma^{-1/2}X_i\) are independent, isotropic, \(K\)-subgaussian, and \(X_i=\Sigma^{1/2}Z_i\). Let \(A\) be the \(m\times n\) matrix with rows \(Z_i^{\mathsf T}\). Then \[ \Sigma_m-\Sigma =\frac1m\sum_{i=1}^m\Sigma^{1/2}Z_iZ_i^{\mathsf T}\Sigma^{1/2}-\Sigma^{1/2}\Sigma^{1/2} =\Sigma^{1/2}\Big(\frac1mA^{\mathsf T}A-I_n\Big)\Sigma^{1/2}, \] and since \(\|\Sigma^{1/2}\|_{\mathrm{op}}^2=\|\Sigma\|_{\mathrm{op}}\), \[ \|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq\|\Sigma\|_{\mathrm{op}}\Big\|\frac1mA^{\mathsf T}A-I_n\Big\|_{\mathrm{op}} . \] The expectation bound is now Corollary 1. For the tail bound, apply Theorem 1 with \(t=\sqrt u\). Then \(\delta=C(\sqrt{n/m}+\sqrt{u/m})\leq2C\sqrt{(n+u)/m}\), and \(K\max\{\delta,\delta^2\}\leq K(\delta+\delta^2)\) gives the claim.
Fix \(\varepsilon\in(0,1)\). The theorem gives relative error \(\mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq\varepsilon\|\Sigma\|_{\mathrm{op}}\) as soon as \[ m\geq\frac{C'K^2n}{\varepsilon^2}. \] So the sample covariance is an accurate estimate of the covariance once the number of samples is a constant multiple of the dimension, and the constant scales like \(\varepsilon^{-2}\), the usual price of a relative error \(\varepsilon\) in a variance estimate. This matches the necessary condition \(m\geq n\) up to the constant.
If \(X\) has an unknown mean \(\mu\), the theorem applies verbatim to the second moment matrix \(\mathbb EXX^{\mathsf T}=\Sigma+\mu\mu^{\mathsf T}\), provided the sub-Gaussian hypothesis is placed on \(\langle X-\mu,x\rangle\) and the mean is not too large relative to \(\Sigma\). Estimating \(\Sigma\) itself requires subtracting the sample mean, which adds an error of order \(\|\bar X_m-\mu\|_2^2\), and this is of smaller order than the main term. We omit the bookkeeping.
6 What the rate cannot see
The bound \(\sqrt{n/m}\) depends on the ambient dimension and on nothing else about \(\Sigma\). Two remarks put this in perspective.
The rate is sharp for isotropic data. If \(X\sim N(0,I_n)\) then \(\Sigma_m=\frac1mA^{\mathsf T}A\) for a Gaussian matrix \(A\), and its extreme eigenvalues are \(s_{\min}(A)^2/m\) and \(s_{\max}(A)^2/m\). Lecture 12 shows that for Gaussian matrices \[ \mathbb Es_{\max}(A)\leq\sqrt m+\sqrt n, \qquad \mathbb Es_{\min}(A)\geq\sqrt m-\sqrt n, \] and both bounds are attained to leading order. The eigenvalues of \(\Sigma_m\) therefore spread over an interval of half-width about \(2\sqrt{n/m}\) around \(1\), and \(\|\Sigma_m-\Sigma\|_{\mathrm{op}}\approx2\sqrt{n/m}\). In fact the whole empirical spectrum of \(\Sigma_m\) has a deterministic limiting shape, the Marchenko–Pastur law (Marchenko and Pastur 1967), supported on exactly this interval.
But the dimension is the wrong parameter for structured data. Suppose the data are nearly low-dimensional: \(\Sigma\) has a few large eigenvalues and many small ones. The theorem still charges for all \(n\) directions, even though the small directions contribute almost nothing to \(\|\Sigma_m-\Sigma\|_{\mathrm{op}}\). The correct parameter is the effective rank \[ r(\Sigma):=\frac{\operatorname{tr}\Sigma}{\|\Sigma\|_{\mathrm{op}}}\in[1,n], \] which equals \(n\) for isotropic data and stays bounded when the spectrum decays fast. For Gaussian data one has \[ \mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq C\Big(\sqrt{\frac{r(\Sigma)}m}+\frac{r(\Sigma)}m\Big)\|\Sigma\|_{\mathrm{op}}, \] a theorem of Koltchinskii and Lounici (Koltchinskii and Lounici 2017). The net argument of this lecture cannot prove it. The proof of Theorem 1 pays \(n\log12\) for the net regardless of the geometry of \(\Sigma\), and no choice of net on the ellipsoid \(\Sigma^{1/2}S^{n-1}\) repairs this. What is needed is a way to control \(\sup_x|\frac1m\|Ax\|_2^2-\|x\|_2^2|\) over the ellipsoid that sees its shape and not only its dimension. That is what the theory of Gaussian processes provides, and Lecture 20 proves the effective-rank theorem with it.
7 Where we go next
Two loose ends of this lecture point in the same direction. The constant in \(\sqrt m\pm CK\sqrt n\) should be \(1\) for Gaussian matrices, and the dimension \(n\) in the covariance bound should be the effective rank. Both statements are about the random function \(x\mapsto\|Ax\|_2\) on a set of vectors, and both are proved by comparing it to simpler Gaussian processes rather than by discretizing the set. The next lectures develop Gaussian processes from the beginning: their construction, their concentration, and the comparison inequalities that make them computable.