Lecture 20: Covariance Estimation for Low-Dimensional Data
Effective rank, ellipsoids, and the matrix deviation inequality at work
Lecture 9 proved that \(m\gtrsim n\) samples estimate an \(n\)-dimensional covariance matrix in operator norm, and observed that the bound charges for all \(n\) directions even when the data are concentrated near a low-dimensional subspace. The net argument could not do better, because the net of the sphere has \(e^{cn}\) points whatever \(\Sigma\) is. The matrix deviation inequality of Lecture 19 replaces the net by the Gaussian width of the relevant index set, and the relevant set is an ellipsoid whose width is \(\sqrt{\operatorname{tr}\Sigma}\) rather than \(\sqrt{n\|\Sigma\|_{\mathrm{op}}}\). The result is the theorem stated without proof in Lecture 9: the sample size that matters is the effective rank \(\operatorname{tr}\Sigma/\|\Sigma\|_{\mathrm{op}}\), not the dimension.
1 The deviation of squares
The covariance error is a quadratic form in the data, so the first step is to pass from \(\|Ax\|_2-\sqrt m\|x\|_2\) to its square. Throughout, \(A\) is an \(m\times n\) matrix with independent, mean-zero, isotropic, \(K\)-subgaussian rows, as in Lectures 9 and 19, and for bounded \(T\subseteq\mathbb R^n\) we write \(w(T)\) for the Gaussian width and \(\operatorname{rad}(T)=\sup_{x\in T}\|x\|_2\).
Proposition 1 (Deviation of squares) For every bounded \(T\subseteq\mathbb R^n\), \[ \mathbb E\sup_{x\in T}\Big|\|Ax\|_2^2-m\|x\|_2^2\Big| \leq CK^2\big(w(T)+\operatorname{rad}(T)\big)^2+CK\sqrt m\operatorname{rad}(T)\big(w(T)+\operatorname{rad}(T)\big). \]
Write \(X_x=\|Ax\|_2-\sqrt m\|x\|_2\) and \(S=\sup_{x\in T}|X_x|\). Since \(a^2-b^2=(a-b)(a+b)\) and \(\|Ax\|_2+\sqrt m\|x\|_2\leq|X_x|+2\sqrt m\|x\|_2\), \[ \Big|\|Ax\|_2^2-m\|x\|_2^2\Big|=|X_x|\big(\|Ax\|_2+\sqrt m\|x\|_2\big) \leq S^2+2\sqrt m\operatorname{rad}(T)\,S \] for every \(x\in T\). By the tail form of the matrix deviation inequality (Lecture 19), \(S\leq CK(1+u)(w(T)+\operatorname{rad}(T))\) with probability at least \(1-2e^{-u^2}\), so \(S\) has a sub-Gaussian tail at scale \(CK(w(T)+\operatorname{rad}(T))\), and integrating gives \(\mathbb ES^2\leq CK^2(w(T)+\operatorname{rad}(T))^2\) and \(\mathbb ES\leq CK(w(T)+\operatorname{rad}(T))\). Take expectations.
The two terms have different origins. The first, of order \(K^2w(T)^2\), is the square of the deviation of the norm process and does not grow with \(m\). The second, of order \(\sqrt m\cdot Kw(T)\), is the cross term between the deviation and the typical size \(\sqrt m\|x\|_2\) of \(\|Ax\|_2\). After division by \(m\), they become the \(r/m\) and \(\sqrt{r/m}\) terms of the covariance bound.
2 Covariance estimation with effective rank
Definition 1 (Effective rank) The effective rank of a nonzero positive semidefinite matrix \(\Sigma\) is \[ r(\Sigma):=\frac{\operatorname{tr}\Sigma}{\|\Sigma\|_{\mathrm{op}}}=\frac{\sum_i\lambda_i(\Sigma)}{\max_i\lambda_i(\Sigma)}\in[1,n]. \]
The effective rank is \(n\) when all eigenvalues are equal and stays bounded when the eigenvalues decay quickly. It is the stable rank of \(\Sigma^{1/2}\) in the sense of Lecture 13, since \(\|\Sigma^{1/2}\|_F^2=\operatorname{tr}\Sigma\) and \(\|\Sigma^{1/2}\|_{\mathrm{op}}^2=\|\Sigma\|_{\mathrm{op}}\), and up to a constant it is the stable dimension of the ellipsoid \(\Sigma^{1/2}B_2^n\), whose width is \(\mathbb E\|\Sigma^{1/2}g\|_2\leq\sqrt{\operatorname{tr}\Sigma}\) and whose radius is \(\|\Sigma\|_{\mathrm{op}}^{1/2}\). Unlike the rank, it is continuous in \(\Sigma\).
Theorem 1 (Covariance estimation for low-dimensional data) Let \(X\) be a mean-zero random vector in \(\mathbb R^n\) with 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 \(\Sigma_m\) be the sample covariance of \(m\) independent copies of \(X\). Then \[ \mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq CK^2\Big(\sqrt{\frac{r(\Sigma)}m}+\frac{r(\Sigma)}m\Big)\|\Sigma\|_{\mathrm{op}} . \]
As in Lecture 9, write \(X_i=\Sigma^{1/2}Z_i\) with \(Z_i\) independent, isotropic, \(K\)-subgaussian, and let \(A\) be the matrix with rows \(Z_i^{\mathsf T}\); if \(\Sigma\) is singular, work in its range. Then \(\Sigma_m-\Sigma=\Sigma^{1/2}(\frac1mA^{\mathsf T}A-I)\Sigma^{1/2}\), and since this matrix is symmetric its operator norm is the supremum of the absolute value of its quadratic form over the sphere: \[ \|\Sigma_m-\Sigma\|_{\mathrm{op}} =\sup_{\|x\|_2=1}\Big|\Big\langle\Big(\frac1mA^{\mathsf T}A-I\Big)\Sigma^{1/2}x,\Sigma^{1/2}x\Big\rangle\Big| =\sup_{y\in\mathcal E}\Big|\frac1m\|Ay\|_2^2-\|y\|_2^2\Big|, \] where \(\mathcal E:=\Sigma^{1/2}S^{n-1}\) is an ellipsoid. This is the whole point: the net argument of Lecture 9 evaluated the same quadratic form on the sphere, where every direction costs the same, whereas here it is evaluated on \(\mathcal E\), whose extent in each direction is the standard deviation of the data in that direction.
Now apply Proposition 1 to \(T=\mathcal E\). Its width is \[ w(\mathcal E)=\mathbb E\sup_{\|x\|_2=1}\langle g,\Sigma^{1/2}x\rangle=\mathbb E\|\Sigma^{1/2}g\|_2 \leq\big(\mathbb E\|\Sigma^{1/2}g\|_2^2\big)^{1/2}=\sqrt{\operatorname{tr}\Sigma}, \] and its radius is \(\|\Sigma^{1/2}\|_{\mathrm{op}}=\|\Sigma\|_{\mathrm{op}}^{1/2}\). Writing \(\sigma=\|\Sigma\|_{\mathrm{op}}\) and \(r=r(\Sigma)\), so that \(\operatorname{tr}\Sigma=r\sigma\) and \(\sqrt{r\sigma}+\sqrt\sigma\leq2\sqrt{r\sigma}\) because \(r\geq1\), \[ \mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}} \leq\frac1m\Big[CK^2\cdot4r\sigma+CK\sqrt m\sqrt\sigma\cdot2\sqrt{r\sigma}\Big] =C'\sigma\Big[K^2\frac rm+K\sqrt{\frac rm}\Big], \] and \(K\leq K^2\) since \(K\geq1\).
Compare with Lecture 9, whose bound was the same with \(n\) in place of \(r(\Sigma)\). Since \(r(\Sigma)\leq n\) always, the new theorem is never worse, and it is much better when the spectrum of \(\Sigma\) decays.
Relative error. For \(\varepsilon\in(0,1)\), \(\mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq\varepsilon\|\Sigma\|_{\mathrm{op}}\) as soon as \(m\geq CK^4r(\Sigma)/\varepsilon^2\). The ambient dimension has disappeared from the sample complexity.
Fast spectral decay. If \(\lambda_i(\Sigma)\leq\lambda_1i^{-2}\), then \(r(\Sigma)\leq\pi^2/6\), and a number of samples independent of \(n\) estimates \(\Sigma\) to any fixed relative accuracy, in dimension \(n=10^6\) as well as in dimension \(10\). The estimate is in operator norm, so it controls the leading eigenvalues and eigenvectors of \(\Sigma\) through the Davis–Kahan theorem, which is what principal component analysis needs.
Spiked covariance. If \(\Sigma=I_n+\theta vv^{\mathsf T}\) for a unit vector \(v\) and \(\theta>0\), then \(r(\Sigma)=(n+\theta)/(1+\theta)\), which is about \(n/\theta\) for large \(\theta\). A strong spike reduces the effective rank, and the theorem says the spike direction is learned from \(m\gtrsim n/\theta\) samples, up to the accuracy needed for Davis–Kahan.
Two regimes. For \(m\geq r\) the bound is \(\sqrt{r/m}\), the parametric rate for \(r\) effective parameters. For \(m<r\) it is \(r/m>1\), which is worse than the trivial bound \(\|\Sigma_m-\Sigma\|_{\mathrm{op}}\leq\|\Sigma_m\|_{\mathrm{op}}+\|\Sigma\|_{\mathrm{op}}\) only by a constant, and reflects that \(\Sigma_m\) has rank at most \(m\).
For Gaussian data the theorem is sharp in both terms: Koltchinskii and Lounici (Koltchinskii and Lounici 2017) proved matching lower bounds, so that \(\mathbb E\|\Sigma_m-\Sigma\|_{\mathrm{op}}\asymp(\sqrt{r/m}+r/m)\|\Sigma\|_{\mathrm{op}}\) for \(X\sim N(0,\Sigma)\), and they showed that \(\|\Sigma_m-\Sigma\|_{\mathrm{op}}\) concentrates around its mean at the scale \(\|\Sigma\|_{\mathrm{op}}(\sqrt{1/m}+1/m)\), much smaller than the mean itself when \(r\) is large. A high-probability version of the theorem above, with \(r+u^2\) in place of \(r\) on an event of probability \(1-2e^{-u^2}\), follows by running the proof with the tail form of the matrix deviation inequality instead of its expectation.
3 Random projections, revisited
The same inequality sharpens the random projection theorem of Lecture 13.
Proposition 2 (Sizes of random projections, sharp constant) Let \(A\) have independent, mean-zero, isotropic, \(K\)-subgaussian rows and let \(T\subseteq\mathbb R^n\) be bounded. Then \[ \mathbb E\operatorname{diam}(AT)\leq\sqrt m\operatorname{diam}(T)+CKw(T). \]
Apply the matrix deviation inequality (Lecture 19) to the set \(T-T\), for which \(\operatorname{rad}(T-T)=\operatorname{diam}(T)\) and \(w(T-T)=2w(T)\): \[ \mathbb E\operatorname{diam}(AT)=\mathbb E\sup_{x\in T-T}\|Ax\|_2 \leq\sqrt m\operatorname{diam}(T)+\mathbb E\sup_{x\in T-T}\Big|\|Ax\|_2-\sqrt m\|x\|_2\Big| \leq\sqrt m\operatorname{diam}(T)+CK\big(2w(T)+\operatorname{diam}(T)\big), \] and \(\operatorname{diam}(T)\leq C'w(T)\): for any \(x,y\in T\), by the symmetrization identity of Lecture 3, \(w(T)=\tfrac12w(T-T)\geq\tfrac12\mathbb E|\langle g,x-y\rangle|=\tfrac12\sqrt{2/\pi}\,\|x-y\|_2\), and one takes the supremum over \(x,y\).
Lecture 13 had \(C(w(T)+\sqrt m\operatorname{diam}(T))\) with an unspecified constant in front of the \(\sqrt m\) term, and only for Gaussian matrices. Here the leading term has constant exactly \(1\), the matrix need only have sub-Gaussian rows, and the price is the width of \(T\) with no factor of \(\sqrt m\). In the regime \(m\gtrsim d(T)\) of Lecture 13, dividing by \(\sqrt m\) shows that \(A/\sqrt m\) changes the diameter of \(T\) by a factor \(1+O(K\sqrt{d(T)/m})\), and the phase transition at the stable dimension is now visible with its constant.
4 Where we go next
Two applications of the matrix deviation inequality have now been computed, and both turned a question about a random matrix into a question about the Gaussian width of a set: an ellipsoid for covariance estimation, a difference set for projections. The final application of the course is sparse recovery. There the set is the intersection of the kernel of a random measurement matrix with a set of sparse or nearly sparse vectors, its width is \(\sqrt{s\log(n/s)}\) by Lecture 13, and the inequality yields the number of linear measurements needed to recover an \(s\)-sparse signal by convex optimization.