What it actually does
Data with many correlated columns does not really occupy all those dimensions — it sits near a lower-dimensional plane inside them. PCA finds that plane: the direction along which the points spread out most, then the direction of most remaining spread perpendicular to it, and so on. Keep the first few and you have a compressed description that loses very little. The catch is in what "loses very little" means, and it is worth being suspicious of: PCA measures loss as squared distance, so it keeps whatever varies most. Variance is not the same as information, and the direction that separates two classes can easily be one of the ones you just discarded.
- The interesting structure is linear — PCA finds a flat subspace, and a curved manifold is fitted badly no matter how many components are kept
- Variance is a reasonable proxy for information, which fails whenever the discriminative direction is a low-variance one
- The data is centred, and scaled comparably if the features are in different units
- Components are orthogonal, which is a constraint imposed for convenience rather than a property real latent factors have
- Outliers are absent or handled, since squared error lets a single extreme point define a component
Reconstruction objective
Find the k-dimensional plane that the data sits closest to, measuring closeness as total squared distance. The second form says the same thing from the other side — minimizing what is lost is identical to maximizing the variance that is kept — and the equivalence is exact, not an approximation. Two consequences follow that are worth holding onto. Squared error means outliers dominate the fit, since a single distant point contributes more than a hundred typical ones. And the orthogonality constraint is an assumption rather than a finding: real latent factors are frequently correlated, and PCA will not return them, it will return an orthogonal rotation of them.
- the d-by-k matrix of component directions, constrained to be orthonormal
- the data projected onto the subspace and mapped back — the reconstruction
- squared Frobenius norm: total squared reconstruction error over every entry
- the covariance of the CENTRED data; skipping the centring silently changes the problem
How it is minimized
Closed form via the singular value decomposition of the centred data matrix
- right singular vectors — the component directions, in the original feature space
- singular values; their squares over n-1 are the variances explained by each component
- the scores — the data expressed in the new coordinates, obtained without ever forming XW
- how many components to keep, which the method itself cannot tell you
The problem is a constrained quadratic maximization whose solution is exactly the top eigenvectors of the covariance, so it has a closed form and no iteration. The only real decision is how to compute it, and the answer is the same one the regression entries reach: factorize the data, never the covariance. Forming X-transpose-X and taking its eigendecomposition squares the condition number and throws away roughly half the available precision, which matters here because the components you care about least — the small ones — are precisely where the precision is lost, and those are the ones a reconstruction-error detector depends on. An SVD of the centred matrix gives the identical answer with the conditioning intact. Power iteration remains useful when only the top one or two components are wanted, and randomized methods when d is large.
k (n_components)
How many components to keep. Not learnable — the reconstruction error falls monotonically with k — so it comes from an explained-variance target, a scree plot, or whatever the downstream model needs (enough for 90-99% explained variance, or chosen by downstream cross-validation)
centring
Mandatory, and part of the model. Skipping it makes the first component point at the mean rather than at the direction of greatest spread
scaling
Standardize when features are in different units, or the component structure is decided by whichever column has the largest numbers. Covariance PCA on raw units versus correlation PCA on standardized ones are different analyses
whitening
Divide scores by their singular values so every component has unit variance. Useful before a distance-based model, and it amplifies the noisiest components, which is exactly the trade
solver
Full SVD, randomized SVD, or power iteration — a cost decision, though randomized methods are approximate on a slowly decaying spectrum
Deterministic and exact: there is no iteration, no seed, and no local optimum, which makes it one of the few things here that gives the same answer every time. The instabilities are elsewhere. Signs are arbitrary — a component and its negation are equally valid solutions — so scores flip between library versions and refits unless a convention is pinned, and anything downstream that keyed on the sign breaks silently. Components with near-equal singular values are not individually identified either; their span is stable but the rotation within it is not, so interpreting either one separately is unfounded. And the whole procedure is dominated by outliers, because squared error is: a single extreme point can define the first component on its own.
Full SVD: O(n·d·min(n, d)) time, O(n·d) memory. Covariance eigendecomposition: O(n·d^2 + d^3), cheaper when n is enormous and d is modest, at a real cost in conditioning. Randomized SVD: O(n·d·k) for k components, which is what makes PCA feasible on wide data. Transform: O(d·k) per sample — a single small matrix product, so inference is essentially free.
How, where, and why it is used
Time-Series Forecasting
AdaptedNot a forecaster — a compression step in front of one. With hundreds or thousands of correlated series, PCA over the cross-section extracts a handful of common factors, and those factors are forecast instead of the individual series, then projected back. That is the dynamic factor model in its simplest form, and it is how forecasting problems with more series than observations are made tractable at all.
- Yield-curve modelling, where three components — level, slope, curvature — reproduce almost all the movement across every maturity
- Many-series demand forecasting, where a few latent factors drive most of the co-movement
- Compressing a wide sensor array before a forecaster, so that the model sees ten factors rather than a thousand correlated channels
- Denoising: reconstruct from the top components and treat the discarded remainder as noise
It attacks the specific problem that makes wide forecasting hard — the series are highly correlated, so treating them independently wastes data and treating them jointly is infeasible. Reducing to a few factors first makes both affordable and often more accurate, because the factors are less noisy than any individual series. The limits are structural: PCA is linear and static, so it captures contemporaneous co-movement and nothing about lead-lag relationships, and the components can rotate between refits, which makes a factor a moving target rather than a stable quantity.
- Difference or detrend before the decomposition, or the first component simply captures the shared trend and tells you nothing else
- Standardize across series when they differ in scale, or the largest-magnitude series dominates every component
- Fit the decomposition on the training window only; refitting it on the full series is a leak that is easy to miss because it feels like preprocessing
- Pin the component signs and ordering across refits, or a downstream forecaster sees factors that silently flip
Rolling-origin backtesting of the final forecast, not the explained variance of the decomposition. A factorization that explains 95% of variance and does not improve forecast error was a description of the data, not a useful compression.
- Fitting the components on the whole series before splitting, which leaks future co-movement into the training window
- Assuming components are stable through time when they rotate as the correlation structure shifts
- Discarding low-variance components that carried the idiosyncratic movement a particular series actually needed
Anomaly Detection
Primary fitFit the subspace on normal data and score each new point by how far it falls from it — the squared reconstruction error, known in process control as the squared prediction error or Q statistic. The complementary score is Hotelling’s T-squared, which measures how extreme the point is *within* the subspace. The two answer different questions: T-squared says the point is an unusual amount of normal behaviour, Q says the point is not normal behaviour at all.
- Multivariate statistical process control in manufacturing and chemical plants, where this pairing has been the standard for decades
- Sensor-array monitoring, where a fault breaks the correlation structure between channels before any single channel leaves its range
- Network and telemetry monitoring, where the residual after removing the dominant traffic pattern is where the incident shows
- Screening inputs to a deployed model for falling off the training manifold
It detects exactly the anomaly a per-channel threshold cannot: a point that is entirely in range on every individual feature but violates the relationships between them. That is the common shape of a real multivariate fault, and it is why this remains standard in process control. It is also cheap, deterministic, and explainable — the contribution of each feature to the residual points directly at the channel that broke. The limits follow from linearity: if normal behaviour is a curved manifold, a linear subspace fits it badly and the residual is large everywhere, and a kernel or autoencoder version is the honest answer.
- Standardize on a confirmed-clean window and persist those statistics; the subspace is defined in scaled space
- Fit on clean data only — squared error means a handful of anomalies in the training set can define a component outright
- Keep enough components that the residual is genuinely noise, and few enough that the fault still lands outside the subspace; that balance is the detector
- Track both Q and T-squared, since a fault that stays inside the subspace only shows in the second
Precision@k against confirmed incidents with control limits set from the residual distribution on a clean period, and PR-AUC rather than ROC-AUC. Check the two statistics separately — a detector where only T-squared ever fires is measuring magnitude, not structure.
- Contaminated training data giving an anomaly its own component, after which that anomaly reconstructs perfectly and scores as normal
- Choosing k by explained variance rather than by detection performance, which is a different objective and frequently picks a worse k
- Non-stationary normal behaviour drifting the subspace until the residual grows everywhere and the detector alarms constantly
Optimization
AdaptedPCA is a constrained optimization problem with a closed-form solution — maximize a quadratic form over the orthonormal matrices — and the fact that the answer is an eigendecomposition is the canonical example of a whole class of problems whose optimum is spectral rather than iterative. In applied pipelines it more often plays the upstream role, reducing a high-dimensional decision space before an optimizer has to search it.
- The reference example of a trace maximization under an orthogonality constraint, where the optimum is the top eigenvectors
- Reducing the dimension of a search space before Bayesian optimization, whose surrogate degrades badly above roughly twenty dimensions
- Decorrelating decision variables so that a downstream solver sees a better-conditioned problem
Worth studying here because it is the clearest case where recognizing the structure of a problem replaces solving it: nobody runs gradient ascent on the Rayleigh quotient, because the constraint set and the objective together hand you the answer. That recognition — this is an eigenproblem, so stop iterating — transfers to spectral clustering, to kernel methods, and to a good deal of numerical linear algebra. As a preprocessing step for search it is genuinely useful and genuinely risky, because a direction discarded for low variance may be the one the objective actually depends on.
- Standardize before decomposing, since the components are otherwise decided by units rather than by structure
- Verify that the retained subspace still contains the variation the objective responds to — variance and relevance are different criteria
Compare the closed-form solution against an iterative eigensolver on the same matrix; agreement to numerical tolerance is the correctness check. For dimension reduction ahead of a search, evaluate the search outcome, not the explained variance.
- Reducing a search space along variance and discarding the direction the objective is most sensitive to
- Forming the covariance matrix to get the eigenvectors, which squares the condition number and damages precisely the small components
- Interpreting individual components when their singular values are close, where only the span is identified and not the rotation within it
Computer Vision
AdaptedHistorically, eigenfaces: treat each image as a vector, take the top components of a face dataset, and represent every face by a few dozen coefficients — recognition then becomes nearest neighbour in that space. Currently, the role is smaller and still real: whitening and decorrelating features, compressing learned embeddings before an index, and as the diagnostic that shows how much of a representation’s variance a handful of directions carries.
- Whitening feature vectors before a distance-based retrieval index, where correlated dimensions distort the metric
- Compressing embeddings ahead of an approximate nearest-neighbour index, trading a little recall for a large memory saving
- Eigenfaces and related classical pipelines, which remain the standard teaching example even though learned features superseded them
On raw pixels it is genuinely obsolete, and the reason is instructive rather than merely historical: PCA is linear, so it cannot represent the fact that a shifted image is the same image, and it spends its components on lighting and pose. That failure is exactly the gap the convolutional prior fills. What survives is the operation rather than the model — as a cheap, deterministic decorrelation and compression step on top of learned features, it is still routinely the right tool.
- Centre over the dataset, not per image, or you remove exactly the between-image variation being modelled
- L2-normalize embeddings before decomposing when the downstream metric is cosine
- Retain enough components that retrieval recall is preserved — measure it rather than choosing by explained variance
Downstream recall@k after compression against the uncompressed baseline, plus the memory saved. Explained variance is not the metric that matters here.
- Compressing embeddings so far that the index becomes fast and wrong, with the loss invisible in explained variance
- Re-fitting the decomposition when the encoder is retrained without re-tuning the downstream thresholds
- Treating components of an image dataset as interpretable features when they mostly encode illumination
Natural Language
ViableLatent Semantic Analysis: build a term-document matrix, take its truncated SVD, and represent both documents and terms in a few hundred latent dimensions. Words that co-occur end up close together, which is what lets the representation match a query to a document that never used the query’s words. It is done without centring — a sparse term-document matrix cannot be centred without destroying the sparsity — which technically makes it truncated SVD rather than PCA proper.
- Semantic search and document similarity in classical information retrieval
- Topic-adjacent exploration, where components approximate themes without a probabilistic topic model
- Dimensionality reduction of tf-idf features before a classifier that cannot handle 50,000 sparse columns
It is the historical answer to vocabulary mismatch and still a sensible baseline when a corpus is small, labels are absent, and a transformer is not justified. Its ceiling is the bag-of-words representation underneath it: no order, no negation, no polysemy — a single component cannot hold two senses of a word, which is precisely what contextual embeddings solved. Reach for it when the corpus is modest and the deployment must stay simple.
- tf-idf weighting before the decomposition, or common words dominate the leading components
- Skip centring and use truncated SVD on the sparse matrix directly; centring densifies it and makes the problem infeasible
- A few hundred components is the usual range — far fewer loses distinctions, far more reintroduces the noise that was being removed
Retrieval metrics on held-out query-document pairs, or downstream classifier accuracy against the uncompressed sparse baseline. Reconstruction error says nothing about whether the semantics survived.
- Densifying the term-document matrix in order to centre it, which turns a tractable problem into an impossible one
- Reading components as topics; they are orthogonal directions of variance, and orthogonality is not a property topics have
- Refitting the decomposition on a corpus with a changed vocabulary, which reindexes every dimension
What it costs to run
A single factorization — seconds on a million rows and a few hundred columns, and randomized SVD extends that to very wide data by computing only the components asked for. There is nothing to tune and no restarts, so the cost is exactly one pass.
A centring subtraction and one small matrix product: O(d·k) per sample, microseconds, and the model is a d-by-k matrix plus a mean vector. Small enough to embed anywhere, including inside a database function.
Periodic, driven by drift in the correlation structure rather than by cost. The important operational fact is that a refit changes the coordinate system: component signs flip, near-degenerate components rotate, and anything downstream trained on the old scores is now reading different axes.
- Track the explained-variance profile across refits — a flattening spectrum means the correlation structure that justified the compression is dissolving
- Monitor mean reconstruction error on new data; rising error means the subspace no longer describes the population, and it moves before downstream accuracy does
- Watch the angle between old and new component subspaces rather than comparing components one by one, which is meaningless under rotation
- Alert on scores far outside the training range, since the projection extrapolates silently and without complaint
- Component signs are arbitrary and flip between refits and library versions — pin a convention, such as forcing the largest-magnitude loading positive, or downstream features silently invert
- The mean vector and the scaler are part of the model; recomputing either at inference changes every score
- Never form the covariance matrix to get the components; the squared condition number damages exactly the small components a reconstruction detector depends on
- A refit rotates the coordinate system, so a downstream model trained on old scores must be retrained together with the decomposition, not separately
- Near-equal singular values mean the individual components are not identified; only the subspace they span is, so interpreting one of them is unfounded
Pros and cons, in context
Deterministic, closed-form, with no hyperparameter that changes the answer
Given centred data there is exactly one solution up to sign — no seed, no restarts, no local optimum. That is why it is trusted as a preprocessing step where a stochastic method would need validating in its own right.
Removes correlation and compresses in one operation
Decisive ahead of distance-based and covariance-based methods, which degrade quickly on correlated high-dimensional input. Worth little when the features were already independent, where it just rotates them.
The reconstruction error is a ready-made anomaly score
One fitted object serves compression, denoising, and detection, with per-feature residual contributions that point at the responsible channel. This is why it has survived in process control for decades.
Cheap to compute and essentially free to apply
A matrix product at inference and a small model artefact. Randomized SVD keeps the fit affordable even on very wide data, which is what makes it usable where a nonlinear method would not be.
Maximizes variance, which is not the same as retaining information
The direction that separates two classes can be low-variance and discarded, and nothing in the method warns you. This is the failure that matters most in practice, and it is why supervised alternatives exist for supervised problems.
Linear, so a curved manifold is approximated badly
Data lying on a curve needs far more components than its intrinsic dimension. Kernel PCA, UMAP, and autoencoders exist for exactly this — at the cost of the determinism that makes PCA trustworthy.
Components are rarely interpretable and often not identified
Each is a dense combination of every original feature, signs are arbitrary, and near-equal singular values leave the rotation undetermined. Component interpretations get published anyway, and most of them are not defensible.
Squared error makes it sensitive to outliers
One extreme point can determine the first component by itself. A serious problem in anomaly-adjacent work, where the outliers are the subject — robust PCA exists precisely because accepting this is not viable there.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""PCA by power iteration and deflation - the objective, transcribed.
Maximizing tr(W^T Sigma W) subject to orthonormality means taking the top
eigenvectors of the covariance. Power iteration finds the top one by repeated
multiplication; deflation removes it and repeats, which is the definition of
"the next direction of most remaining variance" written as an algorithm.
"""
import math
def centre(X):
"""Centring is part of the model, not preprocessing.
Without it the first component points at the mean rather than at the
direction of greatest spread, which is a different problem entirely.
"""
n = len(X)
d = len(X[0])
means = [0.0] * d
for j in range(d):
total = 0.0
for i in range(n):
total += X[i][j]
means[j] = total / n
centred = [[X[i][j] - means[j] for j in range(d)] for i in range(n)]
return centred, means
def covariance(X):
"""Sigma = X^T X / (n - 1) for centred X."""
n = len(X)
d = len(X[0])
sigma = [[0.0] * d for _ in range(d)]
for j in range(d):
for k in range(d):
total = 0.0
for i in range(n):
total += X[i][j] * X[i][k]
sigma[j][k] = total / (n - 1)
return sigma
def power_iteration(matrix, iterations=500, tol=1e-10):
"""Top eigenvector: repeated multiplication converges to it.
Each multiplication scales every eigen-direction by its eigenvalue, so the
largest one dominates more with every step. The convergence rate is the
ratio of the top two eigenvalues - close eigenvalues converge slowly, which
is the same degeneracy that makes those components unidentified.
"""
d = len(matrix)
vector = [1.0 / math.sqrt(d)] * d
for _ in range(iterations):
product = [0.0] * d
for j in range(d):
total = 0.0
for k in range(d):
total += matrix[j][k] * vector[k]
product[j] = total
norm = math.sqrt(sum(value * value for value in product))
if norm == 0.0:
break
product = [value / norm for value in product]
shift = sum(abs(product[j] - vector[j]) for j in range(d))
vector = product
if shift < tol:
break
# Rayleigh quotient gives the eigenvalue for the converged vector.
eigenvalue = 0.0
for j in range(d):
for k in range(d):
eigenvalue += vector[j] * matrix[j][k] * vector[k]
return vector, eigenvalue
def fit(X, k):
centred, means = centre(X)
sigma = covariance(centred)
d = len(sigma)
components = []
variances = []
for _ in range(k):
vector, eigenvalue = power_iteration(sigma)
components.append(vector)
variances.append(eigenvalue)
# Deflation: subtract this component's contribution so the next
# iteration finds the largest direction of what remains.
for j in range(d):
for m in range(d):
sigma[j][m] -= eigenvalue * vector[j] * vector[m]
return components, variances, means
def transform(components, means, x):
"""Project one point onto the retained subspace."""
scores = []
for component in components:
total = 0.0
for j in range(len(x)):
total += (x[j] - means[j]) * component[j]
scores.append(total)
return scoresO(n*d^2) to build the covariance plus O(k*d^2) per power iteration, in interpreter loops over lists of lists.
How this language actually wants it written.
"""PCA by SVD - typed, sign-stabilized, factorizing the data not the covariance."""
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
@dataclass(frozen=True)
class PrincipalComponents:
"""A fitted decomposition. The mean is part of the model, not preprocessing."""
components: Matrix # (k, d), rows are the directions
explained_variance: Vector # (k,)
mean: Vector # (d,)
@property
def explained_variance_ratio(self) -> Vector:
return self.explained_variance / self.explained_variance.sum()
def transform(self, X: Matrix) -> Matrix:
if X.ndim != 2 or X.shape[1] != self.mean.size:
raise ValueError(f"expected (n, {self.mean.size}) input, got {X.shape}")
return (X - self.mean) @ self.components.T
def inverse_transform(self, Z: Matrix) -> Matrix:
return Z @ self.components + self.mean
def reconstruction_error(self, X: Matrix) -> Vector:
"""Squared distance from the subspace - the Q statistic, and the
anomaly score. Computed as a residual rather than by reconstructing
and subtracting, which would materialize a second copy of X."""
residual = (X - self.mean) - self.transform(X) @ self.components
return np.einsum("ij,ij->i", residual, residual)
def fit(X: Matrix, n_components: int, standardize: bool = False) -> PrincipalComponents:
"""Fit by SVD of the centred data. Raises ValueError on malformed input."""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
if not 1 <= n_components <= min(X.shape):
raise ValueError(f"n_components must lie in [1, {min(X.shape)}], got {n_components}")
mean = X.mean(axis=0)
centred = X - mean
if standardize:
scale = centred.std(axis=0)
scale[scale == 0.0] = 1.0
centred = centred / scale
# SVD of the DATA, never an eigendecomposition of X^T X. The two give the
# same answer in exact arithmetic; forming the covariance squares the
# condition number and loses roughly half the available precision - and it
# loses it precisely in the small components, which is where a
# reconstruction-based detector does its work.
_, singular_values, right = np.linalg.svd(centred, full_matrices=False)
components = right[:n_components]
explained = (singular_values[:n_components] ** 2) / (X.shape[0] - 1)
# Signs are arbitrary - a component and its negation are equally valid.
# Pinning a convention here is what stops downstream features silently
# inverting between refits and library versions.
for index, component in enumerate(components):
if component[np.abs(component).argmax()] < 0.0:
components[index] = -component
return PrincipalComponents(
components=components,
explained_variance=explained,
mean=mean,
)What changed: The algorithm changes for a reason that is stated in the code: power iteration on an explicitly-formed covariance matrix is replaced by an SVD of the centred data. The two agree in exact arithmetic, but building X-transpose-X squares the condition number and destroys roughly half the precision — and it destroys it in the small components, which is exactly where a reconstruction-error detector operates. Deflation disappears with it, since the SVD returns every component at once. The other addition is a sign convention: components are only defined up to sign, and pinning one here is what stops downstream features silently inverting between refits.
O(n*d*min(n,d)) for the full SVD, executed in LAPACK. Exact, and it computes every component whether or not you want them.
Faster, and the justification for each change.
"""PCA by randomized SVD - k components without factorizing the whole matrix."""
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
def randomized_pca(
X: Matrix,
n_components: int,
oversampling: int = 10,
power_iterations: int = 2,
seed: int = 0,
) -> tuple[Matrix, Vector, Vector]:
"""Top-k components in O(n*d*k) instead of O(n*d*min(n,d)).
A full SVD computes every component and then throws most of them away. The
randomized method instead finds a small orthonormal basis that captures the
dominant action of the matrix, projects onto it, and factorizes the tiny
result - so the cost scales with the components actually wanted rather than
with the smaller matrix dimension.
It is APPROXIMATE. The approximation is excellent when the spectrum decays
quickly and poor when it is flat, which is exactly when the trailing
components are not identified anyway. Power iterations sharpen the decay at
the cost of one pass over the data each.
"""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
rng = np.random.default_rng(seed)
n, d = X.shape
sketch_width = min(n_components + oversampling, d)
# One contiguous float64 block: every pass below is a GEMM over it, and a
# non-contiguous or float32 input would force LAPACK to copy on each one.
design = np.ascontiguousarray(X, dtype=np.float64)
mean = design.mean(axis=0)
np.subtract(design, mean, out=design) # centre in place; no second copy
# Range finder: Y = X Omega, a random projection whose columns span
# approximately the same space as the top singular vectors.
omega = rng.standard_normal((d, sketch_width))
sketch = design @ omega
# Power iterations pull the sketch toward the dominant subspace. Each one
# is two GEMMs, with a re-orthonormalization between them so the small
# singular directions are not lost to round-off.
for _ in range(power_iterations):
sketch, _ = np.linalg.qr(sketch, mode="reduced")
sketch = design.T @ sketch
sketch, _ = np.linalg.qr(sketch, mode="reduced")
sketch = design @ sketch
basis, _ = np.linalg.qr(sketch, mode="reduced") # (n, sketch_width)
# Project into the small basis and factorize there: the SVD is now
# (sketch_width, d) instead of (n, d), which is the whole saving.
projected = basis.T @ design
_, singular_values, right = np.linalg.svd(projected, full_matrices=False)
components = right[:n_components].copy()
explained = (singular_values[:n_components] ** 2) / (n - 1)
for index, component in enumerate(components):
if component[np.abs(component).argmax()] < 0.0:
np.negative(component, out=components[index])
return components, explained, meanWhat changed: A full SVD computes every component and discards most of them, which is pure waste when k is small and d is large. The randomized method replaces that with a sketch: project the data through a random matrix to find a small basis that captures its dominant action, then factorize inside that basis, so the expensive decomposition is on a k-by-d matrix rather than an n-by-d one. Everything is expressed as GEMMs over one contiguous block, and the centring is done in place because on wide data a second copy of the design matrix is the binding memory constraint. The approximation is stated where it will be read: it is excellent on a decaying spectrum and poor on a flat one.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
The range finder, the power iterations, and the projection are all GEMMs over the full matrix, and the only expensive factorization is on a small projected block.
Trade-off: The result is approximate — on a slowly decaying spectrum the trailing components are noticeably wrong, and no amount of oversampling fixes a genuinely flat spectrum.
Pre-allocate output arrays and use in-place operations
Centring is done in place and the sign flip writes back through the component array, so a wide design matrix is never duplicated at the moment memory is tightest.
Trade-off: The caller’s array is copied once at entry and then mutated, so the function owns that copy and the operation order can no longer be rearranged safely.
Ensure contiguous memory layout and a single dtype
One C-contiguous float64 block means every GEMM and QR reads it directly rather than through an internal LAPACK copy, which would otherwise be paid on each of several passes.
Trade-off: Forces a full copy when the caller passes Fortran-ordered or float32 data, doubling peak memory for the duration — significant precisely on the wide matrices this routine exists for.
O(n*d*k) for k components, versus O(n*d*min(n,d)) for the exact SVD. Illustrative, not a measured benchmark.