Principal Component Analysis (PCA) & Autoencoders

Recap & Introduction Kernel PCA Autoencoders Demo Denoising Autoencoders Manifolds

Recap & Introduction

Regression, classification, and SVMs are all supervised methods. They learn from labeled examples. We now turn to unsupervised learning, where the goal is to discover structure in data without labels.

The central technique is dimensionality reduction: finding low-dimensional representations that capture the essential variation in high-dimensional data. This page builds on eigenvalues and eigenvectors and the SVD from our linear-algebra section, and extends PCA to nonlinear settings using both kernels and neural networks.

Principal Component Analysis (PCA) is one of the most fundamental techniques in machine learning and statistics for dimensionality reduction. It provides a method to reduce the number of variables in high-dimensional datasets while retaining the most meaningful structure and variation present in the original data.

PCA begins by analyzing the covariance structure of the data. Given a dataset, we compute the covariance matrix to understand how different features co-vary. Since the covariance matrix is always symmetric and positive semi-definite, it can be orthogonally diagonalized. This is the finite-dimensional shadow of a more general result. That result is the spectral theorem for compact self-adjoint operators on Hilbert spaces, and it underlies the kernel PCA construction below. The resulting eigenvectors, called principal components (PCs), form an orthonormal basis that captures the directions of maximum variance. The corresponding eigenvalues indicate how much of the total variance is captured by each component.

By projecting the data onto the subspace spanned by the top \(k\) principal components (those associated with the largest eigenvalues), PCA identifies a lower-dimensional representation that retains as much variance as possible. This allows us to reduce dimensionality by discarding less informative directions (that is, those with small variance), thereby simplifying the dataset while minimizing information loss.

For example, if a dataset in \(\mathbb{R}^{10}\) has 3 principal components capturing 95% of the total variance, we can project the original data onto a 3D subspace. This projection preserves the dominant patterns and relationships in the data and filters out noise and redundancy. The transformed vectors in the low-dimensional space are known as latent representations.

While PCA is a powerful tool, it identifies directions of maximum variance using linear combinations of the original features. As a result, it may fail to uncover complex, nonlinear structures in the data.

Kernel PCA

Given an \(N \times D\) data matrix, PCA requires the eigenvectors of the \(D \times D\) covariance matrix \(\frac{1}{N} X^\top X\) of the centered data matrix \(X\). The factor \(\frac{1}{N}\) does not change the eigenvectors. If \(D \gg N\), working with the \(N \times N\) Gram matrix \(K = X X^\top\), whose entries are the inner products \(K_{ij} = x_i^\top x_j\), is much more efficient. In kernel PCA below, it becomes \(K = \Phi \Phi^\top\) with entries \(\langle \phi(x_i), \phi(x_j) \rangle\).

Kernel PCA (kPCA) is a nonlinear generalization of classical PCA that uses the kernel trick, which allows us to replace inner products \(x_i^\top x_j\) with a kernel function \(K_{ij} = \mathcal{K}(x_i, x_j)\). The spectral structure that makes this work is formalized by Mercer's theorem (for continuous positive-definite kernels on a compact domain) and, more generally, by the reproducing kernel Hilbert space (RKHS) construction developed in our functional-analysis section.

Kernel PCA implicitly replaces \(x_i\) with \(\phi(x_i) = \phi_i\). Let \(\Phi\) be the corresponding design matrix. Assuming the features are centered, the covariance matrix in feature space is represented by: \[ S_{\phi} = \frac{1}{N} \sum_i \phi_i \phi_i^\top. \]

The normalized eigenvectors of \(S_{\phi}\) with nonzero eigenvalues are given by: \[ V_{kPCA} = \Phi^\top U \Lambda^{-\frac{1}{2}} \] where the columns of \(U\) are orthonormal eigenvectors of the kernel (Gram) matrix \(K = \Phi \Phi^\top\) for its positive eigenvalues, which form the diagonal matrix \(\Lambda\). Indeed, since \(S_{\phi} = \frac{1}{N} \Phi^\top \Phi\), if \(K u = \lambda u\) with \(\lambda \gt 0\), then \[ \begin{align*} S_{\phi} \Phi^\top u &= \frac{1}{N} \Phi^\top K u \\\\ &= \frac{\lambda}{N} \Phi^\top u, \\\\ \|\Phi^\top u\|^2 &= u^\top K u = \lambda, \end{align*} \] so \(\Phi^\top u / \sqrt{\lambda}\) is a unit eigenvector of \(S_{\phi}\) with eigenvalue \(\lambda / N\).

