PCA and Matrix Factorization

Principal component analysis (PCA) represents observations using directions of large variation. Linear Algebra for Machine Learning introduces covariance, eigenvectors, and singular value decomposition (SVD). Here we connect those objects to PCA scores and reconstruction, then examine scaling, component selection, whitening, and non-negative matrix factorization (NMF). Both PCA and NMF approximate a data matrix with smaller factors, but their constraints give those factors different meanings.

From a data matrix to scores and reconstruction

Let \(X\) have \(n\) rows of observations and \(d\) feature columns. Subtract the feature means to obtain \(X_c=X-\mathbf1\bar x^\top\), where \(\mathbf1\) repeats the mean across rows. Write its SVD as \(X_c=U\Sigma V^\top\), with singular values \(\sigma_j\) in descending order. The columns of \(V\) are unit feature-space directions. The sample covariance is \(X_c^\top X_c/(n-1)=V\Sigma^\top\Sigma V^\top/(n-1)\), so these directions are also covariance eigenvectors, with eigenvalues \(\lambda_j=\sigma_j^2/(n-1)\). Each eigenvalue is the sample variance along its unit direction, so descending eigenvalues order the directions by retained variance. The explained-variance ratio is \(\lambda_j/\sum_l\lambda_l\).

Keep the first \(k\) directions in \(V_k\), a \(d\)-by-\(k\) matrix. The scores \(T_k=X_cV_k=U_k\Sigma_k\) give each observation’s \(k\) coordinates. Reconstruction is \(\hat X=T_kV_k^\top+\mathbf1\bar x^\top\). Thus the centered data is approximated by an \(n\)-by-\(k\) matrix times a \(k\)-by-\(d\) matrix: this is its low-rank factorization. In scikit-learn, components_ stores \(V_k^\top\), transform produces the scores, and inverse_transform restores feature coordinates and the fitted mean.

For a small calculation, take centered rows \((-2,-1),(0,0),(2,1)\). Their only varying direction is \(v=(2,1)/\sqrt5\). Dotting each row with \(v\) gives scores \(-\sqrt5,0,\sqrt5\); multiplying each score by \(v\) restores its two coordinates exactly. One number per row suffices because all three rows lie on one line. More generally, tied eigenvalues allow multiple equally good PCA bases, and each component’s overall sign is arbitrary.

Two definitions, one answer

For centered data and an orthogonal projection onto a \(k\)-dimensional subspace, maximizing retained variance and minimizing squared reconstruction error select the same subspace. Each observation splits into a projection and a perpendicular residual. By the Pythagorean identity, their squared lengths sum to the original centered squared length. Summing over observations gives \(\lVert X_c\rVert_F^2=\lVert T_k\rVert_F^2+\lVert X_c-T_kV_k^\top\rVert_F^2\), where the squared Frobenius norm is the sum of squared matrix entries.

import numpy as np
from sklearn.decomposition import PCA

rng = np.random.default_rng(0)
X = rng.normal(size=(500, 6)) @ rng.normal(size=(6, 6))
X = X - X.mean(0)
total = np.sum(X ** 2)
p = PCA(svd_solver="full").fit(X)

print(f"{'k':>3} {'explained ratio':>17} {'cumulative':>11} {'recon err frac':>15}")
for k in range(1, 7):
    fitted = PCA(k, svd_solver="full").fit(X)
    recon = fitted.inverse_transform(fitted.transform(X))
    print(f"{k:3d} {p.explained_variance_ratio_[k-1]:17.6f}"
          f" {p.explained_variance_ratio_[:k].sum():11.6f}"
          f" {np.sum((X - recon) ** 2) / total:15.6f}")
#   k   explained ratio  cumulative  recon err frac
#   1          0.500593    0.500593        0.499407
#   2          0.271689    0.772282        0.227718
#   3          0.158166    0.930447        0.069553
#   4          0.045756    0.976204        0.023796
#   5          0.023747    0.999951        0.000049
#   6          0.000049    1.000000        0.000000

The cumulative-variance fraction and reconstruction-error fraction add to one in exact arithmetic; displayed values can differ by rounding. For \(k=2\), the table gives \(0.772282+0.227718=1\). “Keeps 90% of variance” means that the training reconstruction error is 10% of the total centered sum of squares. It does not mean 10% error in each feature, a 10% relative error in each observation, or the same error on new data.

