Concentration Bounds for Sample Covariance Error

Why covariance error depends on dimension

A covariance matrix summarizes variability in every direction at once. Estimating it well in a fixed direction is not the same as estimating it well simultaneously across all directions: as the dimension grows, there are more directions in which sampling error can be large. The operator norm captures the largest such error, making it a natural measure of spectral covariance error.

The characteristic high-probability rate has two parts: a square-root term, $\sqrt{p/n}$, and a linear term, $p/n$. The first is the familiar fluctuation scale when dimension is modest relative to sample size. The second matters when dimension is comparable to, or larger than, sample size. Both are multiplied by a scale determined by the covariance and, for sub-Gaussian data, by the strength of the tail bound.

Definition: Let $X_1,\ldots,X_n$ be independent, identically distributed random vectors in $\mathbb{R}^p$, with mean $\mu$ and covariance matrix $\Sigma=E[(X_i-\mu)(X_i-\mu)^T]$. The sample-centered covariance, using divisor $n$, is $$ S=\frac{1}{n}\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^T,\qquad \bar X=\frac{1}{n}\sum_{i=1}^n X_i. $$ Its operator-norm error is $\|S-\Sigma\|_{\mathrm{op}}$, where $\|A\|_{\mathrm{op}}=\sup_{\|u\|_2=1}\|Au\|_2$. For a symmetric matrix, this is also $\sup_{\|u\|_2=1}|u^TAu|$.

We use divisor $n$ so that the concentration argument has a direct empirical-average form. The unbiased covariance with divisor $n-1$ differs by the factor $n/(n-1)$; for $n\ge 2$, this changes the bound only by an additional scale term of order $\|\Sigma\|_{\mathrm{op}}/n$.

Gaussian concentration: the net argument

First suppose $Z_i=X_i-\mu$ are Gaussian with distribution $N(0,\Sigma)$. It is helpful to separate the covariance formed around the true mean, $S_0=n^{-1}\sum_i Z_iZ_i^T$, from the sample-centered covariance. For each fixed unit vector $u$, the projections $u^TZ_i$ are independent $N(0,u^T\Sigma u)$ variables. Thus their squared, normalized sum has a chi-square distribution.

Theorem: There is a universal constant $C$ such that, for every $t\ge 0$, Gaussian observations satisfy $$ P\!\left\{\|S-\Sigma\|_{\mathrm{op}}> C\|\Sigma\|_{\mathrm{op}}\left[ \sqrt{\frac{p+t}{n}}+\frac{p+t}{n} \right]\right\}\le 2e^{-t}. $$ The result holds for all positive integers $p,n$; it does not assert that the error is small when $p$ is much larger than $n$.

Here is a derivation. For a chi-square variable $Q$ with $n$ degrees of freedom, the standard chi-square tail inequality is $P(|Q/n-1|>2\sqrt{x/n}+2x/n)\le 2e^{-x}$ for $x\ge0$. Applying it to each fixed $u$ gives $$ P\!\left\{|u^T(S_0-\Sigma)u|> 2(u^T\Sigma u)\left(\sqrt{\frac{x}{n}}+\frac{x}{n}\right)\right\} \le 2e^{-x}. $$ Since $u^T\Sigma u\le\|\Sigma\|_{\mathrm{op}}$, the scale can be replaced by $\|\Sigma\|_{\mathrm{op}}$.

To control all directions, take a $1/4$-net $\mathcal N$ of the Euclidean unit sphere: every unit vector is within distance $1/4$ of a point in $\mathcal N$. Such a net can be chosen with at most $9^p$ points. For any symmetric matrix $A$, $\|A\|_{\mathrm{op}}\le 2\max_{v\in\mathcal N}|v^TAv|$. Indeed, choose a unit vector $u$ attaining the supremum of $|u^TAu|$, and a net point $v$ with $\|u-v\|_2\le1/4$. Then $|u^TAu-v^TAv|\le 2\|A\|_{\mathrm{op}}\|u-v\|_2\le\|A\|_{\mathrm{op}}/2$; rearranging proves the claim.

Apply the fixed-direction inequality to all net points and use the union bound. Taking $x=p\log 9+t$ makes the failure probability at most $2e^{-t}$. The net inequality then yields the stated form for $\|S_0-\Sigma\|_{\mathrm{op}}$, after adjusting a universal constant.

