Sampling Covariance Matrices for PCA Inference

Gaussian Sampling for PCA

Principal component analysis is often introduced as a decomposition of one observed data matrix. In statistical inference, that matrix is instead viewed as a random sample from a population. Its covariance matrix is then random too: even when the population covariance is fixed, a different sample generally produces a different sample covariance and different estimated principal directions.

This tutorial sets up the fixed-dimensional Gaussian sampling model and establishes an exact distributional representation for the centered sample covariance. The representation is useful because it makes the source of sampling variation explicit. Throughout, $p$ is the fixed number of variables, $n\geq 2$ is the number of independent observations, and the population covariance matrix is assumed positive definite.

Definition: Let $X_1,\ldots,X_n$ be independent random vectors in $\mathbb{R}^p$ with common distribution $N_p(\mu,\Sigma)$, where $\mu\in\mathbb{R}^p$ and $\Sigma$ is a positive-definite $p\times p$ covariance matrix. Define the sample mean and centered sample covariance by $$ \bar X=\frac{1}{n}\sum_{i=1}^n X_i, \qquad S=\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^{\mathsf T}. $$ The divisor $n-1$ makes $S$ unbiased for $\Sigma$.

The Gaussian assumption specifies the sampling distribution, not a requirement that the observed data look exactly Gaussian. The model says that each observation has the same mean and covariance and that observations are independent. In the fixed-dimensional setting, $p$ is held fixed when sample size is considered; this tutorial uses the model for finite-sample distributional results, without requiring a limiting approximation.

The Population Principal Components

The population covariance matrix has an orthonormal eigendecomposition $$ \Sigma=U\Lambda U^{\mathsf T}, \qquad \Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_p), \qquad \lambda_1\geq\cdots\geq\lambda_p>0, $$ where the columns $u_1,\ldots,u_p$ of $U$ are orthonormal eigenvectors. The population principal component score in direction $u_j$ is $u_j^{\mathsf T}(X-\mu)$; its variance is $\lambda_j$. Thus, the population eigenvalues describe variation along population principal directions, while $S$ provides a sample-based estimate of the covariance matrix.

If an eigenvalue is repeated, its individual eigenvectors need not be uniquely determined; the corresponding eigenspace is the meaningful object. If an eigenvalue is simple, its eigenvector is still defined only up to sign. These are properties of the population decomposition, not sampling errors. Sampling uncertainty arises because $S$ varies from sample to sample.

Worked Example: Population Axes and Variances

Suppose $X\sim N_2(\mu,\Sigma)$, with $\mu=(10,20)^{\mathsf T}$ and $$ \Sigma= \begin{pmatrix} 2&1\\ 1&2 \end{pmatrix}. $$ Find the population eigenvalues and principal directions.

The characteristic polynomial is $$ \det(\Sigma-\lambda I)=(2-\lambda)^2-1 =(1-\lambda)(3-\lambda), $$ so the eigenvalues are $\lambda_1=3$ and $\lambda_2=1$.

For $\lambda_1=3$, the equation $(\Sigma-3I)u=0$ gives $u_1$ proportional to $(1,1)^{\mathsf T}$. For $\lambda_2=1$, it gives $u_2$ proportional to $(1,-1)^{\mathsf T}$. Normalizing, $$ u_1=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\1\end{pmatrix}, \qquad u_2=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\-1\end{pmatrix}. $$ The variance of $u_1^{\mathsf T}(X-\mu)$ is $3$, and the variance of $u_2^{\mathsf T}(X-\mu)$ is $1$. The total variance is $3+1=4$, equal to $\operatorname{tr}(\Sigma)=2+2=4$.

The Wishart Representation

For a positive-definite matrix $\Sigma$ and positive integer $\nu$, the Wishart distribution $W_p(\Sigma,\nu)$ is the distribution of a sum of $\nu$ independent outer products of $N_p(0,\Sigma)$ random vectors. In particular, if $Z_1,\ldots,Z_\nu$ are independent $N_p(0,\Sigma)$, then $\sum_{a=1}^{\nu}Z_aZ_a^{\mathsf T}\sim W_p(\Sigma,\nu)$. This definition fixes both the scale matrix and the degrees of freedom.

Theorem: Under the Gaussian sampling model above, the centered scatter matrix and sample covariance satisfy $$ (n-1)S=\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^{\mathsf T} \sim W_p(\Sigma,n-1), \qquad S\sim \frac{1}{n-1}W_p(\Sigma,n-1). $$