The truncated SVD attains the smallest training squared error among rank-at-most-\(k\) approximations to \(X_c\); its residual sum of squares is \(\sum_{j>k}\sigma_j^2\). Restoring the mean gives the best affine \(k\)-dimensional reconstruction of the original data under this loss. This optimization uses the inputs alone and does not guarantee a good prediction of a separate target. It can also be dominated by outliers because large residuals are squared.

Scaling decides the components

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(1)
n = 800
X = np.column_stack([rng.normal(0, 1, n),
                     rng.normal(0, 1, n) * 1000,     # independent draw, scaled by 1000
                     rng.normal(0, 1, n)])

for label, M in (("raw", X), ("standardized", StandardScaler().fit_transform(X))):
    p = PCA().fit(M)
    print(f"{label:14s} explained {np.round(p.explained_variance_ratio_, 4)}")
    print(f"{'':14s} PC1 abs coefficients {np.round(np.abs(p.components_[0]), 4)}")
p_raw = PCA().fit(X)
print(f"raw PC1 ratio before rounding: {p_raw.explained_variance_ratio_[0]:.9f}")
# raw            explained [1. 0. 0.]
#                PC1 abs coefficients [0. 1. 0.]
# standardized   explained [0.35   0.3437 0.3063]
#                PC1 abs coefficients [0.5936 0.2385 0.7686]
# raw PC1 ratio before rounding: 0.999998187

The raw first-component variance fraction is 0.999998187, which rounds to 1.0000. Its direction is nearly the second coordinate axis, not exactly equal to it. The fit retains all three components; projecting onto only the first would discard the other directions. The columns are independent Gaussian draws whose population standard deviations are 1, 1000, and 1, so the raw variance objective is dominated by the second column’s scale.

After standardization, the variance fractions are about 0.35, 0.34, and 0.31. The population covariance of the standardized variables is the identity, which has no preferred direction: any orthonormal basis is a population PCA basis. The particular mixed directions printed here come from finite-sample correlations. The absolute coefficients show how strongly each feature contributes to a direction while suppressing the arbitrary overall sign.

Standardization measures reconstruction error relative to each feature’s sample standard deviation. It is useful when raw units would overwhelm the comparison, but it can also give a noisy, low-variance feature more influence. Shared physical units do not settle the choice, and low variance does not imply little predictive information. Choose the scaling to reflect the intended error measure and available measurement knowledge, then evaluate it for the task.

How many components

import numpy as np
from sklearn.datasets import make_classification
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

X, y = make_classification(n_samples=2000, n_features=40, n_informative=5,
                           n_redundant=15, random_state=0)
cum = np.cumsum(PCA().fit(StandardScaler().fit_transform(X)).explained_variance_ratio_)

for k in (2, 5, 10, 20, 40):
    score = cross_val_score(Pipeline([("s", StandardScaler()), ("p", PCA(k)),
                                      ("m", LogisticRegression(max_iter=2000))]),
                            X, y, cv=5).mean()
    print(f"k={k:3d}  cumulative variance {cum[k-1]:.4f}   cv accuracy {score:.4f}")
# k=  2  cumulative variance 0.3398   cv accuracy 0.7535
# k=  5  cumulative variance 0.5045   cv accuracy 0.8390
# k= 10  cumulative variance 0.6456   cv accuracy 0.8450
# k= 20  cumulative variance 0.8932   cv accuracy 0.8450
# k= 40  cumulative variance 1.0000   cv accuracy 0.8475

Ten components retain 64.56% of the variance and yield mean cross-validation accuracy 0.8450. Twenty retain 89.32% and give the same mean accuracy in this run. That is no observed gain for this classifier and these folds; it does not establish that the added directions contain no useful information. The variance fractions describe a separate full-data fit, while each scored pipeline fits its scaler and PCA using only that fold’s training rows.

A variance threshold can be a valid choice for a specified reconstruction-error budget. For prediction, compare candidate component counts using the downstream metric, including an appropriate no-reduction baseline. Scaling, centering, and PCA must all be fitted inside each training fold and then applied unchanged to validation rows. If these scores select \(k\), use a fresh test set or nested cross-validation to estimate the selected procedure’s performance. A curve fitted to all observations should not guide choices before an evaluation that claims those observations were held out.