Centering at $\bar X$ subtracts a rank-one matrix: $$ S=S_0-\bar Z\bar Z^T,\qquad \bar Z=\frac{1}{n}\sum_{i=1}^n Z_i. $$ For Gaussian data, $\bar Z\sim N(0,\Sigma/n)$. Applying the same net and chi-square argument to $\|\bar Z\|_2^2=\sup_{\|u\|_2=1}(u^T\bar Z)^2$ gives, with failure probability at most $2e^{-t}$, $\|\bar Z\|_2^2\le C\|\Sigma\|_{\mathrm{op}}(p+t)/n$. Since $\|\bar Z\bar Z^T\|_{\mathrm{op}}=\|\bar Z\|_2^2$, this term is absorbed by the displayed rate. Combining the two bounds, with a harmless adjustment of constants and tail parameter, proves the theorem. $\square$

Worked Example: Gaussian observations with moderate dimension

Suppose $X_i\sim N(\mu,4I_{100})$, independently, and $n=1000$. The scale is $\|\Sigma\|_{\mathrm{op}}=4$. At tail parameter $t=\log 20$, the dimension-plus-tail quantity is $100+\log 20$, approximately $103.0$. Therefore $$ \sqrt{\frac{p+t}{n}}\approx\sqrt{0.1030}\approx0.321,\qquad \frac{p+t}{n}\approx0.103. $$

The theorem bounds the error, except on an event of probability at most $2e^{-t}=0.1$, by a universal constant times $4(0.321+0.103)$, or approximately $1.70C$. This is a bound, not a prediction that the realized error equals that number. The calculation shows how both dimension-dependent terms enter, and why the covariance scale multiplies them.

Sub-Gaussian observations

Gaussianity is not essential. For a centered random variable $W$, define its sub-Gaussian norm by $\|W\|_{\psi_2}=\inf\{c>0:E[\exp(W^2/c^2)]\le2\}$. Assume that the centered vector $Z=X-\mu$ obeys the covariance-relative condition $\|u^TZ\|_{\psi_2}\le K\sqrt{u^T\Sigma u}$ for every $u\in\mathbb{R}^p$, where $K\ge1$. This condition controls the tails of every linear projection, not just the coordinate marginals.

Theorem: Under the preceding condition, there is a universal constant $C$ such that for all $t\ge0$, $$ P\!\left\{\|S-\Sigma\|_{\mathrm{op}}> CK^2\|\Sigma\|_{\mathrm{op}}\left[ \sqrt{\frac{p+t}{n}}+\frac{p+t}{n} \right]\right\}\le 2e^{-t}, $$ after increasing $C$ if necessary. The observations must be independent and identically distributed, with finite covariance as specified.

To derive the bound, fix a unit vector $u$ and set $W_i=u^TZ_i$. The sub-Gaussian assumption gives $\|W_i\|_{\psi_2}\le K\sqrt{u^T\Sigma u}$. A standard norm relation is $\|W_i^2-EW_i^2\|_{\psi_1}\le C_0\|W_i\|_{\psi_2}^2$, where $\|V\|_{\psi_1}=\inf\{c>0:E[\exp(|V|/c)]\le2\}$. This follows by bounding the exponential tail of $W_i^2$ from the definition of $\psi_2$, and centering changes the bound by at most a universal factor.

Bernstein's inequality for independent, centered sub-exponential variables states that if $\|Y_i\|_{\psi_1}\le b$, then $P(|n^{-1}\sum_iY_i|>C_1b(\sqrt{x/n}+x/n))\le2e^{-x}$. Apply it to $Y_i=W_i^2-EW_i^2$, with $b=C_0K^2u^T\Sigma u$. This gives the fixed-direction covariance bound with scale at most $C_0K^2\|\Sigma\|_{\mathrm{op}}$. The same $1/4$-net, cardinality bound, and union bound used above yield the claimed rate for $S_0-\Sigma$.

For centering, each projection $u^T\bar Z$ is sub-Gaussian with norm at most $K\sqrt{u^T\Sigma u/n}$. The sub-Gaussian tail inequality and a $1/4$-net imply $\|\bar Z\|_2^2\le CK^2\|\Sigma\|_{\mathrm{op}}(p+t)/n$, except on an event with probability bounded by a constant times $e^{-t}$. Adding this rank-one correction proves the sample-centered result; constants can be enlarged to obtain the displayed tail form. $\square$

Worked Example: Reading the two terms when dimension exceeds sample size

Suppose $p=800$, $n=200$, $\|\Sigma\|_{\mathrm{op}}=3$, and the relative sub-Gaussian constant is $K=1.5$. Ignoring the tail-parameter addition for a rate comparison, $$ \sqrt{\frac{p}{n}}=\sqrt{4}=2,\qquad \frac{p}{n}=4,\qquad K^2\|\Sigma\|_{\mathrm{op}}=2.25(3)=6.75. $$

