Minimax Lower Bounds for Covariance Estimation

Why Lower Bounds Matter

An estimator may be improved by a clever formula, but no formula can overcome a fundamental lack of information in the data. A minimax lower bound makes this limitation precise: it identifies a class of covariance matrices on which every estimator must incur a certain amount of error. The argument used here is to construct many Gaussian covariance models that are difficult to distinguish from their samples, yet are separated under the loss of interest.

The dimension matters sharply. Estimating all entries of an unrestricted covariance matrix creates many potential directions of error. In contrast, operator loss asks only for the largest directional error. These losses lead to different lower-bound rates, even on the same class of covariance matrices.

Model, Class, and Losses

Definition: Let $X_1,\ldots,X_n$ be independent $p$-dimensional observations with $X_i\sim N_p(0,\Sigma)$. For a fixed $p$, define the bounded-spectrum class $$ \mathcal{C}_p=\{\Sigma=\Sigma^\mathsf{T}: \tfrac12 I_p\preceq\Sigma\preceq\tfrac32 I_p\}. $$ Here $A\preceq B$ means that $B-A$ is positive semidefinite. An estimator $\hat{\Sigma}$ is assessed under squared Frobenius loss $L_F(\hat{\Sigma},\Sigma)=\|\hat{\Sigma}-\Sigma\|_F^2$ or squared operator loss $L_{\mathrm{op}}(\hat{\Sigma},\Sigma)=\|\hat{\Sigma}-\Sigma\|_{\mathrm{op}}^2$. The minimax risk is the infimum, over estimators, of the largest expected loss over $\Sigma\in\mathcal{C}_p$.

The eigenvalue restriction prevents the lower bound from being driven by arbitrarily large or nearly singular covariances. The data distribution is fully specified: the mean is known to be zero and the observations are Gaussian. The resulting lower bounds therefore apply even in this relatively favorable setting.

For a matrix $A$, $\|A\|_F^2$ is the sum of the squares of its entries, while $\|A\|_{\mathrm{op}}=\sup_{\|v\|_2=1}\|Av\|_2$. All rates below suppress numerical constants. For the unsquared norm losses, the corresponding lower-bound scales are the square roots of the squared-loss rates; the testing arguments also give these directly through the probability of a testing error.

The Testing Ingredient

Write $P_\Sigma$ for the joint distribution of the $n$ observations under covariance $\Sigma$. For two positive-definite covariance matrices, the Gaussian Kullback–Leibler divergence is

D(P_\Sigma\|P_\Gamma) =\frac{n}{2}\left\{\operatorname{tr}(\Gamma^{-1}\Sigma)-p-\log\det(\Gamma^{-1}\Sigma)\right\}.

In particular, if the eigenvalues of both matrices lie in $[1/2,3/2]$, then

D(P_\Sigma\|P_\Gamma)\le Cn\|\Sigma-\Gamma\|_F^2

for a numerical constant $C$. To see why, set $A=\Gamma^{-1/2}\Sigma\Gamma^{-1/2}$, whose eigenvalues lie in a fixed positive interval. The divergence is $n/2$ times the sum of $\lambda-1-\log\lambda$ over the eigenvalues $\lambda$ of $A$. On that interval, Taylor's theorem gives $\lambda-1-\log\lambda\le C_1(\lambda-1)^2$. Also, $\|A-I_p\|_F\le\|\Gamma^{-1/2}\|_{\mathrm{op}}^2\|\Sigma-\Gamma\|_F$, and $\|\Gamma^{-1/2}\|_{\mathrm{op}}^2\le2$. Combining these inequalities proves the bound.

We use the following form of Fano's inequality. Suppose $M$ distributions $P_1,\ldots,P_M$ correspond to parameters separated by at least $s$ in a metric, and a reference distribution $P_0$ satisfies $M^{-1}\sum_{j=1}^M D(P_j\|P_0)\le\alpha\log M$. Any rule for identifying which $P_j$ generated the data then has average error probability at least $1-\alpha-\log(2)/\log M$. To turn this into an estimation lower bound, decode by choosing the closest parameter to the estimator. If the estimation error is less than $s/2$, this decoder is correct. Consequently, a non-vanishing lower bound on decoding error forces estimation error of order at least $s$ with non-vanishing probability.

Worked Example: A Rank-One Family and Its Testing Difficulty

Let $u$ and $v$ be unit vectors and consider $\Sigma_u=I_p+\delta uu^\mathsf{T}$ and $\Sigma_v=I_p+\delta vv^\mathsf{T}$, with $0<\delta\le1/2$. Both matrices belong to $\mathcal{C}_p$. For unit vectors, direct multiplication gives

