⏱️ Lectura: 15 min
Reduce a matrix with a million columns down to ten, and you lose less information than you’d expect: that’s how SVD and PCA work, the linear algebra duo behind recommendation systems, image compression, and much of modern machine learning.
📑 En este artículo
- TL;DR
- What SVD and PCA Are, and Why They Matter
- How the Decomposition Works Under the Hood
- Practical Examples: From a Minimal Matrix to a Real Dataset
- How to Start Using PCA on Your Own Dataset
- Real-World Use Cases
- Common Mistakes and Best Practices
- Comparison: SVD/PCA vs. Other Reduction Techniques
- Going Deeper: What Happens in High Dimensions
- Frequently Asked Questions
- References
A book published on July 11, 2026 on arXiv by Afonso Bandeira, Amit Singer, and Thomas Strohmer, titled “Mathematics of Data Science,” devotes an entire chapter to this technique alongside 15 others covering the mathematical foundations of data science: from high-dimensional geometry to deep neural networks. This article takes that chapter as a starting point and explains it from scratch, with Python code.
TL;DR
- You’ll understand what singular value decomposition (SVD) does mathematically, step by step.
- You’ll be able to implement PCA from scratch with NumPy, without relying on a black box.
- You’ll know how to choose how many principal components to keep using explained variance.
- You’ll be able to tell when to use PCA/SVD and when t-SNE, UMAP, or an autoencoder is a better fit.
- You’ll understand why the curse of dimensionality makes distances lose meaning.
- You’ll be able to compress a matrix or image by computing its low-rank approximation with SVD.
- You’ll know how to verify your PCA is computed correctly by checking explained_variance_ratio_.
What SVD and PCA Are, and Why They Matter
Singular value decomposition (SVD) is a factorization that exists for any matrix, regardless of its shape or whether it’s square. It takes a matrix A of size m x n and rewrites it as the product of three matrices: U, Sigma, and V transpose. A doesn’t need to represent anything special: it can be a customer table, the pixels of an image, or user ratings for movies.
Principal Component Analysis (PCA) is, at its core, a direct application of SVD: it’s used to find the directions where data varies the most and project the information onto those directions, discarding the noise. When a dataset has 500 columns but only 8 explain 95% of the behavior, PCA is the tool that reveals it.
Why it matters: almost all machine learning works with high-dimensional data matrices (images, vectorized text, sensor signals, genes). SVD and PCA are the mathematical foundation that makes it possible to compress that information without destroying it, speed up training, remove noise, and visualize data that would otherwise be impossible to plot on two or three axes. The Eckart-Young theorem, one of the results the Bandeira, Singer, and Strohmer book develops in detail, proves something very strong: keeping the first k singular values gives you the best possible rank-k approximation of your original matrix, measured in Frobenius norm. It’s not a heuristic, it’s a mathematical optimum.
How the Decomposition Works Under the Hood
Every matrix A with m rows and n columns can be written as A = U · Sigma · VT, where U is an orthogonal m x m matrix (its columns are the left singular vectors), V is another orthogonal n x n matrix (columns: right singular vectors), and Sigma is a diagonal matrix with the singular values, always positive and ordered from largest to smallest.
The singular values measure how much “energy” or variance each direction contributes. The first singular value corresponds to the direction where data varies the most; the second, to the next orthogonal direction with the most variance, and so on. If the first k singular values concentrate almost all the energy, you can discard the rest without losing relevant information: that’s the mathematical basis of compression with SVD.
flowchart TD
A["Matrix A (m rows x n columns)"] --> U["U: left singular vectors"]
A --> S["Sigma: ordered singular values"]
A --> V["V transpose: right singular vectors"]
U --> R["A = U x Sigma x V-transpose"]
S --> R
V --> R
PCA reuses this exact decomposition, but first centers the data (subtracts the mean of each column) so the first singular value captures real variance rather than the offset from the origin. The right singular vectors (rows of V transpose) become the principal components: new directions, linear combinations of the original columns, ordered by how much variance they explain. A key property: U and V are orthogonal, which means they rotate the data without distorting it, only changing the reference frame to align the axes with the real variance.
Practical Examples: From a Minimal Matrix to a Real Dataset
Let’s start with the simplest possible case: a 2×2 matrix, to see what numpy.linalg.svd returns without any noise involved.
import numpy as np
A = np.array([
[4, 0],
[3, -5]
])
U, S, Vt = np.linalg.svd(A)
print("U:\n", U)
print("Singular values:", S)
print("Vt:\n", Vt)
This block decomposes a 2×2 matrix into its three pieces. S returns a vector with the singular values (not the full diagonal matrix, NumPy simplifies it), ordered from largest to smallest. You can reconstruct the original matrix with U @ np.diag(S) @ Vt and confirm it gives the same result as A: that exact reconstruction is what later gets truncated to the first k terms to compress.
Now a more realistic case: a simulated dataset of 200 customers with 50 numeric variables (purchase behavior, for example), reduced to 2 dimensions with PCA implemented manually on top of SVD:
import numpy as np
def pca_manual(X, k):
X_centrado = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(X_centrado, full_matrices=False)
componentes = Vt[:k]
X_reducido = X_centrado @ componentes.T
varianza_explicada = (S[:k] ** 2) / np.sum(S ** 2)
return X_reducido, varianza_explicada
rng = np.random.default_rng(42)
datos_clientes = rng.normal(size=(200, 50))
datos_reducidos, varianza = pca_manual(datos_clientes, k=2)
print("Original shape:", datos_clientes.shape)
print("Reduced shape:", datos_reducidos.shape)
print("Variance explained per component:", varianza)
The function centers the data, computes the SVD, and projects onto the first k components. The result, datos_reducidos, goes from 50 columns to 2, something you can already plot on a scatter plot. varianza_explicada tells you what percentage of the original information each component retained: if you add up the two values and get something low, as tends to happen with purely random data like in this simulated example, that’s expected, because we didn’t put any real structure into the data, so any direction explains roughly the same amount.
💡 Tip: if your matrix is sparse, like a term-document matrix for text, usescipy.sparse.linalg.svdsorsklearn.decomposition.TruncatedSVDinstead ofnumpy.linalg.svd: the latter densifies the matrix in memory and can crash your process from running out of RAM.
A third example, closer to a real use case: compressing a matrix (for example, the pixels of a grayscale image) by keeping only the first k singular values and measuring the error of that approximation.
import numpy as np
def comprimir_svd(matriz, k):
U, S, Vt = np.linalg.svd(matriz, full_matrices=False)
aproximacion = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]
error_relativo = np.linalg.norm(matriz - aproximacion) / np.linalg.norm(matriz)
return aproximacion, error_relativo
imagen = np.random.default_rng(0).normal(size=(300, 200))
_, error_k10 = comprimir_svd(imagen, k=10)
_, error_k50 = comprimir_svd(imagen, k=50)
print(f"Relative error with k=10: {error_k10:.4f}")
print(f"Relative error with k=50: {error_k50:.4f}")
With real images (photos, scans, X-rays) most of the energy tends to concentrate in the first singular values, because there’s repeated structure (edges, textures, flat areas). In a matrix of pure noise like the one in this example, on the other hand, the error drops much more slowly as you increase k, precisely because there’s no structure to compress.
How to Start Using PCA on Your Own Dataset
To reproduce these examples you need Python 3.10 or higher and two packages: NumPy and, for the production version, scikit-learn.
Windows (PowerShell):
python -m pip install numpy scikit-learn
macOS:
python3 -m pip install numpy scikit-learn
Linux:
pip install numpy scikit-learn
With scikit-learn, the same manual PCA from above shrinks down to these lines, with the advantage that the class already includes validations, different solvers, and additional options:
from sklearn.decomposition import PCA
import numpy as np
rng = np.random.default_rng(42)
datos_clientes = rng.normal(size=(200, 50))
pca = PCA(n_components=2)
datos_reducidos = pca.fit_transform(datos_clientes)
print("Explained variance:", pca.explained_variance_ratio_)
print("Cumulative sum:", pca.explained_variance_ratio_.sum())
The exact parameter to configure is n_components: you can pass it an integer (a fixed number of components) or a float between 0 and 1, for example n_components=0.95, in which case scikit-learn automatically chooses how many components it needs to retain that percentage of variance.
To verify the PCA is computed correctly, check that pca.explained_variance_ratio_.sum() gets close to 1.0 when you use all possible components, and compare it against np.linalg.matrix_rank(datos_clientes) to confirm how many dimensions carry real information (the rank of the matrix) before reducing.
Real-World Use Cases
Collaborative filtering recommendation systems, the kind popularized by the Netflix Prize, factorize the user-item matrix with a variant of SVD to predict what rating a user would give a movie they’ve never seen. The original matrix has millions of empty cells (users who haven’t watched almost anything), and SVD finds latent factors that generalize across those gaps.
Image compression uses SVD directly: each grayscale image is a matrix of pixels, and keeping only the first k singular values reconstructs a smaller version with controlled loss. The more structure the original image has (well-defined edges, repeated areas), the smaller the k you need for a visually acceptable reconstruction.
In bioinformatics, PCA reduces gene expression datasets with tens of thousands of columns (one per gene) down to a handful of components that group patients by biological similarity. In finance, PCA on interest rate series identifies the “factors” (level, slope, curvature) that explain almost all the movement of the entire curve, something risk models use to simulate scenarios with far fewer variables than the individual rates.
Common Mistakes and Best Practices
The most common mistake is applying PCA without standardizing variables that are on different scales: if one column measures age (0-100) and another measures income (0-1,000,000), the income variance dominates the first component without that reflecting any real relationship. The fix is to scale each column to mean 0 and standard deviation 1 before running PCA, using sklearn.preprocessing.StandardScaler.
Another mistake is interpreting principal components as if they had direct business meaning: a component is a linear combination of the original variables, not a brand-new named variable. You need to inspect the weights (pca.components_) to understand which original variables weigh the most in each component before giving it a name.
It’s also common to pick k arbitrarily. The recommended practice is to plot the cumulative explained variance against the number of components (the “scree plot”) and cut where the curve flattens out, or set a threshold like 90-95% of retained variance instead of guessing a round number.
⚠️ Watch out: PCA only captures linear relationships. If your data has curved structure or forms nonlinear clusters, you’ll need t-SNE, UMAP, or an autoencoder: PCA will flatten them into a projection that loses that structure.
Comparison: SVD/PCA vs. Other Reduction Techniques
| Method | When to Use It | Advantage | Limitation |
|---|---|---|---|
| SVD / PCA | Data with linear relationships, when you need speed and interpretability | Deterministic, fast, mathematically optimal | Doesn’t capture nonlinear structure |
| t-SNE | Exploratory 2D/3D visualization of nonlinear clusters | Preserves local neighborhoods | Slow on large datasets, non-deterministic, doesn’t project new data |
| UMAP | Visualization and reduction at larger scale than t-SNE | Faster than t-SNE, better preserves global structure | Sensitive hyperparameters, results vary between runs |
| Autoencoder | Data with highly nonlinear relationships and enough volume to train a network | Can learn any reduction function | Needs more data, GPU, and fine-tuning than PCA |
Going Deeper: What Happens in High Dimensions
When the number of columns grows much faster than the number of rows, for example 20,000 genes measured in just 200 patients, you run into what the Bandeira, Singer, and Strohmer book discusses under the heading “Curses, Blessings, and Surprises in High Dimensions”: the geometric intuition we have in 2 or 3 dimensions stops working.
One concrete example: in high dimensions, almost all the volume of a sphere concentrates near its surface, not at the center. And Euclidean distances between random points tend to become nearly identical to each other, which breaks algorithms that rely on “nearest neighbor,” like k-NN, if the dimensionality isn’t reduced before applying them.
flowchart LR
D1["Few dimensions: varied distances"] --> D2["More dimensions: distances become similar"]
D2 --> D3["Many dimensions: almost all volume is at the surface"]
D3 --> D4["Without reduction: k-NN and clustering lose meaning"]
SVD and PCA don’t solve the curse of dimensionality by magic: what they do is identify how many dimensions actually contribute information, the effective rank of the matrix, and discard the rest, which is usually noise. That’s different from blindly reducing dimensions, and it’s why preprocessing with PCA tends to improve the performance of distance-based algorithms.
The complete flow of a PCA pipeline in production, from raw data to components ready for a downstream model, looks like this:
flowchart TD
A["Raw data"] --> B["Standardize columns"]
B --> C["Center data (subtract the mean)"]
C --> D["Compute SVD of the centered matrix"]
D --> E["Choose k components by cumulative variance"]
E --> F["Project data onto k components"]
F --> G["Reduced data ready for the model"]
Another advanced detail: when the matrix is too large to compute the full SVD, for example with millions of rows, randomized SVD is used, an algorithm that approximates the first k singular values by projecting the matrix onto a random reduced-dimension subspace before decomposing it. sklearn.decomposition.PCA activates it automatically with svd_solver="randomized" when it detects large datasets, and the computational cost shifts from depending on the matrix’s full dimensions to depending only on k, the number of components you requested.
💭 Key takeaway: the effective rank of a real matrix is almost never its full mathematical rank. Most real-world datasets have far less “information” than their number of columns suggests, and that’s why SVD works so well as a compression technique.
📖 Summary on Telegram: View summary
Your next step: take your own dataset in a pandas DataFrame, run PCA(n_components=0.95).fit(df_estandarizado), and plot the cumulative explained variance to see how many real dimensions your problem has.
Frequently Asked Questions
What’s the difference between SVD and PCA?
SVD is the general mathematical operation that factorizes any matrix into three pieces. PCA is a specific application of SVD on centered data, designed to find directions of maximum variance. In practice, PCA is computed using SVD under the hood.
How many principal components should I choose?
There’s no universal number. The standard practice is to set a cumulative explained variance threshold, 90-95% is common, and use pca.explained_variance_ratio_.cumsum() to find how many components reach it.
Does PCA work with categorical data?
Not directly. PCA assumes continuous numeric variables. For categorical data there are variants like MCA (Multiple Correspondence Analysis), or you need to encode the categories numerically first, being careful not to introduce false ordinal relationships.
Is SVD the same as eigendecomposition?
They’re related but not identical. Eigendecomposition only exists for square, diagonalizable matrices. SVD exists for any matrix, and the singular values of A are the square root of the eigenvalues of A transpose times A.
Which Python library should I use in production?
For datasets that fit in memory, numpy.linalg.svd or sklearn.decomposition.PCA are enough. For sparse matrices from text or graphs, scipy.sparse.linalg.svds or sklearn.decomposition.TruncatedSVD avoid densifying the matrix.
Does SVD work for sparse data, like a text matrix?
Yes, with the truncated variant. Computing the full SVD of a giant sparse matrix turns it dense in memory, which can consume unnecessary gigabytes: that’s why TruncatedSVD only computes the first k singular values without densifying the full matrix.
References
- Bandeira, Singer, and Strohmer, “Mathematics of Data Science” (arXiv:2607.11938): the complete book with all 16 chapters on the mathematical foundations of data science.
- Wikipedia: Singular Value Decomposition: formal definition, properties, and proofs of the factorization.
- Official scikit-learn documentation on PCA: class parameters, available solvers, and examples.
- numpy.linalg.svd documentation: reference for the function used in this article’s examples.
📱 Do you like this content? Follow @programacion on Telegram for daily tech content in Spanish: quick summaries, fresh content every day.
Imagen destacada: Foto de Logan Voss en Unsplash
0 Comments