# How many samples are needed to estimate a p-dimensional covariance matrix?

In general, how many points are needed to estimate a p-dimensional covariance matrix? Does it depend on how the data are spread out across the different dimensions? Does it depend on the true distribution of the data? Thank you!

• Two cases will do it, period: this can easily be proven by examining the formula for the usual sample covariance estimate. If you want a more precise answer you will need to stipulate the level of accuracy needed in the estimate.
– whuber
Mar 14, 2014 at 18:13
• Thank you. I see what you mean, since having two points always allows you to compute a sample mean. Then, one can take differences between the original points and sample means to compute covariances. When I try simulating in MATLAB, however, the covariance matrix I end up with is singular. For example, if I do: x = [1 2 3; 4 3 6; 7 8 10]; eig(cov(x)), one of the eigenvalues is 0. On the other hand, if I do x = [1 2 3; 4 3 6; 7 8 10; -1 2 5], none of the eigenvalues are zero. How many points do I need to ensure that the matrix I estimate is not singular? Thanks! Mar 14, 2014 at 18:51
• It depends on how you simulate the matrix. Your matrices do not appear to be simulated randomly, nor according to any clear rule, so it would be impossible to answer that question. What, then, is the question you really want to ask? Are you interested in creating random data with nonsingular covariance matrices? For what purpose?
– whuber
Mar 14, 2014 at 19:03
• Not for any particular purpose. I just want to understand some principles behind estimating covariance matrices. In both my examples above, the dimensions are not highly correlated (i.e. there is not a deterministic relationship between any of the dimensions). When I have three data points and three dimensions, the covariance matrix is singular; when I have four data points and three dimensions, the covariance matrix is not. Does that mean that as long as the dimensions aren't correlated that having p + 1 data points will guarantee a non-singular covariance matrix? Thank you! Mar 14, 2014 at 19:17
• You are exploring two different situations. Estimated covariance matrices will be singular if and only if the data are collinear. With $p$ data points in $p$ dimensions, if their values are chosen independently and randomly according to some common underlying continuous distribution, the chance is extremely high that the covariance matrix will not be singular: this merely generalizes the fact that a randomly chosen number will not be zero with high probability. If the data are generated by some other process they can have a singular covariance matrix no matter how much there are.
– whuber
Mar 14, 2014 at 19:35

As @whuber has commented (+1 to his comment) you need "$p$ data points in $p$ dimensions" under the assumption those points are independent. Nevertheless as you correctly recognise this will depended on the underlying distribution of your data as well as your sample size. That is because in finite samples sizes you have finite sample effects. That means that, in very approximate manner, your random sample exhibits properties (usually regularities in the form of collinearity) that should not be there. This is not something new: Finite or small sample corrections are something that is done ubiquitously in Statistics, for example the very popular Akaike Information Criterion (AIC) has a very easy to compute version that is corrected for finite samples: the AICc (that is unfortunately underused - AICc corrects for deviations from normalities, not collinearities). So how bad things might be?

[A quick note: A singular covariance matrix is essentially one that is not positive definite (PD). You can check if a matrix is PD by checking if it has a Cholesky decomposition. That is much faster than using eigendecomposition or SVD.]

Let's say is looking at a random sample $S_1$ such that $s_1 \sim U[0,1]$ and another sample $S_2$ such that $s_2 \sim T_{\nu=1}$ ($\nu$ being the degrees of freedom). How would thing go with a relatively small (<200) sample? Well... not splendid! Let's simulate something like this (in MATLAB):

rng(1234)
N = 200;
M = 200;
FailsU = zeros(N,M);
FailsT = zeros(N,M);

for i = 1:M
for j = 1:N;
K = cov(rand(i));
[L,p] = chol(K,'lower');
if(p)
FailsU(j,i) = 1;
end
K = cov(random('t',1,[i,i]));
[L,p] = chol(K,'lower');
if(p)
FailsT(j,i) = 1;
end
end
end


About 50% of the time you will get a non-PD matrix. That's definitely not good. What about if we had $p+1$ samples though?

% Same initializations as above

for i = 1:M
for j = 1:N;
K = cov(rand(i+1,i));
[L,p] = chol(K,'lower');
if(p)
Fails1(j,i) = 1;
end
K = cov(random('t',1,[i+1,i]));
[L,p] = chol(K,'lower');
if(p)
FailsT1(j,i) = 1;
end
end
end


Yeap we are fine! Why all this trouble though?

Let see the $p \times p$ version first: When you are saying that a matrix is non-invertible or it has zero eigenvalues you are saying it is rank deficient. That is because the rank of a covariance matrix can be thought as the number of non-zero eigenvalues. In addition, when you compute the product $S_x^T S_x$ you can get at most a matrix of rank $p$ because the rank of a product of matrices is at most equal to the rank of any matrix in the product. That means that the covariance matrix of your $p$-dimensional sample has at most rank $p$. Therefore for any number of samples $N$, the rank of $S_x$ is at most $p$. From this we can assume that even small deviations from complete randomness would make our covariance matrix rank deficient (as seen in the first simulation). OK, fine why does $p+1$ seems to work so nicely?

Now let's see the $p+1 \times p$ version: Think of how you estimate a sample covariance matrix: While in a quick manner we can write : $K = \frac{1}{N-1} S_x^T S_X$ (because we assumed $S_x$ to have mean 0, we should properly write things as: $K = \frac{1}{N-1} \Sigma_{i=1}^N (S_{x(i)} - \hat{\mu})(S_{x(i)} - \hat{\mu})^T$ where $\hat{\mu}$ is the sample mean. But what about the rank of this $(S_{x(i)} - \hat{\mu})(S_{x(i)} - \hat{\mu})^T$ matrix? Well... that is 1! Not only that but exactly because we subtracted $\hat{\mu}$ we reduced the original rank of the matrix $S_x$ to begin with! So as we add up points ($N$ gets larger) we have more chances to have a full-rank matrix. For the current case, exactly the covariance in our sample is completely diagonal, just adding another point ensure that it was extremely unlikely (not guaranteed) to have a rank deficiency.

By construction, a covariance matrix of dimension $n$, computed with $n$ data points is singular. Consider the following $2 \times2$ example.

Data: $x_1, x_2$ and $y_1, y_2$.

Let

$x_1 - \bar{x} = -(x_2 - \bar{x}) = a$

$y_1 - \bar{y} = -(y_2 - \bar{y}) = b$

The covariance matrix (ignoring bias/ML issues) is:

$$\pmatrix{2 a^2 & 2ab \\ 2 ab & 2b^2}$$

which is singular.

• I'm not sure that this answers the question. May 17, 2015 at 18:17
• This seems more like part of an answer to this question than a complete answer. Can you expand this to provide a complete answer? May 17, 2015 at 18:20

The short answer to the question is: it depends on the distribution of the data.

Also, it depends on your definition of estimation: should they be close in the l_2 norm or do we use other metrics?

As specific case, in the case of the operator norm, we can ask how many samples do I need if the data is a p-dimensional Gaussian with variace $\Sigma$. the answer is fortunately $O(p)$. You might need more samples if the data has other distributions.

Best source of information that I was able to find online this paper, along with the presentation of the same author here.

• Could you please provide some basis for these assertions, in light of the fact that they contradict other answers? You seem to have in mind a more specific version of the question in which the covariance matrix needs to be estimated with a given degree of precision and accuracy.
– whuber
Sep 11, 2018 at 14:02