Answer

SVD and PCA: how the singular value decomposition does the same job

Inspired by a question on Cross Validated · November 4, 2026

pcadimensionality reduction

Two routes to the same answer

Textbooks usually introduce PCA through the covariance matrix: compute it, find its eigenvectors and eigenvalues, and project the data onto the eigenvectors. (The plain-language PCA page explains what those eigenvectors and eigenvalues mean.) The singular value decomposition, or SVD, skips the covariance matrix and works on the data directly. Any n-by-p matrix with n at least p can be written as a product of three pieces:

• U (n by p): one row per observation, with orthonormal columns.
• D (p by p): a diagonal matrix of singular values, sorted from largest to smallest, all zero or positive.
• V (p by p): one row per variable, with orthonormal columns. The decomposition is X = U D V', where V' is the transpose of V.

The link to PCA takes one line of algebra. If X has centered columns, its covariance matrix is X'X / (n - 1). Substituting the SVD gives X'X = V D U'U D V' = V D² V', because U'U is the identity. So the covariance matrix equals V (D² / (n - 1)) V'. That is an eigen-decomposition: the eigenvectors are the columns of V, and the eigenvalues are the squared singular values divided by n - 1.

What each SVD piece means in PCA terms

• V = the principal directions. Each column is a set of weights, one per variable, defining a component. This is what prcomp() calls rotation.
• U D = the component scores. Since X V = U D V'V = U D, each observation's coordinates on the new axes are its row of U scaled by the singular values. U on its own holds the scores rescaled to unit length.
• D² / (n - 1) = the component variances. These are the eigenvalues. prcomp() reports their square roots as sdev, which equal d / sqrt(n - 1).
• Proportion of variance explained is each d² divided by the sum of all d², so the n - 1 cancels.
• Dimensionality reduction is truncation: keep the first k columns of U and V and the first k singular values. The product of these truncated pieces is the closest rank-k matrix to X in the least-squares sense (the Eckart-Young theorem), so dropping the small singular values throws away as little variance as possible.

Centering is essential. An SVD of the raw, uncentered matrix gives a first direction that mostly points towards the column means, which is not PCA. If you want PCA on the correlation matrix instead, also divide each centered column by its standard deviation before the SVD.

Checking the equivalence in R and Python

The example below simulates four variables for 200 observations, centers them, and does PCA both ways: eigen() on the covariance matrix and svd() on the data. It then checks the scores, compares with prcomp(), and builds a rank-2 approximation. The Python tab does the same with NumPy.

set.seed(134282)
n <- 200
# Four measurements driven by one strong and one weaker hidden factor
f1 <- rnorm(n); f2 <- rnorm(n)
X <- cbind(a = 2 * f1 + rnorm(n, sd = 0.5),
           b = 2 * f1 + f2 + rnorm(n, sd = 0.5),
           c = f2 + rnorm(n, sd = 0.5),
           d = rnorm(n))

Xc <- scale(X, center = TRUE, scale = FALSE)  # PCA needs centered columns

# Route 1: eigen-decomposition of the covariance matrix
e <- eigen(cov(Xc))
round(e$values, 4)

# Route 2: SVD of the centered data matrix, Xc = U D V'
s <- svd(Xc)
round(s$d, 4)                   # singular values
round(s$d^2 / (n - 1), 4)       # squared and divided by n - 1: the eigenvalues
round(abs(s$v) - abs(e$vectors), 10)  # V = eigenvectors (up to sign)

# Scores: U D equals Xc V
scores <- s$u %*% diag(s$d)
max(abs(abs(scores) - abs(Xc %*% s$v)))
round(apply(scores, 2, var), 4)  # variances again equal the eigenvalues

# prcomp() itself uses the SVD
p <- prcomp(X)
round(p$sdev^2, 4)
round(s$d^2 / sum(s$d^2), 4)    # proportion of variance explained

# Rank-2 reconstruction keeps the two largest singular values
X2 <- s$u[, 1:2] %*% diag(s$d[1:2]) %*% t(s$v[, 1:2])
round(1 - sum((Xc - X2)^2) / sum(Xc^2), 4)  # share of variance kept

import numpy as np

