Principal Component Analysis

Definition Let $x_1, \dots, x_N \in \mathbb{R}^d$ be $N > 1$ observations, each a column of $d$ real measurements. Their mean is the vector

computed coordinate by coordinate. The centered data matrix is the real $N \times d$ matrix $X$ whose $i$th row is $(x_i - \bar{x})^{\mathsf{T}}$. So rows are observations and columns are measured coordinates. If $X_{\text{raw}}$ has rows $x_i^{\mathsf{T}}$ and $\mathbf{1} \in \mathbb{R}^N$ is the vector of ones, then $X = X_{\text{raw}} - \mathbf{1}\bar{x}^{\mathsf{T}}$. Every column of $X$ sums to zero, $\mathbf{1}^{\mathsf{T}}X = 0$, so the rows are dependent and $\operatorname{rank} X \le \min(N - 1, d)$. The mean is the best single center. For every $c \in \mathbb{R}^d$,

because the cross term $2\sum_i (x_i - \bar{x})^{\mathsf{T}}(\bar{x} - c)$ vanishes. So $c = \bar{x}$ is the unique minimizer of the total squared distance.

Definition For a centered data matrix $X \in \mathbb{R}^{N \times d}$ with $N > 1$, the sample covariance matrix is the $d \times d$ matrix

It is real symmetric and positive semidefinite, since $v^{\mathsf{T}}Cv = \lVert Xv \rVert^2/(N - 1) \ge 0$, and $\operatorname{rank} C = \operatorname{rank} X \le N - 1$. The diagonal entry $C_{jj}$ is the sample variance of coordinate $j$, and for $j \ne l$ the entry $C_{jl}$ is the sample covariance of coordinates $j$ and $l$. The divisor $N - 1$ is a convention. Dividing by $N$ instead multiplies $C$ by the positive constant $(N - 1)/N$, which scales every variance by the same factor and leaves every eigenvector and every comparison unchanged.

Definition Let $v \in \mathbb{R}^d$ be a unit vector. The scores of the centered observations along $v$ are the scalar projections $(x_i - \bar{x})^{\mathsf{T}}v$, the entries of $Xv \in \mathbb{R}^N$. They sum to zero, because $\mathbf{1}^{\mathsf{T}}Xv = 0$, so their sample variance is

the projected variance along $v$. This is the Rayleigh quotient of the symmetric matrix $C$ at a unit vector. In this lesson the eigenvalues of $C$ are listed in decreasing order, $\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_d$, to match the order of the singular values. By the bounds on the Rayleigh quotient, the maximum of the projected variance is $\lambda_1$, attained exactly at the unit vectors of the eigenspace of $\lambda_1$. For a $k$-dimensional subspace spanned by the orthonormal columns of $V_k \in \mathbb{R}^{d \times k}$, the retained variance is

the sum of the projected variances along the columns of $V_k$. It depends only on the subspace, and its largest possible value is $\lambda_1 + \cdots + \lambda_k$, as the reconstruction error below shows.

Definition Let $X \in \mathbb{R}^{N \times d}$ be a centered data matrix with $N > 1$ and real full singular value decomposition $X = U\Sigma V^{\mathsf{T}}$, with orthogonal $U$ of size $N$, orthogonal $V$ of size $d$ and singular values $\sigma_1 \ge \cdots \ge \sigma_p \ge 0$, where $p = \min(N, d)$. Since $X^{\mathsf{T}}X = V\Sigma^{\mathsf{T}}\Sigma V^{\mathsf{T}}$, the right singular vectors $v_1, \dots, v_d$ are orthonormal eigenvectors of $C$, with eigenvalues $\sigma_i^2/(N - 1)$ for $i \le p$ and $0$ for any remaining ones. The principal directions, or principal axes, are $v_1, v_2, \dots$ in this order, and the $i$th principal component score vector is $Xv_i \in \mathbb{R}^N$, whose entries are the coordinates of the centered observations along $v_i$. For $i \le p$ the variance along $v_i$ is $\sigma_i^2/(N - 1)$ and $Xv_i = \sigma_iu_i$. The first direction maximizes projected variance, and each later $v_j$ maximizes it among unit vectors orthogonal to $v_1, \dots, v_{j-1}$. The directions live in the measurement space $\mathbb{R}^d$ and the scores in the observation space $\mathbb{R}^N$. Books differ on whether “principal component” names the direction or the score vector, so say which one you mean. The directions are not unique. Replacing $v_i$ by $-v_i$ flips the sign of its scores but leaves the reconstruction $Xv_iv_i^{\mathsf{T}}$ unchanged, and when singular values tie, any orthonormal basis of the tied eigenspace of $C$ serves equally well.

Definition Let $V_k \in \mathbb{R}^{d \times k}$ have orthonormal columns spanning a subspace $S$ of dimension $k$, and let $P = V_kV_k^{\mathsf{T}}$ be the orthogonal projector onto $S$. The reconstruction of the centered data in $S$ is $XP$, whose rows are the orthogonal projections of the centered observations onto $S$. Its reconstruction error is

Each row splits into perpendicular parts, so

The left side does not depend on $S$, and $\lVert XP \rVert_F^2 = \lVert XV_k \rVert_F^2$ is $N - 1$ times the retained variance. So maximizing retained variance and minimizing reconstruction error are the same problem. When $k \le p$ and $S$ is spanned by the top $k$ principal directions, $XP = \sum_{i \le k} \sigma_iu_iv_i^{\mathsf{T}} = X_k$ is the truncated SVD, with error $\sigma_{k+1}^2 + \cdots + \sigma_p^2$. Every $XP$ has rank at most $k$, so the Eckart-Young theorem in the Frobenius norm shows that no $k$-dimensional subspace does better. The largest retained variance is $(\sigma_1^2 + \cdots + \sigma_k^2)/(N - 1)$, and when $k < p$ and $\sigma_k > \sigma_{k+1}$ the best subspace is unique. For the raw observations and an affine subspace $a + S$, the same splitting around the mean gives

so a best affine fit of dimension $k$ can pass through the mean, and $\bar{x} + \operatorname{span}(v_1, \dots, v_k)$ is one.

Definition The answer of principal component analysis depends on two choices made before any computation, the center and the inner product on the measurement space, which the units of the coordinates fix. Without centering, the SVD of $X_{\text{raw}}$ finds best subspaces through the origin, and since $\mathbf{1}^{\mathsf{T}}X = 0$,

A nonzero mean adds a rank-one term that can dominate the leading direction. Centered principal component analysis, by contrast, is unchanged when every observation is translated by the same vector. Changing units multiplies each coordinate $j$ by a positive constant $c_j$, which replaces $X$ by $XD$ and $C$ by $DCD$, where $D = \operatorname{diag}(c_1, \dots, c_d)$, and this generally changes the principal directions. When every coordinate has nonzero variance, standardizing divides each coordinate by its sample standard deviation. That is the choice $c_j = 1/\sqrt{C_{jj}}$, and it turns $C$ into the correlation matrix. Standardizing is a modeling decision with mathematical consequences, not a step the mathematics makes for you. The optimality theorem guarantees a best subspace for the chosen center, units and dimension, and nothing more. A direction of large variance is not automatically important, causal or relevant to a prediction.