Proof. Arrange the observations as rows of an $n\times p$ matrix $X$, and let $\mathbf{1}$ be the $n$-vector of ones. The centering matrix is $H=I_n-\mathbf{1}\mathbf{1}^{\mathsf T}/n$; multiplying $X$ by $H$ on the left subtracts the sample mean from every row. Choose an $n\times(n-1)$ matrix $Q$ whose columns form an orthonormal basis for the subspace perpendicular to $\mathbf{1}$. Then $Q^{\mathsf T}Q=I_{n-1}$ and $QQ^{\mathsf T}=H$. Consequently, the centered scatter matrix is $$ X^{\mathsf T}HX=(Q^{\mathsf T}X)^{\mathsf T}(Q^{\mathsf T}X). $$

Each column of $X$ is Gaussian, and $Q^{\mathsf T}X$ is a linear transformation of the jointly Gaussian entries of $X$. Its mean is zero because $Q^{\mathsf T}\mathbf{1}=0$. The covariance between transformed rows $a$ and $b$, in variable coordinates, is $$ \operatorname{Cov}\big((Q^{\mathsf T}X)_{a,\cdot}^{\mathsf T},(Q^{\mathsf T}X)_{b,\cdot}^{\mathsf T}\big) =(Q^{\mathsf T}Q)_{ab}\Sigma =\begin{cases} \Sigma,&a=b,\\ 0,&a\ne b. \end{cases} $$ Since these rows are jointly Gaussian, zero cross-covariance makes them independent. There are $n-1$ such rows, each distributed as $N_p(0,\Sigma)$. Their outer products therefore sum to $W_p(\Sigma,n-1)$, proving the scatter-matrix statement. Dividing by $n-1$ gives the result for $S$. $\square$

The centering is essential to the degrees of freedom: one direction in the $n$-dimensional observation index is removed when the sample mean is estimated. This is why the centered scatter has $n-1$, rather than $n$, independent Gaussian contributions. The theorem also shows that the distribution of $S$ depends on $\mu$ only through the fact that it has been removed by centering.

Worked Example: Exact Entrywise Uncertainty

Let $n=5$ and suppose $$ \Sigma= \begin{pmatrix} 4&0\\ 0&1 \end{pmatrix}. $$ The Wishart representation gives $4S\sim W_2(\Sigma,4)$. The entries of $S$ have expectation $E[S_{ij}]=\Sigma_{ij}$. Their exact covariances, derived below, are $$ \operatorname{Cov}(S_{ij},S_{kl}) =\frac{\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}}{n-1}. $$

For the first diagonal entry, $$ \operatorname{Var}(S_{11}) =\frac{4\cdot4+4\cdot4}{4}=8. $$ For the second, $$ \operatorname{Var}(S_{22}) =\frac{1\cdot1+1\cdot1}{4}=\frac12. $$ For the off-diagonal entry, $$ \operatorname{Var}(S_{12}) =\frac{\Sigma_{11}\Sigma_{22}+\Sigma_{12}\Sigma_{21}}{4} =\frac{4\cdot1+0}{4}=1. $$ Finally, $$ \operatorname{Cov}(S_{11},S_{22}) =\frac{\Sigma_{12}\Sigma_{12}+\Sigma_{12}\Sigma_{12}}{4}=0. $$ Thus $S_{11}$ and $S_{22}$ are individually unbiased for $4$ and $1$, but their sampling variances are $8$ and $1/2$; the sample off-diagonal entry has mean zero and variance $1$. A zero population covariance does not force a sample covariance to be zero.

Means and Covariances of the Sample Covariance

Theorem: Under the stated Gaussian model, for $1\leq i,j,k,l\leq p$, $$ E[S_{ij}]=\Sigma_{ij}, \qquad \operatorname{Cov}(S_{ij},S_{kl}) =\frac{\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}}{n-1}. $$ In particular, $\operatorname{Var}(S_{ij})=(\Sigma_{ii}\Sigma_{jj}+\Sigma_{ij}^2)/(n-1)$.