\|\Sigma_u-\Sigma_v\|_F^2 =\delta^2\|uu^\mathsf{T}-vv^\mathsf{T}\|_F^2 =2\delta^2\{1-(u^\mathsf{T}v)^2\}.

The operator distance is $\|\Sigma_u-\Sigma_v\|_{\mathrm{op}}=\delta\sqrt{1-(u^\mathsf{T}v)^2}$. Thus, for vectors with $|u^\mathsf{T}v|\le1/2$, the operator distance is at least $\sqrt{3}\delta/2$. The divergence from the reference model $I_p$ is also easy to check: $\Sigma_u$ has one eigenvalue $1+\delta$ and $p-1$ eigenvalues equal to 1, so

D(P_{\Sigma_u}\|P_{I_p}) =\frac n2\{\delta-\log(1+\delta)\} \le \frac{n\delta^2}{4}.

The inequality follows from $\log(1+\delta)\ge\delta-\delta^2/2$ for $\delta\ge0$. A large collection of separated directions can therefore be hard to distinguish when $n\delta^2$ is no larger than the logarithm of the collection's size. This is the mechanism behind the operator-norm lower bound.

Frobenius Loss: Many Independently Perturbed Entries

A dense covariance matrix has on the order of $p^2$ distinct entries. To expose this dimension, let $q=p(p-1)/2$, and let $B$ range over symmetric matrices with zero diagonal and entries $B_{ij}=B_{ji}\in\{-1,1\}$ for $i<j$. There is a collection of at least $\exp(cq)$ such matrices with pairwise Hamming distance at least $cq$ and operator norms at most $C\sqrt p$, for numerical constants $c,C>0$.

Here is why such a collection exists. If the upper-triangular signs are chosen independently and uniformly, then, for any fixed unit vector $x$, the random quantity $x^\mathsf{T}Bx=2\sum_{i<j}B_{ij}x_ix_j$ has sub-Gaussian tails with a numerical scale, since $4\sum_{i<j}x_i^2x_j^2\le2$. A $1/4$-net of the unit sphere has at most $9^p$ elements, and the operator norm is at most twice the maximum of $|x^\mathsf{T}Bx|$ over this net. A union bound with a sufficiently large multiple of $\sqrt p$ shows that at least a fixed fraction of all sign matrices have operator norm at most $C\sqrt p$. A greedy coding argument retains exponentially many of these matrices while discarding any pair closer than a fixed fraction of the $q$ signs.

For the retained matrices set

\Sigma_B=I_p+tB,\qquad t=\frac{a}{\sqrt{n\vee p}},

where $n\vee p=\max(n,p)$ and $a>0$ is a sufficiently small numerical constant. The perturbation has operator norm at most $aC\sqrt{p/(n\vee p)}\le aC$, so choosing $a$ small ensures $\Sigma_B\in\mathcal{C}_p$. If two sign matrices differ in $h$ upper-triangular positions, their squared Frobenius distance is $8h$. The packing therefore has squared separation at least

