PCA From the Singular Value Decomposition

From a Centered Expression Matrix to PCA Coordinates

PCA represents cells using a smaller set of coordinates while retaining directions of variation in the input data. The singular value decomposition (SVD) gives a direct algebraic route to those coordinates and to the corresponding directions across genes. Understanding which SVD factor represents cells, which represents genes, and how the singular values connect them prevents a common source of confusion: treating scores and loadings as if they were interchangeable.

Let \(X\) be a centered matrix with \(n\) cells in rows and \(p\) genes in columns. Here, “centered” means that each gene column has mean zero across the cells. Each entry \(X_{ij}\) is the centered value for cell \(i\) and gene \(j\), on whatever expression scale was chosen before centering. No probability distribution for the entries is required for the matrix identities below. The values must be finite real numbers, and the same \(X\) must be used throughout the calculation.

The columns of \(U\) are cell-space directions, and the columns of \(V\) are gene-space directions. The \(j\)th column \(v_j\) of \(V\) gives one direction in the \(p\)-gene coordinate system. The corresponding \(j\)th column \(u_j\) of \(U\) gives a unit-length pattern across the \(n\) cells. The singular value \(d_j\) sets the scale linking these two directions.

Scores, Loadings, and the SVD

For PCA, the gene-space directions are commonly called loading directions, or simply loadings. The coordinates of cells along those directions are the PCA scores. With the SVD convention above, the score matrix is \(T=XV\). Substituting the SVD shows directly that \(T=UD\). Thus, scores are not \(U\) alone unless the singular values happen to equal one.

The dimensions make the roles concrete. The score matrix \(T\) is \(n\) by \(r\): one row per cell and one column per retained component. The loading matrix \(V\) is \(p\) by \(r\): one row per gene and one column per component. A cell’s score on component \(j\) is a single entry \(T_{ij}\); a gene’s loading on component \(j\) is \(V_{gj}\). The reconstruction identity is \(X=TV^T\), since \(UDV^T=TV^T\).

If only the first \(k\) components are retained, with \(k\le r\), write \(U_k\), \(V_k\), and \(D_k\) for the corresponding first \(k\) columns or singular values. The retained score matrix is \(T_k=U_kD_k=XV_k\), with dimensions \(n\) by \(k\). The retained loadings \(V_k\) have dimensions \(p\) by \(k\), and \(T_kV_k^T\) is an \(n\) by \(p\) matrix on the original centered scale. This identity describes how the retained coordinates map back to gene space; it does not, by itself, prescribe how many components to keep.

Singular Values and the Covariance Matrix

The SVD also explains why its right singular vectors are the PCA loading directions. For the usual sample covariance matrix of the centered genes, \(C=X^TX/(n-1)\), the SVD gives:

To verify the second identity, use \(X=UDV^T\), \(U^TU=I\), and \(V^TV=I\):

Here \(e_j\) is the \(j\)th coordinate vector. In particular, \(v_j\) is an eigenvector of the gene covariance matrix and its eigenvalue is \(d_j^2/(n-1)\). The corresponding score vector has squared length \(t_j^Tt_j=d_j^2\), because \(t_j=d_ju_j\) and \(u_j^Tu_j=1\). These relationships use the sample-covariance denominator \(n-1\); using a different covariance convention changes that denominator, not the SVD itself.

Centering also constrains the possible rank. Since every column of \(X\) sums to zero, its columns lie in the \((n-1)\)-dimensional subspace perpendicular to the all-ones vector. Consequently, \(r\le\min(n-1,p)\). If some genes or cells are linearly dependent after centering, the rank can be smaller still. Components beyond the rank have zero singular value and do not provide additional nonzero score directions.

Truncation, Sign, and Numerical Checks

A truncated SVD keeps \(k\) components and forms \(T_kV_k^T\). This produces a matrix with the same dimensions as \(X\), while \(T_k\) still has one row per cell and \(V_k\) one row per gene. In applications with many more genes than cells, the centered matrix has at most \(n-1\) nonzero components. A computation may therefore report only a selected subset rather than constructing every column of a full SVD.

SVD vectors have a sign ambiguity. If \(u_j\) and \(v_j\) are both multiplied by \(-1\), then \(d_ju_jv_j^T\) is unchanged, as are the reconstructed data and the component’s underlying direction. Scores and loadings for that component both change sign. Therefore, a positive score or loading has no universal meaning without specifying the chosen sign convention and examining the associated gene direction.

If two or more singular values are equal, the individual vectors within their shared subspace are not uniquely determined: different orthonormal bases of that subspace give valid SVDs. Comparisons of such components across software runs or data perturbations should not assume that each basis vector must match. The subspace can be stable even when its particular basis changes.

Useful checks include confirming that the input columns are centered, that score and loading matrices have the expected dimensions, and that the factors reproduce the input to numerical precision. For a compact SVD, verify \(U^TU=I\), \(V^TV=I\), and \(X\) approximately equals \(UDV^T\). Also check that singular values are nonnegative and ordered, and that the number of reported nonzero values does not exceed \(\min(n-1,p)\) for centered data. Small numerical deviations from exact zero or exact orthogonality can arise from floating-point computation.

Choosing a Computational Route and Avoiding Misinterpretation

A direct SVD factors \(X\) and yields scores and loading directions without first forming a covariance matrix. An alternative is to eigendecompose \(X^TX\), a \(p\) by \(p\) matrix: its eigenvectors give the loading directions and its nonnegative eigenvalues are the squared singular values. When \(p\) is large, forming this matrix can be costly; it can also magnify numerical difficulties associated with small singular values. An eigendecomposition of \(XX^T\), an \(n\) by \(n\) matrix, can be more convenient when there are many genes and fewer cells, after which loading directions can be recovered from \(X^Tu_j/d_j\) for nonzero \(d_j\). These are algebraically related routes, but their computational costs differ.

The SVD is exact as a factorization, but the scientific interpretation depends on the input matrix. Centering gene columns defines the reference level for each gene, and prior transformations or scaling determine which differences contribute to the decomposition. A large loading magnitude indicates a substantial coefficient in that component’s gene direction; it does not alone establish a gene’s biological importance, a causal role, or a statistically significant effect. Scores are coordinates relative to the chosen directions, not cell labels.

A frequent mistake is to call the columns of \(U\) the PCA scores and stop there. Those columns are unit-length cell-space vectors; the scores are \(UD\), so each column is multiplied by its singular value. A second mistake is to interpret a loading as a cell coordinate. Check dimensions: rows of the score matrix index cells, and rows of the loading matrix index genes. Keeping these roles distinct makes both computation and interpretation clearer.