Eigenvalues of correlation matrices exhibit exponential decay

I have a data-set of $P$ samples of size $N$, and noticed that the eigenvalues of the correlation matrices $A^TA$, when presented in descending order, can in many cases be described as an exponential decaying function. That is, there is a good fit linear of from $i=1..N$ to $\log|\lambda_i|$. Moreover, for several data sets I found that the exponent of the decay is fairly constant.

Is this a well-known fact or just one's tendency for finding patterns?

Obviously from PCA / SVD it tells me something about the ability to approximate the data using low-dimensional matrix.

Are there any solid mathematical results on the size of the exponent of this decay?

• To support what @random said, consider the fact that you can begin with literally any non-increasing sequence of $N$ non-negative values, form an $N$ by $N$ matrix with their square roots on its diagonal, and right- and left-multiply it by suitable orthogonal matrices to produce an $N$ by $M$ data matrix whose PCA has the original sequence for its eigenvalues. Thus there cannot be a universal law describing the decay for all datasets. In the spirit of Benford's Law, one might nevertheless contemplate the possibility that many actual datasets could exhibit such behavior.
– whuber
Feb 18 '14 at 17:24
• @whuber: and they do! I observed it many times myself, and have always assumed that for a random covariance (i.e. positive semi-definite) matrix the spectrum will show exponential decay (maybe in the limit of large $N$?), but I am not aware of precise mathematical results about it. Would be very interesting to know. Feb 18 '14 at 21:26
• @Amoeba and Uri: We need to be careful about is meant by "exponential decay." The PCA results will order the eigenvalues by decreasing magnitude, so of course the resulting scree plot will look initially like some kind of decay. It would be of interest to formally test the exponential model of decay for goodness of fit. I'm pretty sure that in most applications of PCA the goodness of fit will be poor overall, bringing us to a crucial issue: if you are allowed to "eyeball" the plot and perform the fit only for the first "few" eigenvalues, exponential decay might be a self-fulfilling prophecy.
– whuber
Feb 18 '14 at 21:49
• It may also be worth mentioning that the physicists have been interested in questions about the distributions of eigenvalues of random symmetric (or hermitean) matrices and have obtained quantitative (asymptotic) results. Google "distribution eigenvalues random matrices" or "Wigner surmise" for more.
– whuber
Feb 18 '14 at 21:53
• @whuber: The first few eigenvalues actually tend to deviate from the overall exponential decay, it is the "bulk" of the spectrum that fits to it very well. Right now I generated a random matrix (each value from standard Gaussian) of 5000*500 size and computed the spectrum of its covariance matrix. Here is the result. So I strongly suspect that the reason we see this behaviour in real datasets is simply due to noise. Feb 18 '14 at 22:01

Everything has already been pretty much figured out in the comments, thanks to @AndyW, @whuber, and @UriCohen, but I would still like to write it up as a coherent answer.

First, let me illustrate the original question. Here is the eigenspectrum of some actual real data (neural recordings) that I happen to work with right now. First few (~20-30) PCs obviously carry some signal, but after that the eigenvalues start slowly decreasing in what does seem like an exponential fashion: note that the middle part of the spectrum is almost a straight line on this log-plot. I am not showing the last part of the spectrum, because there the eigenvalues decrease pretty much to 0, due to some temporal smoothing that I used before PCA. Question is: why exponential decay?

The answer is, I believe, that any high-dimensional real data are highly contaminated by noise, so the bulk of the eigenspectrum shows the spectral behaviour of pure noise. What is the spectrum of a random covariance matrix? Turns out, there is a nice asymptotic result given by Marchenko–Pastur distribution, see the pdf of the original 1967 paper in Russian if you like.

Marchenko and Pastur tell us to consider a random data matrix of $N\times D$ size filled with independent Gaussian random values from $\mathcal{N}(0,\sigma^2)$. If $\sigma^2=1$ and $N=D$, then in the limit $N \to \infty$ the distribution of eigenvalues of its covariance matrix is given by $$\mu(x)=\frac{\sqrt{4x-x^2}}{2\pi x}.$$

Let us verify. I generated a random matrix of the $1000 \times 1000$ size, computed its covariance matrix and then calculated the eigenspectrum. The first subplot below shows the covariance matrix. The second shows the distribution (histogram) of the eigenvalues and the Marchenko-Pastur function given above. It does fit perfectly. But we are interested not so much in the distribution of eigenvalues, but in the eigenspectrum itself. If we draw 1000 values from the Marchenko-Pastur distribution (forming the spectrum) and sort them in decreasing order, then the resulting function will be given by $S(x)=(1-M(x))^{-1}$ rescaled to $[1, 1000]$, where $M(x)$ is the Marchenko-Pastur cumulative distribution function, i.e. $M(x) = \int_0^x \mu(t) dt$. The third subplot on the figure above shows the empirical spectrum vs Marchenko-Pastur fit.

It is quite a mess to compute $M(x)$, here is Wolfram Alpha's attempt. But we can note that $\mu(x)$ in the middle of its domain (around $x\approx 2$) is very well approximated by a straight line. This means that $M(x)$ will be approximately quadratic, and so its inverse $S(x) \sim \mathrm{const}-\sqrt{x}$.

In other words, the decay is not exponential at all, it is a square-root decay! However, funny enough, it is close enough to the exponential shape so that on the log-plot (see the fourth subplot above) the middle part of the spectrum looks pretty straight.

• +1 Very nice answer. The part of your contribution--which is informative, well explained, and well illustrated--seems to boil down to the fact that any reasonably smooth function (such as the plot of eigenvalues versus their rank for a large correlation matrix) will, by the very definition of "smooth," look close to linear when one focuses on any short section of it. You have established that the linearity does not hold throughout, thereby putting an important qualification on the O.P.'s observation, and have characterized the general behavior very well.
– whuber
Feb 20 '14 at 14:58
• @whuber: Thank you, but I am not sure I fully agree with your statement (or maybe I misunderstood it). Note that the plot of eigenvalues vs. rank (I call this "spectrum") is certainly smooth, but decidedly not linear: it forms a nice curl (see the low-left subplot on my main figure). What does look linear though, is the middle section of its logarithm (shown on the low-right subplot). Feb 20 '14 at 15:07
• The definition of "smooth" means at least differentiable and that, in turn, means the graph is very close to linear when looked at over small enough sections. That's all I'm saying--and that's the entire content of the OP's observation, when all is said and done.
– whuber
Feb 20 '14 at 15:17
• @whuber: But the main point of the OP's observation is that log(spectrum) is much more "linear" than raw spectrum, even though both are "smooth" and so your reasoning could be applied to either of them. Feb 20 '14 at 16:24
• You're right: my reasoning can and does apply to both observations. Its conclusion is that the OP's observation is practically without content in either case. It was only when you provided a clear quantitative assessment of the situation that any meaningful conclusions could be derived.
– whuber
Feb 20 '14 at 16:37