Since \(\phi_i\) can be infinite dimensional, we cannot compute \(V_{kPCA}\) directly. Instead, we express the projection of a test vector \(x_*\) entirely in terms of kernel evaluations:

Proposition: Kernel PCA Projection

\[ \phi_*^\top V_{kPCA} = \phi_*^\top \Phi^\top U \Lambda^{-\frac{1}{2}} = k_*^\top U \Lambda^{-\frac{1}{2}} \] where \(\phi_* = \phi(x_*)\) and \(k_* = \left[\mathcal{K}(x_*, x_1), \ldots, \mathcal{K}(x_*, x_N)\right]^\top\) is the vector of kernel evaluations between the test point and all training points.

Note that using \(K = \Phi \Phi^\top\) is valid only if the features have zero empirical mean, \(\frac{1}{N} \sum_{i=1}^N \phi_i = 0\). However, the feature space can be infinite dimensional, so we cannot subtract off the mean. Here, we introduce the double centering trick.

Let the centered feature vector be \[ \tilde{\phi}_i = \phi(x_i) - \frac{1}{N} \sum_{j =1}^N \phi(x_j). \] Its Gram matrix is given by \[ \tilde{K}_{ij} = \tilde{\phi}_i^\top \tilde{\phi}_j. \]

Proposition: Double Centering Trick

The centered Gram matrix is computed as: \[ \begin{align*} \tilde{K} &= C_N K C_N \\\\ &= K - \frac{1}{N}JK - \frac{1}{N}KJ + \frac{1}{N^2}J K J \end{align*} \] where \[ C_N = I_N - \frac{1}{N}1_N 1_N^\top \] is the centering matrix, \(1_N \in \mathbb{R}^N\) is the all-ones vector, and \(J = 1_N 1_N^\top\) is the \(N \times N\) all-ones matrix.

Proof.

We compute the entries of \(\tilde{K}\) directly: \[ \begin{align*} \tilde{K}_{ij} &= \tilde{\phi}_i^\top \tilde{\phi}_j \\\\ &= \left(\phi_i - \frac{1}{N} \sum_{k=1}^N \phi_k \right)^{\!\top} \left(\phi_j - \frac{1}{N} \sum_{l=1}^N \phi_l \right) \\\\ &= \phi_i^\top \phi_j - \frac{1}{N} \sum_{l=1}^N \phi_i^\top \phi_l - \frac{1}{N} \sum_{k=1}^N \phi_k^\top \phi_j + \frac{1}{N^2} \sum_{k=1}^N \sum_{l=1}^N \phi_k^\top \phi_l. \end{align*} \]

Every inner product on the right is an entry of the uncentered Gram matrix, since \(\phi_i^\top \phi_j = K_{ij}\): \[ \tilde{K}_{ij} = K_{ij} - \frac{1}{N} \sum_{l=1}^N K_{il} - \frac{1}{N} \sum_{k=1}^N K_{kj} + \frac{1}{N^2} \sum_{k=1}^N \sum_{l=1}^N K_{kl}. \] Every entry of \(J = 1_N 1_N^\top\) equals 1, so each of these three sums is itself an entry of a matrix product: \[ (KJ)_{ij} = \sum_{l=1}^N K_{il}, \quad (JK)_{ij} = \sum_{k=1}^N K_{kj}, \quad (JKJ)_{ij} = \sum_{k=1}^N \sum_{l=1}^N K_{kl}. \] Substituting these identities gives \[ \tilde{K} = K - \frac{1}{N}JK - \frac{1}{N}KJ + \frac{1}{N^2}J K J, \] and expanding \(C_N K C_N = \left(I_N - \frac{1}{N}J\right) K \left(I_N - \frac{1}{N}J\right)\) produces exactly these four terms.