\|\Sigma_B-\Sigma_{B'}\|_F^2 \ge c t^2p^2 = c a^2\frac{p^2}{n\vee p}.

The divergence from $I_p$ is at most $Cn\|\Sigma_B-I_p\|_F^2\le Cna^2p^2/(n\vee p)$. Since the packing has logarithmic size of order $p^2$, choosing $a$ small makes its average divergence a small fraction of that logarithmic size. Fano's inequality then gives the following result.

Theorem: For $p\ge2$, there is a numerical constant $c>0$ such that $$ \inf_{\hat{\Sigma}}\sup_{\Sigma\in\mathcal{C}_p} E_\Sigma\|\hat{\Sigma}-\Sigma\|_F^2 \ge c\min\left\{\frac{p^2}{n},p\right\}. $$

Indeed, the squared separation just derived is of order $p^2/(n\vee p)=\min\{p^2/n,p\}$, and the testing error stays bounded away from zero. This proves the asserted lower bound. For the finitely many small dimensions where the packing constants can be adjusted, a two-point variance perturbation gives the same bound up to a universal change in $c$. $\square$

Worked Example: How the Frobenius Rate Changes With Dimension

Take $p=20$ and $n=200$. Then $p^2/n=400/200=2$, while $p=20$, so the theorem gives a squared Frobenius minimax lower bound of order $\min(2,20)=2$, up to its constant. If instead $p=400$ and $n=100$, then $p^2/n=1600$ but the bounded-spectrum class yields a lower-bound scale of $\min(1600,400)=400$. The risk scale cannot grow indefinitely like $p^2/n$ in this class: the covariance matrices themselves have squared Frobenius diameter of order $p$.

Operator Loss: Many Possible Directions

For operator loss, it suffices to vary a covariance matrix along a single direction, provided there are many possible directions to test. Unit vectors on the sphere contain a collection of size at least $\exp(cp)$ with $|u^\mathsf{T}v|\le1/2$ for distinct members. One way to obtain such a collection is to take a maximal separated set on the sphere; the volume of a spherical cap of fixed angular radius is at most $\exp(-c p)$, so exponentially many separated directions fit.

Use the rank-one covariances $\Sigma_u=I_p+\delta uu^\mathsf{T}$ from the worked example, with

\delta=a\sqrt{\min\{p/n,1\}},

where $a$ is a small fixed constant. The eigenvalues are 1 and $1+\delta$, so the class constraint holds when $a\le1/2$. Pairwise operator distances are at least $\sqrt{3}\delta/2$. Each covariance has divergence from the reference $I_p$ bounded by $n\delta^2/4\le a^2p/4$, which is a small fraction of the logarithm of the packing size when $a$ is sufficiently small. Applying Fano's inequality proves the operator result.

Theorem: For $p\ge2$, there is a numerical constant $c>0$ such that $$ \inf_{\hat{\Sigma}}\sup_{\Sigma\in\mathcal{C}_p} E_\Sigma\|\hat{\Sigma}-\Sigma\|_{\mathrm{op}}^2 \ge c\min\left\{\frac{p}{n},1\right\}. $$

The parameter separation is of order $\delta$, so a testing error bounded away from zero gives squared estimation error of order $\delta^2=\min\{p/n,1\}$. This completes the derivation. $\square$

Worked Example: Reading the Operator Bound

For $p=16$ and $n=64$, choose $a=0.1$. The rank-one construction uses $\delta=0.1\sqrt{16/64}=0.05$. Distinct directions in the packing give operator separation at least $\sqrt{3}(0.05)/2$, approximately $0.043$, while the divergence to $I_p$ is at most $64(0.05)^2/4=0.04$. The number of separated directions grows exponentially with $p=16$, so for a sufficiently small constant in the construction the data cannot reliably identify the direction. The theorem expresses this obstruction as a squared-error rate of order $16/64=0.25$, with the universal constant in the lower bound not specified by this numerical illustration.

Implications for Unregularized Estimation

When $p$ is smaller than $n$, the lower-bound scales are $p^2/n$ under squared Frobenius loss and $p/n$ under squared operator loss. Their square roots are $p/\sqrt n$ and $\sqrt{p/n}$, respectively. These rates quantify the cost of estimating an unrestricted covariance matrix: even with Gaussian observations and known mean, no estimator can uniformly improve the dimension dependence over $\mathcal{C}_p$.

The sample covariance is a natural unregularized estimator. Its error concentration, covered in the earlier material on concentration bounds for sample covariance error, can be compared with these lower bounds. The bounds here do not by themselves prove that the sample covariance attains the minimax rate in every regime. They do show that an estimator cannot remove the unrestricted problem's fundamental dimension-dependent difficulty simply by using a different unregularized formula.

When $p$ is at least as large as $n$, the lower-bound scales reach order $p$ for squared Frobenius loss and order one for squared operator loss. There is then no uniform convergence to zero in these absolute losses over the full class. This is consistent with a basic geometric fact: the sample covariance based on $n$ centered observations has rank at most $n$, so it is singular when $p>n$. The lower bounds say more than that this particular estimator is singular: they rule out uniformly accurate unrestricted covariance estimation in these losses without further assumptions.

These conclusions depend on the parameter class. A structural restriction can reduce the number of plausible covariance directions and change the minimax rates. The bounds above do not rule out improved performance on such smaller classes; they show why performance claims must state what covariance structure is assumed.

Common Pitfall

Common mistake: “The covariance matrix has about $p^2$ entries, so the operator-norm risk must also be of order $p^2/n$.” This confuses counting parameters with measuring the largest directional error. The dense packing gives a squared Frobenius lower bound of order $\min\{p^2/n,p\}$, whereas the rank-one directional packing gives a squared operator lower bound of order $\min\{p/n,1\}$. The two losses detect different features of the error matrix, so their dimension factors need not agree.

Check Your Understanding

Use the definitions and testing arguments above to answer the following questions.

Why is a bounded eigenvalue range useful when constructing a covariance-estimation lower bound?
For $p=30$ and $n=300$, what are the Frobenius and operator squared-risk lower-bound scales, ignoring constants?
In the rank-one construction, why does packing many directions lead to an operator-norm lower bound even though each perturbation has rank one?
What role does the reference distribution $P_{I_p}$ play in applying Fano's inequality?
Why does a lower bound over $\mathcal{C}_p$ not rule out better estimation on a structured covariance class?