rng = np.random.default_rng(134282)
n = 200
# Four measurements driven by one strong and one weaker hidden factor
f1 = rng.normal(0, 1, n)
f2 = rng.normal(0, 1, n)
X = np.column_stack([
    2 * f1 + rng.normal(0, 0.5, n),        # a
    2 * f1 + f2 + rng.normal(0, 0.5, n),   # b
    f2 + rng.normal(0, 0.5, n),            # c
    rng.normal(0, 1, n),                   # d
])

Xc = X - X.mean(axis=0)  # PCA needs centered columns

# Route 1: eigen-decomposition of the covariance matrix
values, vectors = np.linalg.eigh(np.cov(Xc, rowvar=False))
order = np.argsort(values)[::-1]          # eigh sorts ascending; flip
values, vectors = values[order], vectors[:, order]
print(np.round(values, 4))

# Route 2: SVD of the centered data matrix, Xc = U D V'
U, d, Vt = np.linalg.svd(Xc, full_matrices=False)
print(np.round(d, 4))                     # singular values
print(np.round(d**2 / (n - 1), 4))        # the eigenvalues
print(np.round(np.abs(Vt.T) - np.abs(vectors), 10))  # V = eigenvectors (up to sign)

# Scores: U D equals Xc V
scores = U * d
print(np.max(np.abs(np.abs(scores) - np.abs(Xc @ Vt.T))))
print(np.round(scores.var(axis=0, ddof=1), 4))
print(np.round(d**2 / np.sum(d**2), 4))   # proportion of variance explained

# Rank-2 reconstruction keeps the two largest singular values
X2 = U[:, :2] * d[:2] @ Vt[:2, :]
print(round(1 - np.sum((Xc - X2)**2) / np.sum(Xc**2), 4))  # share of variance kept

The figures below come from one seeded run of the R code. Python uses a different random number generator, so its figures differ slightly, but the pattern is the same, and every identity holds exactly in both.

• Same eigenvalues both ways. eigen() gives 9.5367, 1.3978, 1.0486 and 0.2182. The singular values are 43.5637, 16.6781, 14.4451 and 6.5901; squaring them and dividing by 199 gives exactly the same four numbers, and so does prcomp().
• Same directions. The absolute values of V and of the covariance eigenvectors agree to 10 decimal places. Only the signs can differ, and a flipped sign changes nothing.
• Same scores. U D and Xc V differ by at most about 1e-14, which is rounding error, and the variances of the scores equal the eigenvalues.
• Reduction. The components explain 78.2%, 11.5%, 8.6% and 1.8% of the variance. Rebuilding the data from the first two singular values keeps 89.6% of the total variance, which is exactly the sum of the first two shares.

Why software prefers the SVD route

• Accuracy. Forming X'X squares the condition number of the problem, so small components can lose most of their digits to rounding. The SVD works on X itself and avoids that. R's prcomp() uses the SVD for this reason, while the older princomp() uses eigen() on the covariance matrix (and divides by n rather than n - 1).
• Wide data. When there are more variables than observations, the p-by-p covariance matrix is huge and singular. The SVD of the n-by-p data matrix is cheaper, and after centering it shows directly that at most n - 1 components have nonzero variance.
• Only a few components needed. Truncated or randomized SVD algorithms compute just the first k singular vectors, which is how PCA is run on very large data sets.
• Loadings. Multiply each column of V by its singular value and divide by sqrt(n - 1) to get loadings in the covariance metric. With standardized variables these are the correlations between variables and components.

How to report PCA variance explained in APA style (7th edition)

Readers care about the components, not the algorithm, so report the analysis as a PCA and mention the SVD only in the method section if at all. Say whether the variables were only centered (covariance matrix) or also standardized (correlation matrix), because the eigenvalues mean different things. Using the example above:

"A principal component analysis of the four centered variables (N = 200), computed by singular value decomposition, showed that the first component (eigenvalue = 9.54) explained 78.2% of the variance and the second (eigenvalue = 1.40) a further 11.5%. The two retained components together explained 89.6% of the variance."

List the weights or loadings for each retained component in a table, with one row per variable, and add rows for the eigenvalues and percentages of variance. A table note should say whether the values are eigenvector weights or loadings.

Related tools and guides

• PCA, eigenvectors and eigenvalues explained
• Factor analysis vs PCA
• Which statistical test should I use?
• Singular value decomposition (Wikipedia)

More answered questions

• Why Euclidean distance breaks down in high dimensions
• Factor analysis vs PCA: what is the difference?
• PCA, eigenvectors and eigenvalues explained in plain language

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.