The same conclusion follows in one step at the matrix level. With the feature vectors stacked as the rows of \(\Phi\), so that \(K = \Phi \Phi^\top\), the centered design matrix is \(\tilde{\Phi} = C_N \Phi\). Since \(C_N\) is symmetric, \(\tilde{K} = \tilde{\Phi} \tilde{\Phi}^\top = C_N \Phi \Phi^\top C_N^\top = C_N K C_N\).

Autoencoders

Data reconstruction serves as the primary quality control mechanism for dimensionality reduction. When we compress data from the original high-dimensional feature space \(\mathbb{R}^D\) to low-dimensional space \(\mathbb{R}^L\) (where \(L \lt D\)), we need to ensure that the essential structure and information of the original data are preserved. The reconstruction error provides a quantitative measure of information loss. If reconstruction is poor, the learned representation is inadequate for the task at hand.

Reconstruction ensures that the learned latent representation \(z = f_e(x)\) captures the most relevant and meaningful features of the data. The function \(f_e : \mathbb{R}^D \to \mathbb{R}^L\) is called the encoder. If the decoder \(f_d : \mathbb{R}^L \to \mathbb{R}^D\) can successfully reconstruct the original input \(x\) from \(z\), it demonstrates that \(z\) contains sufficient information about the underlying data structure.

The reconstruction function is \(r(x) = f_d \left(f_e(x) \right)\), and the encoder-decoder pair is trained to minimize the negative log-likelihood \(\mathcal{L}(\theta) = -\log p(x | r(x))\). For a Gaussian likelihood with fixed variance, this is the squared reconstruction error up to a positive factor and an additive constant. PCA is the case in which \(f_e\) and \(f_d\) are linear maps.

We can implement \(f_e\) and \(f_d\) by neural networks. This is called an autoencoder. In particular, a linear autoencoder is equivalent to PCA:

Theorem: Linear Autoencoder
  • Input: \(x \in \mathbb{R}^D\)
  • Hidden units: \(z = W_1 x, \quad W_1 \in \mathbb{R}^{L \times D}, \quad L \lt D\)
  • Output: \(\hat{x} = W_2 z = W_2 W_1 x = Wx, \quad W_2 \in \mathbb{R}^{D \times L}\)

The model is trained by minimizing the squared reconstruction error \[ \mathcal{L}(W) = \sum_{n=1}^N \|x_n - Wx_n\|_2^2. \] If the data are centered, that is, \(\sum_{n=1}^N x_n = 0\), then a minimizer \(\hat{W} = \hat{W}_2 \hat{W}_1\) is the orthogonal projection onto the span of \(L\) orthonormal eigenvectors of the empirical covariance matrix belonging to its \(L\) largest eigenvalues. When the \(L\)-th largest eigenvalue is strictly greater than the \((L+1)\)-th, every minimizer reconstructs each data point \(x_n\) as this projection of \(x_n\).

We state this classical result without proof. It rests on the best low-rank approximation property of the truncated SVD, which lies beyond this page. Without centering, the optimal subspace is determined by the second-moment matrix \(\frac{1}{N} \sum_{n} x_n x_n^\top\) instead of the covariance matrix.

Introducing nonlinearities lets the autoencoder represent curved low-dimensional structure that no linear projection can capture. Moreover, deep autoencoders scale to large datasets far better than kernel PCA, whose \(N \times N\) Gram matrix grows quadratically with the number of data points. In practical applications they have therefore become the more popular choice for nonlinear mapping.

Demo: One Dataset, Three Lenses

This demo passes a single dataset through the three constructions of this page: closed-form linear PCA, kernel PCA computed exactly as derived above, and autoencoders. The point colors mark the generating class only. None of the three methods ever sees a label. All of them must discover the structure from coordinates alone.

The Kernel PCA Tab