A scree plot displays ordered eigenvalues and looks for a bend; the bend need not be clear or unique. Parallel analysis compares observed eigenvalues with a reference distribution from simulated null data or repeated independent permutations of each feature column. The latter preserves each column’s marginal values while breaking its associations with the others. A chosen reference quantile provides a retention threshold. This asks which directions exceed that null model’s variation; it does not directly minimize reconstruction error or predict the best supervised component count.

Whitening

Whitening rescales retained PCA scores by \(z_j=t_j/\sqrt{\lambda_j}\), using the fitted component variance \(\lambda_j\). For positive retained eigenvalues, the whitened training scores have identity sample covariance under the same variance convention. This means unit variance and zero cross-covariance; it does not imply independent coordinates or an identity covariance on new data. With the fitted scales retained, whitening is invertible within the retained subspace.

Whitening changes the metric seen by a distance-based or penalized model and can improve numerical conditioning. Whether it improves predictive performance depends on the data and downstream model. Low-variance directions can contain noise or useful signal, so treating all retained directions equally is a modeling choice.

Small eigenvalues require care: their training scores shrink along with their standard deviations, so whitening does not automatically create enormous training values. A new displacement along such a direction can be strongly amplified. Exactly zero-variance directions cannot be whitened by this formula; near-zero directions call for a conditioning check, truncation, or a regularized scaling rule. There is no universally safe count such as ten components.

import numpy as np
from sklearn.decomposition import PCA

rng = np.random.default_rng(4)
X_white = rng.normal(size=(1000, 2)) * np.array([1.0, 0.001])
p_white = PCA(2, whiten=True, svd_solver="full").fit(X_white)
Z = p_white.transform(X_white)
small_sd = np.sqrt(p_white.explained_variance_[-1])
probe = p_white.mean_ + 0.01 * p_white.components_[-1]
print("training covariance after whitening:", np.round(np.cov(Z, rowvar=False), 6))
print(f"last component standard deviation: {small_sd:.6f}")
print(f"new displacement 0.01 along that direction: {p_white.transform(probe[None, :])[0, -1]:.4f}")
# training covariance after whitening: [[ 1. -0.]
#  [-0.  1.]]
# last component standard deviation: 0.000983
# new displacement 0.01 along that direction: 10.1720

The training covariance is identity to the printed precision, even though the second component’s standard deviation is only about 0.000983. A new displacement of 0.01 along that direction becomes about 10.17 whitened units. This is sensitivity relative to the fitted scale, not an automatic failure: whether such a displacement is signal, noise, or distribution shift depends on the application.

Non-negative matrix factorization

NMF approximates a non-negative data matrix by \(X\approx WH\), with \(W\in\mathbb R_{\geq0}^{n\times k}\) and \(H\in\mathbb R_{\geq0}^{k\times d}\). A row of \(H\) is a component pattern; a row of \(W\) gives the amounts used to build an observation. For example, two patterns \((1,0,1)\) and \((0,2,0)\), used in amounts 2 and 1, reconstruct \(2(1,0,1)+(0,2,0)=(2,2,2)\). Negative cancellation is unavailable.

The squared-error version minimizes \(\tfrac12\lVert X-WH\rVert_F^2\) subject to non-negativity. Numerical methods alternate updates of the two factors; the joint problem is non-convex, so a fitted solution need not be globally optimal. Unlike ordinary PCA, this formulation works on the non-negative levels without subtracting the feature means. The constraint can support an additive interpretation, but does not guarantee sparse, unique, or meaningful parts. Different losses, such as generalized KL divergence, express different reconstruction objectives; the NMF documentation specifies the supported losses and factor conventions. Exercise 1 measures how well additive parts are separated in one constructed example.

Factorization variants

Explicitly subtracting feature means usually fills a sparse matrix’s zeros, but PCA does not always have to materialize that dense centered matrix. Supported sparse solvers can account for centering implicitly, and a covariance-based route instead forms a potentially large feature-by-feature matrix. TruncatedSVD deliberately optimizes an uncentered approximation and is commonly used for latent semantic analysis of sparse text matrices. That changes the objective: its leading direction can reflect the mean level as well as variation. Choose between centered and uncentered models based on that meaning and the available solver, not a universal rule that sparse input requires one method.

How many components, on real data

The handwritten digit images bundled with scikit-learn contain 8 by 8 pixels, so each image has 64 features. The following fit treats those pixel intensities on their original common scale. It describes reconstruction of this dataset; no digit labels are used to choose the directions.

import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA

X, y = load_digits(return_X_y=True)
print("data shape:", X.shape)

ev = np.cumsum(PCA().fit(X).explained_variance_ratio_)
for k in (2, 10, 21, 41, 64):
    print(f"  {k:2d} components  cumulative variance {ev[k - 1]:.4f}")
for target in (0.80, 0.90, 0.95, 0.99):
    print(f"  reach {target:.0%}: {int(np.searchsorted(ev, target)) + 1} components")
# data shape: (1797, 64)
#    2 components  cumulative variance 0.2851
#   10 components  cumulative variance 0.7382
#   21 components  cumulative variance 0.9032
#   41 components  cumulative variance 0.9901
#   64 components  cumulative variance 1.0000
#   reach 80%: 13 components
#   reach 90%: 21 components
#   reach 95%: 29 components
#   reach 99%: 41 components

Two components retain 28.51% of this dataset’s centered variation. A two-dimensional scatter plot can still reveal some relationships, but distances and overlap in that plot omit the remaining 71.49%. Apparent overlap in the plot does not establish that the images are inseparable using more features.

Retaining 80%, 90%, 95%, and 99% requires 13, 21, 29, and 41 components respectively. These values quantify a compression tradeoff; the few printed thresholds do not establish that the complete eigenvalue curve has no elbow. A visualization, an image reconstruction, and a classifier may justify different choices of \(k\).

Exercises

1. Comparing additive factors with PCA directions. Generate four disjoint non-negative parts. Compare component concentration and reconstruction error for PCA, uncentered truncated SVD, and NMF.

A component’s peak block identifies where it is strongest, but does not establish exact recovery. Measure how much of its total mass lies in that block.

import numpy as np
from sklearn.decomposition import PCA, NMF

rng = np.random.default_rng(0)
parts = np.zeros((4, 20))
for i in range(4):
    parts[i, i * 5:(i + 1) * 5] = 1          # four disjoint blocks of five
X = rng.random((600, 4)) @ parts + 0.05 * rng.random((600, 20))

pca = PCA(4).fit(X)
nmf = NMF(4, init="nndsvda", max_iter=1000, random_state=0).fit(X)

peak = lambda C: [int(np.abs(c).reshape(4, 5).sum(1).argmax()) for c in C]
print(f"PCA components: negative entries {np.mean(pca.components_ < 0):.2%}"
      f"   range [{pca.components_.min():.3f}, {pca.components_.max():.3f}]")
print(f"NMF components: negative entries {np.mean(nmf.components_ < 0):.2%}"
      f"   range [{nmf.components_.min():.3f}, {nmf.components_.max():.3f}]")
print(f"NMF component -> block it peaks in: {peak(nmf.components_)}")
print(f"PCA component -> block it peaks in: {peak(pca.components_)}")
block_mass = nmf.components_.reshape(4, 4, 5).sum(axis=2)
purity = block_mass.max(axis=1) / block_mass.sum(axis=1)
print("NMF fraction of component mass in its peak block:", np.round(purity, 4))
U, singular, Vt = np.linalg.svd(X, full_matrices=False)
recons = {"centered PCA": pca.inverse_transform(pca.transform(X)),
          "uncentered SVD": (U[:, :4] * singular[:4]) @ Vt[:4],
          "NMF": nmf.transform(X) @ nmf.components_}
for name, reconstructed in recons.items():
    print(f"{name:14s} RMSE {np.sqrt(np.mean((X-reconstructed)**2)):.6f}")
# PCA components: negative entries 43.75%   range [-0.257, 0.364]
# NMF components: negative entries 0.00%   range [0.000, 3.928]
# NMF component -> block it peaks in: [0, 1, 2, 3]
# PCA component -> block it peaks in: [0, 2, 3, 1]
# NMF fraction of component mass in its peak block: [0.9472 0.9621 0.9241 0.9888]
# centered PCA   RMSE 0.012885
# uncentered SVD RMSE 0.012898
# NMF            RMSE 0.012922

2. A low-variance predictive direction. Build data where the label depends on a small-variance feature. Compare PCA with a supervised projection, keeping both inside the validation pipeline.

Compare the one-component results with the all-feature baseline. Check what standardization fixes and what it cannot tell an unsupervised method about the labels.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(1)
n = 2000
loud = rng.normal(0, 10, (n, 1))               # huge variance, no signal
signal = rng.normal(0, 1, (n, 1))              # small variance, all the signal
y = (signal.ravel() > 0).astype(int)
X = np.hstack([loud, signal, rng.normal(0, 3, (n, 8))])

