Why a Bulk Limit Is Needed
When the number of variables is comparable to the number of observations, sample covariance matrices have many eigenvalues, and their collective pattern matters as much as any one eigenvalue. The Marchenko–Pastur law describes that pattern in a basic reference model: independent observations with identity covariance. Its limiting density gives a benchmark for interpreting sample principal components, while its support endpoints describe the edges of the bulk.
Let \(X\) be a \(p\times n\) matrix whose entries are independent \(N(0,1)\) random variables, and use the uncentered sample covariance matrix \(S=XX^{\mathsf T}/n\). Using a known zero mean rather than estimating the mean keeps the normalization explicit; centering changes the matrix by a rank-one correction and does not change this proportional-growth bulk limit. Let both \(p\) and \(n\) tend to infinity with \(p/n\to c\), where \(0<c<\infty\).
The Marchenko–Pastur Law
The support of the continuous part is \([a,b]\). Thus its upper edge is always positive, while the lower edge is positive for \(c\ne1\) and equals zero for \(c=1\). These endpoints describe the support of the limiting bulk, not a claim that every finite-sample eigenvalue lies inside it.
Worked Example: A Dimension Ratio Below One
Suppose \(p/n\to c=1/4\). Then \(\sqrt c=1/2\), so \(a=(1-1/2)^2=1/4\) and \(b=(1+1/2)^2=9/4\). Since \(c<1\), the limiting law has no atom at zero. Its density is $$ f(x)=\frac{\sqrt{(9/4-x)(x-1/4)}}{(\pi/2)x}, \qquad 1/4\le x\le 9/4. $$ The normalization in the denominator is \(2\pi c=\pi/2\). In this regime there are asymptotically more observations than variables, so the Gaussian data matrix has full row rank with probability one; correspondingly, there is no positive fraction of zero eigenvalues.
Worked Example: A Dimension Ratio Above One
Suppose \(p/n\to c=4\). The endpoints are \(a=(1-2)^2=1\) and \(b=(1+2)^2=9\). The atom at zero has mass \(1-1/c=3/4\). The continuous density is $$ f(x)=\frac{\sqrt{(9-x)(x-1)}}{8\pi x}, \qquad 1\le x\le 9. $$ The continuous part therefore has total mass \(1/4\), with the remaining \(3/4\) at zero. This agrees with rank: when \(p>n\), \(S\) has at least \(p-n\) zero eigenvalues, a fraction approaching \((p-n)/p=1-1/c=3/4\).
Deriving the Limit from Moments
A moment calculation reveals why the ratio \(c\) controls the bulk. For a positive integer \(k\), $$ \frac{1}{p}\operatorname{tr}(S^k) =\frac{1}{p n^k} \sum_{\substack{i_1,\ldots,i_k\\j_1,\ldots,j_k}} X_{i_1j_1}X_{i_2j_1}X_{i_2j_2}X_{i_3j_2}\cdots X_{i_kj_k}X_{i_1j_k}. $$ The Gaussian moment formula (Wick's rule) expresses the expectation of each product as a sum over pairings of its \(2k\) entries. Each pairing identifies some row indices and some column indices. If a pairing leaves \(r\) free row-index classes and \(s\) free column-index classes, its contribution has order \(p^r n^s/(p n^k)\). The leading pairings satisfy \(r+s=k+1\); the others are lower order as \(p,n\) grow proportionally. The leading pairings are counted by the Narayana number \(\frac1k\binom{k}{r}\binom{k}{r-1}\). Consequently, for fixed \(k\ge1\), $$ E\left[\frac{1}{p}\operatorname{tr}(S^k)\right]\longrightarrow m_k=\sum_{r=1}^k \frac{1}{k}\binom{k}{r}\binom{k}{r-1}c^{r-1}. $$
For completeness, these expected moments also concentrate around their limits. Apply the Gaussian Poincaré inequality, which bounds the variance of a differentiable function of independent standard Gaussian variables by the expected squared Euclidean norm of its gradient. For \(T_k=p^{-1}\operatorname{tr}(S^k)\), matrix differentiation gives \(\nabla_X T_k=2kS^{k-1}X/(pn)\). Thus $$ \operatorname{Var}(T_k) \le E\|\nabla_X T_k\|_F^2 =\frac{4k^2}{p^2n}\,E[\operatorname{tr}(S^{2k-1})]. $$ The same fixed-moment expansion shows \(E[\operatorname{tr}(S^{2k-1})]=O(p)\) when \(p/n\) stays bounded above and away from zero. Therefore \(\operatorname{Var}(T_k)=O((pn)^{-1})\to0\). The moments converge in probability. Their limiting sequence determines a unique probability measure: it is the Marchenko–Pastur measure, whose support is compact. This proves the bulk limit by the moment method. \(\square\)
The first few moments offer useful checks. Setting \(m_0=1\), the formula gives \(m_1=1\), \(m_2=1+c\), and \(m_3=1+3c+c^2\). In particular, the limiting mean eigenvalue is one, as expected from the identity population covariance.
Worked Example: Checking the First Moments at \(c=2\)
For \(c=2\), the first three positive moments are $$ m_1=1,\qquad m_2=1+2=3,\qquad m_3=1+3(2)+2^2=11. $$ The second-moment formula has terms \(r=1,2\): \(1+ c\). For the third, the terms are \(1+3c+c^2\). These checks are useful because they include the zero eigenvalues in the \(p\)-normalized ESD. In particular, the positive eigenvalues must carry enough mass and size for the overall mean still to equal one despite the atom at zero.
From Moments to Density and Edges
Let \(M(t)=1+\sum_{k\ge1}m_k t^k\) be the moment generating power series near \(t=0\). The Narayana-number formula above yields \(M(t)=1+tM(t)\{1-c+cM(t)\}\). Define the Stieltjes transform using the convention \(m(z)=\int (x-z)^{-1}\,d\mu_c(x)\), for \(z\) outside the support. Since \(M(1/z)=-z\,m(z)\), substitution gives the quadratic equation $$ c z m(z)^2+(z+c-1)m(z)+1=0. $$ The solution that satisfies \(m(z)\sim-1/z\) as \(|z|\to\infty\) is $$ m(z)=\frac{-(z+c-1)+\sqrt{(z-a)(z-b)}}{2cz}, $$ where the square-root branch is chosen to behave like \(z\) at infinity. This condition selects a probability-measure transform rather than the other algebraic root.
For \(x\) strictly between \(a\) and \(b\), the upper-half-plane boundary value of the square root has imaginary part \(\sqrt{(b-x)(x-a)}\). The Stieltjes inversion formula states that the continuous density is \(\pi^{-1}\lim_{\eta\downarrow0}\operatorname{Im}m(x+i\eta)\). Substitution gives exactly \(\sqrt{(b-x)(x-a)}/(2\pi c x)\). The possible point mass at zero can be checked from rank: an iid Gaussian \(p\times n\) matrix has rank \(\min(p,n)\) almost surely. If \(p>n\), the fraction of zero eigenvalues is exactly \((p-n)/p\); if \(p\le n\), it is zero. These fractions converge to \((1-1/c)_+\).
The square root in the density produces a soft, square-root edge whenever the endpoint is positive. For example, as \(x\uparrow b\), the density is asymptotic to \(\sqrt{b-a}\sqrt{b-x}/(2\pi c b)\). If \(c\ne1\), the same square-root form holds as \(x\downarrow a>0\), with leading term \(\sqrt{b-a}\sqrt{x-a}/(2\pi c a)\). At \(c=1\), however, \(a=0\) and the denominator \(x\) changes the behavior: $$ f(x)=\frac{1}{2\pi}\sqrt{\frac{4-x}{x}} \sim\frac{1}{\pi\sqrt{x}} \qquad\text{as }x\downarrow0. $$ This divergence is integrable and is not an atom at zero. It is called a hard edge: the density grows near zero, but a point mass would require a nonzero fraction of eigenvalues to be exactly zero in the limit.
Worked Example: The Hard Edge at \(c=1\)
When \(p/n\to1\), the endpoints are \(a=0\) and \(b=4\), and there is no zero atom because \((1-1/c)_+=0\). The density reduces to \(f(x)=\sqrt{(4-x)/x}/(2\pi)\) on \(0<x<4\). Near zero, \(\sqrt{4-x}\to2\), so \(f(x)\sim1/(\pi\sqrt{x})\). The density is unbounded but integrable: \(\int_0^\epsilon x^{-1/2}\,dx=2\sqrt{\epsilon}\), which tends to zero with \(\epsilon\). Thus a large density close to zero does not imply positive mass at zero.
Common Pitfalls and Interpretation
- Using the wrong normalization. Here the ESD divides by \(p\), and the covariance matrix uses \(1/n\). With the alternative \(1/p\) scaling, the eigenvalue scale and formula change.
- Confusing a zero atom with a small eigenvalue. At \(c>1\), exactly \(p-n\) eigenvalues are zero almost surely in this Gaussian model. At \(c=1\), there is no limiting atom, even though the density diverges near zero.
- Forgetting the atom when \(c>1\). The density alone has mass \(1/c\), not one. The missing mass is \((1-1/c)\delta_0\).
- Treating the support edges as an outlier analysis. The endpoints \(a,b\) describe the identity-covariance bulk. They do not, by themselves, characterize isolated eigenvalues under a spiked population covariance.
The law matters because it separates collective sampling variation from structure in the population. In high-dimensional PCA, a sample eigenvalue is not automatically evidence of a distinct signal: even identity-covariance data produce a nontrivial eigenvalue spread. The ratio \(p/n\) determines that spread and whether rank deficiency creates an exact zero mass. A reliable interpretation starts by identifying the normalization, the limiting ratio, and whether an apparent feature belongs to the bulk, an atom, or a finite-sample fluctuation.
Check Your Understanding
Use the definitions and derivation above to answer these questions.
- For \(p/n\to c=1/9\), calculate both support endpoints and state whether the limiting law has an atom at zero.
- When \(p/n\to c=3\), what fraction of the ESD is at zero, and what fraction is in the continuous part?
- Why does an unbounded density near zero at \(c=1\) not imply a point mass at zero?
- Using the first three moment formulas, calculate \(m_1,m_2,m_3\) when \(c=1/2\).
- Explain why the support endpoints depend on the ratio \(p/n\), rather than on \(p\) or \(n\) alone.