Pressing Compute builds the Gram matrix \(K\), centers it by the double centering trick \(\tilde{K} = C_N K C_N\), and extracts the top four eigenpairs of \(\tilde{K}\) by subspace iteration. The solver does not declare convergence on a heuristic. It reports the eigenpairs as converged only when the residual \(\|\tilde{K}u - \lambda u\|\) is at most \(10^{-8} \max(1, \lambda_1)\), where \(\lambda_1\) is the leading eigenvalue. If its iteration cap is reached first, the eigenpairs are reported as uncertified. The certificate (residual, iteration count, converged yes/no) is displayed in the metrics panel.

Training scores are then computed in the projection form derived above, \(\tilde{k}^\top U \Lambda^{-1/2}\), and the explained variance ratios divide each \(\lambda_j\) by \(\operatorname{tr} \tilde{K}\), which is \(N\) times the total variance in feature space.

A useful sanity check is built into the controls. With the linear kernel, kernel PCA reproduces ordinary PCA exactly. The two score plots agree up to a sign, and the kernel eigenvalues are \(N\) times the covariance eigenvalues. The kernel trick changes nothing until the kernel does.

The Autoencoder Tab

The linear autoencoder is precisely the architecture of the theorem above: \(z = W_1 x\), \(\hat{x} = W_2 z\), no biases, trained by plain full-batch gradient descent on the mean squared reconstruction error. The data are centered (and rescaled) before training, as the theorem requires. Watching it train is watching the theorem. The learned one-dimensional subspace (solid line) rotates onto the first principal component (dashed line), and the loss curve settles onto the dashed reference level \(\lambda_2\). That level is the variance discarded by a one-component PCA, and it is exactly the best any linear reconstruction can do.

Switching to the nonlinear model (a small \(\tanh\) network with the same one-dimensional bottleneck) breaks that floor. On the circles dataset its loss dives well below \(\lambda_2\), because a curve through the data can follow the rings where a line cannot. The one-dimensional latent panel shows what the bottleneck learned. Stopping a training run midway discards the partial model entirely. Nothing partially trained is ever displayed as a result.

Experiment Suggestions

  1. Circles, RBF kernel, \(\gamma = 3\). The two rings, hopelessly concentric for PCA, already fall into clearly distinct groups in the plane of the first two kernel principal components, although not always on opposite sides of a straight line. Then push \(\gamma\) to 30. The kernel becomes so local that the global ring structure fragments and the separation degrades.
  2. Two moons, RBF kernel, \(\gamma = 10\). A single kernel principal component already separates the classes far better than either linear principal component.
  3. Linear kernel, any dataset. Confirm the PCA equivalence visually. The two score panels should show the same picture up to reflection.
  4. Blobs. Here linear PCA is already adequate. That is a reminder that nonlinear machinery earns its cost only when the structure is nonlinear.
  5. Linear autoencoder on any dataset. Verify that the final loss matches \(\lambda_2\) in the metrics panel to several decimal places, then train the nonlinear model on circles and watch it go below.

Denoising Autoencoders

Denoising autoencoders (DAE) are a more regularized variant of standard autoencoders that add noise to the input during training, then learn to reconstruct the original, uncorrupted data. The theoretical consequence, stated below, is that the optimal denoiser implicitly learns the score function of the data distribution.

The training process involves corrupting the input \(x\) to produce \(\tilde{x}\), typically using Gaussian noise: \[ p_c (\tilde{x} | x) = \mathcal{N}(\tilde{x} | x, \sigma^2 I) \]

The model then minimizes the reconstruction error between its output \(r(\tilde{x})\) and the clean input \(x\): \[ \ell (x, r(\tilde{x})) = \| e \|_2^2 \] where \(e(x) = r(\tilde{x}) - x \) is the residual error for a sample \(x\).

Theorem: Denoising and Score Estimation

As \(\sigma \to 0\), the optimal reconstruction function \(r^*\), the minimizer of the expected reconstruction error above, satisfies: \[ r^*(x) - x = \sigma^2 \nabla_x \log p(x) + o(\sigma^2). \] The denoising autoencoder implicitly learns the score function \(\nabla_x \log p(x)\) (gradient of log-density). This identity is due to Vincent (2011) and Alain & Bengio (2014). The statement is taken on faith here, as the underlying score-matching theory sits beyond the scope of this page.