Proof. By the Wishart representation, write $(n-1)S=\sum_{a=1}^{n-1} Z_aZ_a^{\mathsf T}$, where the $Z_a$ are independent $N_p(0,\Sigma)$. For one such vector, $E[Z_iZ_j]=\Sigma_{ij}$, and the Gaussian fourth-moment identity is $$ E[Z_iZ_jZ_kZ_l] =\Sigma_{ij}\Sigma_{kl}+\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}. $$ This identity follows by differentiating the Gaussian moment-generating function $E[\exp(t^{\mathsf T}Z)]=\exp(t^{\mathsf T}\Sigma t/2)$ four times at $t=0$. Subtracting $E[Z_iZ_j]E[Z_kZ_l]$ gives $$ \operatorname{Cov}(Z_iZ_j,Z_kZ_l) =\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}. $$ Independence across $a$ makes the covariance of the corresponding sums equal to $n-1$ times this quantity. Dividing each sum by $n-1$ yields the covariance formula for $S$. Likewise, each summand has mean $\Sigma_{ij}$, so the divided sum has mean $\Sigma_{ij}$. Setting $k=i$ and $l=j$ gives the variance formula. $\square$

These formulas quantify matrix-entry uncertainty without a large-sample approximation. For example, they also give the expected squared Frobenius error: $$ E\!\left[\|S-\Sigma\|_F^2\right] =\sum_{i=1}^p\sum_{j=1}^p\operatorname{Var}(S_{ij}) =\frac{(\operatorname{tr}\Sigma)^2+\operatorname{tr}(\Sigma^2)}{n-1}. $$ The final equality follows by summing $(\Sigma_{ii}\Sigma_{jj}+\Sigma_{ij}^2)/(n-1)$ over all $i,j$. It displays the roles of sample size and population scale: for fixed $\Sigma$, the expected squared matrix error decreases in proportion to $1/(n-1)$.

Worked Example: A Concrete Centered Covariance

Consider three observed two-variable vectors $x_1=(1,0)^{\mathsf T}$, $x_2=(-1,0)^{\mathsf T}$, and $x_3=(0,1)^{\mathsf T}$. Their mean is $$ \bar x=\frac{x_1+x_2+x_3}{3} =\begin{pmatrix}0\\1/3\end{pmatrix}. $$ The centered vectors are $(1,-1/3)^{\mathsf T}$, $(-1,-1/3)^{\mathsf T}$, and $(0,2/3)^{\mathsf T}$.

Their outer products sum to $$ \sum_{i=1}^3(x_i-\bar x)(x_i-\bar x)^{\mathsf T} = \begin{pmatrix} 1+1+0&-1/3+1/3+0\\ -1/3+1/3+0&1/9+1/9+4/9 \end{pmatrix} = \begin{pmatrix}2&0\\0&2/3\end{pmatrix}. $$ Because $n-1=2$, the centered sample covariance is $$ S=\begin{pmatrix}1&0\\0&1/3\end{pmatrix}. $$ This calculation illustrates the divisor and centering directly. If the observations are viewed as a realization from a Gaussian model, the Wishart theorem describes the distribution across repeated samples; it does not say that this particular matrix equals the population covariance.

Divisors and a Common Pitfall

Some definitions use the centered scatter divided by $n$, rather than $n-1$. If $S_n=n^{-1}\sum_i(X_i-\bar X)(X_i-\bar X)^{\mathsf T}$, then $S_n=((n-1)/n)S$. Its distribution is therefore the same Wishart scatter with the corresponding scale factor, and $E[S_n]=((n-1)/n)\Sigma$. The divisor changes the scale and bias, not the underlying centered scatter.

Common mistake: Treating the $n$ centered observation vectors as independent $N_p(0,\Sigma)$ vectors and concluding that the scatter has $n$ Wishart degrees of freedom. Centered observations are constrained to sum to zero, so they are dependent. The orthogonal transformation in the proof isolates exactly $n-1$ independent Gaussian rows; that is the reason the degrees of freedom are $n-1$.

The Wishart representation and entrywise covariance formulas provide a finite-sample baseline for PCA inference: they describe the randomness in the covariance estimate from which sample principal components are computed. They do not, by themselves, give the distribution of a sample eigenvalue or eigenvector. Matrix entries can fluctuate jointly, and translating that fluctuation into eigenvalue or direction uncertainty requires additional analysis.

Check Your Understanding

Use the definitions and results above to answer each question.

Under the Gaussian sampling model, what are the degrees of freedom in the Wishart representation of the centered scatter matrix, and why?
If $n=6$, what is the exact Wishart distribution of $5S$?
For $p=2$, express $\operatorname{Var}(S_{12})$ in terms of the entries of $\Sigma$.
Why can the sample covariance between two variables be nonzero even when their population covariance is zero?
If the centered scatter is divided by $n$ instead of $n-1$, what is the expectation of the resulting matrix?