lda = cross_val_score(Pipeline([("l", LDA(n_components=1)),
                                ("m", LogisticRegression(max_iter=2000))]),
                      X, y, cv=5).mean()
for k in (1, 2, 3, 5, 10):
    raw = cross_val_score(Pipeline([("p", PCA(k)),
                                    ("m", LogisticRegression(max_iter=2000))]),
                          X, y, cv=5).mean()
    scaled = cross_val_score(Pipeline([("s", StandardScaler()), ("p", PCA(k)),
                                       ("m", LogisticRegression(max_iter=2000))]),
                             X, y, cv=5).mean()
    extra = f"   LDA(1) {lda:.4f}" if k == 1 else ""
    print(f"k={k:3d}  PCA raw {raw:.4f}   PCA scaled {scaled:.4f}{extra}")
print(f"no reduction, all 10 features: "
      f"{cross_val_score(LogisticRegression(max_iter=2000), X, y, cv=5).mean():.4f}")
# k=  1  PCA raw 0.4995   PCA scaled 0.5645   LDA(1) 0.9840
# k=  2  PCA raw 0.4805   PCA scaled 0.6480
# k=  3  PCA raw 0.5025   PCA scaled 0.6865
# k=  5  PCA raw 0.4985   PCA scaled 0.7145
# k= 10  PCA raw 0.9900   PCA scaled 0.9900
# no reduction, all 10 features: 0.9900

3. Randomized SVD for a large matrix. Compare full SVD with two randomized settings. Measure fit time, reconstruction error, and agreement of the retained directions and subspaces.

Run the same matrix through each solver. Better approximation may take more work; measure that tradeoff instead of assuming the leading vectors will match.

import time
import numpy as np
from sklearn.decomposition import PCA
from threadpoolctl import threadpool_limits

rng = np.random.default_rng(0)
X = rng.normal(size=(40_000, 400)) @ rng.normal(size=(400, 400))
configs = [("full", dict(svd_solver="full")),
           ("randomized", dict(svd_solver="randomized", random_state=0)),
           ("refined", dict(svd_solver="randomized", random_state=0,
                            n_oversamples=30, iterated_power=10))]
fits = {}
with threadpool_limits(limits=1):
    PCA(10, svd_solver="full").fit(X[:1000])
    for name, params in configs:
        times = []
        for _ in range(3):
            start = time.perf_counter()
            p = PCA(n_components=10, **params).fit(X)
            times.append(time.perf_counter() - start)
        fits[name] = p
        print(f"{name:12s} median fit {np.median(times):.2f} s")
    exact = fits["full"]
    Xc = X - exact.mean_
    total = np.sum(Xc**2)
    for name, fitted in fits.items():
        scores = fitted.transform(X)
        error_fraction = (total - np.sum(scores**2)) / total
        print(f"{name:12s} reconstruction error fraction {error_fraction:.6f}")
    for name in ("randomized", "refined"):
        candidate = fits[name]
        align = np.abs(np.sum(exact.components_ * candidate.components_, axis=1))
        subspace_cosines = np.linalg.svd(exact.components_ @ candidate.components_.T,
                                         compute_uv=False)
        print(f"{name:12s} corresponding |cosines| {np.round(align, 4)}")
        print(f"{name:12s} minimum subspace cosine {subspace_cosines.min():.4f}")
# full         median fit 1.28 s # varies by machine
# randomized   median fit 0.63 s # varies by machine
# refined      median fit 1.26 s # varies by machine
# full         reconstruction error fraction 0.907986
# randomized   reconstruction error fraction 0.909193
# refined      reconstruction error fraction 0.908008
# randomized   corresponding |cosines| [0.9886 0.9687 0.931  0.3644 0.1115 0.324  0.692  0.7089 0.8125 0.2713]
# randomized   minimum subspace cosine 0.3222
# refined      corresponding |cosines| [1.     1.     0.9998 0.9998 0.9998 0.9998 0.9992 0.9986 0.9982 0.993 ]
# refined      minimum subspace cosine 0.9922

References

• Lee, D. D., & Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature.
• Halko, N., Martinsson, P.-G., & Tropp, J. A. (2011). Finding Structure with Randomness. SIAM Review.