This score function forms a vector field over the entire feature space. At each point, this vector field indicates the direction and magnitude to move toward regions of higher data density. The reconstruction process follows these vectors, effectively "flowing" corrupted points back to the data manifold along the steepest ascent of the probability landscape.

Data Manifolds in Our Demo

The autoencoder tab in our demo visualizes data manifolds. When the 2D data is compressed through a 1D bottleneck, every reconstruction lies on a single curve, the image of the decoder, so the network must fit the data with a 1-D curve. The Data & reconstruction panel shows each point mapped onto this learned curve. This is the autoencoder discovering and representing the intrinsic lower-dimensional structure of the data.

Note. Lipschitz continuity is a common regularity condition. A function \(f: (X, d_X) \to (Y, d_Y)\) is Lipschitz with constant \(L \geq 0\) if \(d_Y(f(x_1), f(x_2)) \leq L \, d_X(x_1, x_2)\) for all \(x_1, x_2 \in X\). The full development appears in our continuity page in the analysis section. It includes the hierarchy Lipschitz \(\Rightarrow\) uniformly continuous \(\Rightarrow\) continuous.

This means the function's output changes at most linearly with respect to the input. In autoencoders, Lipschitz continuity in the reconstruction map \(r(x)\) ensures stability. Small changes in input lead to small changes in reconstruction.

In the context of denoising autoencoders, enforcing or assuming Lipschitz continuity makes the learned vector field well-behaved. Under a Lipschitz condition, the flow of this vector field exists and is unique (the Picard-Lindelöf theorem, which we do not prove here), and the reconstructions cannot change sharply or unstably. This is essential when approximating the gradient \(\nabla_x \log p(x)\).

Manifolds

Intuitively, a manifold is a topological space that locally resembles Euclidean space near each point. Imagine a curved surface like a sphere. Zooming in on any small patch makes it look flat, like \(\mathbb{R}^2\). In machine learning, we often assume data lies on a low-dimensional manifold embedded in high-dimensional space.

Definition: Manifold (Intuitive Form)

An \(n\)-dimensional manifold \(\mathcal{M}\) is a topological space where every point \(p \in \mathcal{M}\) has an open neighborhood \(U\) that is homeomorphic to \(\mathbb{R}^n\). That is, there exists a continuous bijection \(\phi : U \to \mathbb{R}^n\) with continuous inverse.

This is the intuitive form sufficient for the manifold-hypothesis discussion that follows. The full differential structure (charts, atlases, smoothness) is developed rigorously in our smooth-manifold series.

The Manifold Hypothesis. Real-world high-dimensional data (for example, images and speech) tends to lie on or near a low-dimensional manifold embedded in the ambient space.

For example, each \(64 \times 64\) grayscale face image can be represented as a point in \(\mathbb{R}^{4096}\), but the set of "realistic" face images occupies only a small, structured region of this space. That region is likely a nonlinear manifold of much lower dimension, governed by factors such as pose, lighting, expression, and identity.

Dimensionality reduction reveals low-dimensional structure in data, but does not assign data points to groups. The clustering page addresses this complementary unsupervised task: partitioning data into meaningful groups without labels. We will see how K-means minimizes distortion in the original feature space, while spectral clustering leverages the graph Laplacian to discover clusters that respect the geometry of the data manifold.

Looking further ahead, the manifold hypothesis introduced here is the entry point to a much broader program. Tangent spaces, Riemannian metrics, and the maps between manifolds await in our manifold series, which generalizes the intuitive picture above. On the discrete side, sampled data on graph-structured domains motivate graph neural networks, whose spectral variants use the eigenstructure of the graph Laplacian as a learnable signal-processing primitive. Both threads converge in Geometric Deep Learning (GDL), where neural-network architectures are designed to respect the symmetries and geometric structure of the data they operate on.