The bound's rate factor is $6.75(2+4)=40.5$, up to the universal constant. Its size is not a claim that the error must exceed the scale of $\Sigma$; it signals that this general concentration bound does not guarantee small operator-norm error in this regime. In particular, retaining only $\sqrt{p/n}$ would omit a term that is larger here.

What sets the scale?

The natural scale in the covariance-relative formulation is $\|\Sigma\|_{\mathrm{op}}$, the largest variance over unit-length linear combinations of coordinates. It is not always the largest coordinate variance. If the sub-Gaussian assumption is instead given in marginal-scale form, $\|u^TZ\|_{\psi_2}\le K\tau\|u\|_2$, where $\tau^2=\max_j\operatorname{Var}(Z_j)$, the same proof gives the bound with $K^2\tau^2$ in place of $K^2\|\Sigma\|_{\mathrm{op}}$. This is a stronger tail assumption when directions can combine many coordinates.

A marginal variance bound by itself does not generally give that replacement: a covariance matrix can have every diagonal entry at most $\tau^2$ while its largest eigenvalue is as large as $p\tau^2$. For Gaussian data, the covariance-relative result remains valid with $\|\Sigma\|_{\mathrm{op}}$; if $\Sigma$ is diagonal and its entries are at most $\tau^2$, then $\|\Sigma\|_{\mathrm{op}}\le\tau^2$, so the marginal scale suffices.

Worked Example: Why marginal variances alone can miss a direction

Let $Z=\tau G\mathbf{1}$, where $G\sim N(0,1)$ and $\mathbf{1}$ is the $p$-vector of ones. Every coordinate has variance $\tau^2$, so the largest marginal variance is $\tau^2$. But $\Sigma=\tau^2\mathbf{1}\mathbf{1}^T$, whose only nonzero eigenvalue is $p\tau^2$. Indeed, in the unit direction $u=\mathbf{1}/\sqrt p$, $$ \operatorname{Var}(u^TZ)=\operatorname{Var}(\tau\sqrt p\,G)=p\tau^2. $$

Thus the largest variance in a direction is $p$ times the marginal variance. A rate scaled only by $\tau^2$ would overlook this direction unless an additional assumption controls dependence or all directional tails.

Why both dimension terms matter

The net proof explains the appearance of dimension: controlling a fixed direction gives a tail probability, and a net of the sphere has exponentially many points in $p$. Setting the tail parameter at order $p+t$ makes the union bound work. The square-root term comes from the square-root part of one-dimensional concentration; the linear term comes from its large-deviation part. For $p/n$ small, the square-root term is larger. When $p/n$ is large, the linear term cannot be dropped.

These are upper bounds for broad distribution classes, not exact descriptions of every covariance estimator's error. They also concern the unregularized sample covariance. Whether a structured or regularized estimator can do better depends on additional assumptions; those assumptions are not part of these concentration statements.

Common mistake

Dropping $p/n$ without checking the regime. The reasoning “covariance error is always of order $\sqrt{p/n}$” keeps only the small-deviation term. The correct general bound includes $\sqrt{(p+t)/n}+(p+t)/n$; the second term is material when $p+t$ is not small relative to $n$.
Using the largest marginal variance as the scale automatically. Marginal variances control coordinate directions, not necessarily all linear combinations. Use $\|\Sigma\|_{\mathrm{op}}$, or state a stronger directional sub-Gaussian assumption that justifies a marginal-scale bound.
Ignoring sample centering. The sample-centered matrix is $S_0-\bar Z\bar Z^T$, not $S_0$. Its correction has operator norm $\|\bar Z\|_2^2$, which is of the same dimensional order and must be accounted for.

Key takeaway: Under Gaussian or suitably sub-Gaussian tails, sample covariance operator error is bounded at scale $\text{covariance scale}\times[\sqrt{(p+t)/n}+(p+t)/n]$ with high probability. The square-root term describes moderate fluctuations; the linear term protects the bound outside that regime.

Check Your Understanding

Use the definitions and bounds above to answer the following questions.

Why does a net over the unit sphere introduce a quantity proportional to $p$ into the tail parameter?
For $p/n=1/16$, which of $\sqrt{p/n}$ and $p/n$ is larger, and by what factor?
Explain why sample centering subtracts a rank-one matrix and state its operator norm.
Can a bound involving only the largest coordinate variance be used for every covariance matrix? Give the directional reason for your answer.
When $p$ is much larger than $n$, what does the linear term imply about the usefulness of this general upper bound?