PCA with SVD and Reconstruction

Definition Let $X$ be an $n \times p$ data matrix whose rows are the observations $\mathbf{x}_1^{\top}, \dots, \mathbf{x}_n^{\top}$, let $\bar{\mathbf{x}}$ be their mean, and let $X_c$ be the centered matrix whose $i$th row is $(\mathbf{x}_i - \bar{\mathbf{x}})^{\top}$. Put $q = \min(n, p)$. A thin singular value decomposition of $X_c$ is a factorization

in which $U$ is $n \times q$ with orthonormal columns, so $U^{\top}U = I_q$, the matrix $S = \operatorname{diag}(\sigma_1, \dots, \sigma_q)$ is $q \times q$ with $\sigma_1 \ge \dots \ge \sigma_q \ge 0$, and $V$ is $p \times q$ with orthonormal columns, so $V^{\top}V = I_q$. Every real matrix has one. Written column by column, it is a sum of rank-one matrices,

where $\mathbf{u}_j$ and $\mathbf{v}_j$, the $j$th columns of $U$ and $V$, are the left and right singular vectors. In PCA the columns of $V$ are the principal axes, also called loadings. The singular values are determined by $X_c$, but the vectors are not. Each pair $(\mathbf{u}_j, \mathbf{v}_j)$ can change sign together. When a positive singular value is repeated, any orthonormal basis of the space spanned by its right singular vectors can serve, with $\mathbf{u}_j = X_c\mathbf{v}_j / \sigma_j$. The vectors that belong to zero singular values can be any orthonormal vectors that complete the others.

Definition The singular values $\sigma_1 \ge \dots \ge \sigma_q \ge 0$ of $X_c$ are the diagonal entries of $S$. Because $U^{\top}U = I_q$,

so each $\mathbf{v}_j$ is an eigenvector of $X_c^{\top}X_c$ with eigenvalue $\sigma_j^2$. For $n \ge 2$ the sample covariance matrix is $C = X_c^{\top}X_c / (n - 1)$, so $\mathbf{v}_j$ is also an eigenvector of $C$, with eigenvalue

When $p > q$, completing the columns of $V$ to an orthonormal basis of $\mathbb{R}^p$ adds $p - q$ eigenvectors of $C$ with eigenvalue $0$. The rank of $X_c$ is the number of positive singular values, and the squares of the singular values add up to the sum of the squared entries of $X_c$, written $\lVert X_c \rVert_F^2$ for the squared Frobenius norm. The rows of $X_c$ add to the zero vector, so the rank of $X_c$ is at most $\min(n - 1, p)$.

Definition Given the loadings $V$, the principal component score matrix is

an $n \times q$ matrix whose entry $t_{ij} = (\mathbf{x}_i - \bar{\mathbf{x}})^{\top}\mathbf{v}_j$ is the coordinate of observation $i$ along axis $j$. The second equality holds because $V^{\top}V = I_q$. Each score column has mean $0$ and sample variance $\lambda_j = \sigma_j^2/(n - 1)$, and different score columns are uncorrelated, since $T^{\top}T = S^2$ is diagonal. A new observation $\mathbf{x}$ gets the scores $\mathbf{t} = V_k^{\top}(\mathbf{x} - \bar{\mathbf{x}})$ on the first $k$ axes, where $V_k$ holds the first $k$ columns of $V$, and the mean and loadings are the ones fitted on the training data, never refitted on the new observation. A score changes sign when its loading does. On this page a loading is a unit vector. Some texts call $\sqrt{\lambda_j}\, \mathbf{v}_j$ the loading, so check the scaling before comparing numbers.

Definition The variance explained by component $j$ is $\lambda_j = \sigma_j^2/(n - 1)$, the sample variance of its scores. The total variance is the sum of the $p$ feature variances, which equals the trace of $C$ and the sum of all its eigenvalues, with $\lambda_l = 0$ for $l > q$. When $X_c \ne 0$, the explained variance ratio of component $j$ is

and the cumulative ratio of the first $k$ components is the sum of their ratios. When every feature is constant, $X_c = 0$ and the ratio is undefined rather than zero. The ratios depend on the scaling of the features, so PCA on centered raw features and PCA on standardized features give different ratios, and they describe the data the PCA was fitted on.

Definition For $1 \le k \le q$, the rank-$k$ reconstruction of the centered matrix keeps the first $k$ terms of the decomposition,

where $U_k$, $V_k$ and $T_k$ hold the first $k$ columns of $U$, $V$ and $T$, and $S_k = \operatorname{diag}(\sigma_1, \dots, \sigma_k)$. Row $i$ is the orthogonal projection of $\mathbf{x}_i - \bar{\mathbf{x}}$ onto the span of $\mathbf{v}_1, \dots, \mathbf{v}_k$. In the original units the reconstructed observation is

and a new observation is reconstructed by the same formulas with the training mean and loadings. If the features were standardized before PCA, each coordinate of $V_k \mathbf{t}_i$ is multiplied by its training standard deviation before the mean is added. By the Eckart-Young theorem, among all $n \times p$ matrices $B$ of rank at most $k$, the matrix $\hat{X}_{c,k}$ makes $\lVert X_c - B \rVert_F^2$ as small as possible. For $k < q$ this minimizer is unique when $\sigma_k > \sigma_{k+1}$.

Definition The total squared reconstruction error of the rank-$k$ reconstruction is

the sum of the squared singular values that were dropped, which is $0$ for $k = q$. Adding the mean back does not change any residual. Dividing $E_k$ by $n$ gives the mean squared error per observation, and dividing by $np$ gives the mean squared error per matrix entry, so a report should state its denominator. When $X_c \ne 0$, the fraction $E_k / \lVert X_c \rVert_F^2$ equals one minus the cumulative explained variance ratio of the first $k$ components. For a new observation the error is $\lVert \mathbf{x} - \hat{\mathbf{x}} \rVert_2^2$ with the training mean and loadings, and nothing bounds it by the training figures. After standardization the errors are in standardized units unless the scaling is reversed.