Computing Leading PCs With Truncated SVD

Why Compute Only the Leading Components?

For a single-cell expression matrix, the number of cells and genes can make a full singular value decomposition unnecessarily expensive. PCA usually needs only a modest number of leading components, not every singular vector. A truncated SVD aims to compute that leading subspace directly, avoiding the work and storage associated with the rest of the decomposition. This matters especially when a matrix has many rows, many columns, or both.

Let \(A\) be the \(n \times p\) matrix supplied to PCA after the centering and gene-scaling decisions described in the earlier material on gene scaling and PCA. Rows represent cells and columns represent genes. Its singular value decomposition is \(A=U\Sigma V^T\), with singular values in nonincreasing order. The first \(k\) left singular vectors describe the leading cell-space directions, while the corresponding right singular vectors describe gene-space directions. The rank-\(k\) approximation is \(U_k\Sigma_kV_k^T\). The computational task here is to obtain those leading directions efficiently and to check the quality of an approximation.

The Frobenius norm satisfies \(\|A\|_F^2=\sum_{i,j} A_{ij}^2=\sum_i \sigma_i^2\). For the exact leading \(k\)-dimensional left singular subspace, the squared residual is \(\sum_{i>k}\sigma_i^2\). An approximate subspace may leave a larger residual. These criteria assess reconstruction of the supplied matrix; they do not decide how many components should later be used for visualization.

Truncated and Randomized Computation

An exact full SVD computes all singular values and vectors. A truncated iterative method, commonly based on Lanczos bidiagonalization, instead builds approximations to the leading singular vectors through repeated matrix-vector products. It can be useful when only a few components are required. A randomized method follows a different route: it uses random projections to find a relatively small space that captures the important directions, then performs a much smaller SVD within that space.

In a randomized range-finding method, choose a target rank \(k\), an oversampling amount \(s\), and let \(\ell=k+s\). Generate a random matrix \(\Omega\) with \(p\) rows and \(\ell\) columns. Form \(Y=A\Omega\), and compute an orthonormal basis \(Q\) for the columns of \(Y\). The range of \(Q\) is a candidate approximation to the leading left singular subspace. Compute a small SVD of \(B=Q^TA\), then retain the leading \(k\) components of that smaller decomposition. The final output has rank \(k\), even though the search space had dimension \(\ell\).

Oversampling gives the random search space extra dimensions in which to capture directions that might otherwise be missed. It is a safeguard against an unlucky projection, not an increase in the number of components ultimately retained. A larger \(s\) generally increases the cost of matrix multiplications and the storage for \(Q\), so it should not be increased without reason. Values such as 5 to 20 are common starting points, but the appropriate choice depends on the matrix, target accuracy, and computational constraints.

When the singular values decline slowly, a basic random projection can mix important and less important directions. Power iterations can sharpen the separation. For an integer \(q\geq 0\), construct a basis from

A singular direction associated with singular value \(\sigma_i\) is weighted by \(\sigma_i^{2q+1}\) in this transformed matrix. Thus, directions with larger singular values become more dominant relative to smaller ones as \(q\) increases. In practice, orthonormalize during the iterations to avoid numerical instability. Each extra power iteration requires additional passes through the matrix, so it improves subspace separation at a computational cost.

The cost is usually dominated by products involving \(A\) or \(A^T\), rather than the final small SVD. This is one reason randomized methods are useful for large matrices. Centering can make a sparse expression matrix dense if performed by explicitly subtracting gene means. Implementations can often apply the centering operation implicitly in matrix-vector products, preserving the computational advantages of the original sparse representation. The details depend on the software and preprocessing, so verify what matrix operation the implementation actually uses.

Checking Approximation Error

For an orthonormal basis \(Q\), the projection \(QQ^TA\) is the closest matrix to \(A\) with columns in the span of \(Q\), under the Frobenius norm. The residual \(R=A-QQ^TA\) records what that subspace fails to represent. Its relative Frobenius residual is scale-free, making it easier to compare matrices with different overall magnitudes. It is important to distinguish this reconstruction check from a test of whether particular cell groups or biological signals are preserved.

Diagnostics, Alternatives, and Limitations

A small residual is useful evidence that a subspace reconstructs much of the input matrix, but it should be considered alongside other checks. Confirm that the reported basis is close to orthonormal by inspecting \(Q^TQ\); its diagonal should be near one and its off-diagonal entries near zero. For individual approximate singular triplets \((u_i,\sigma_i,v_i)\), check both \(\|Av_i-\sigma_i u_i\|\) and \(\|A^Tu_i-\sigma_i v_i\|\). Large values indicate that the vectors do not closely satisfy the singular-vector equations.

When an exact truncated decomposition is feasible, compare its leading singular values or subspace with the approximate result on a representative matrix. Signs of singular vectors are arbitrary, so compare subspaces or reconstructed projections rather than requiring matching signs. If repeated or nearly repeated singular values occur, individual vectors can rotate substantially within their shared subspace even when the subspace itself is accurate. Repeating a randomized computation with different random seeds can reveal sensitivity; materially different results merit investigation of rank, oversampling, power iterations, or numerical stability.

The principal alternatives have different costs and uses:

A common mistake is to set the search dimension equal to the desired number of components and assume the resulting approximate subspace must be adequate. The correction is to treat \(k\) as the output rank, use oversampling when needed to make subspace capture more reliable, and assess the result with residual and stability checks. Another mistake is to interpret a small global residual as proof that rare cell states or low-variance biological signals have been preserved. The Frobenius norm is dominated by contributions to total squared magnitude; it does not guarantee accuracy for every subset or feature.

Truncated computation also does not settle how many components should be used for a downstream display, nor does it remove the effects of earlier preprocessing choices. It produces a computationally efficient approximation to the PCA of the matrix supplied. The meaning of its components remains tied to that matrix and to the scaling and centering used to construct it.