Dimensionality reduction

The mixture models of the previous note explained each point by a discrete latent variable, the cluster it came from. Much data is better explained by a few continuous ones. An $8 \times 8$ image of a written digit has 64 pixel values, but the images vary along far fewer directions: which digit, its slant, its width, the thickness of the stroke. The points fill a thin region of $\mathbb R^{64}$, close to a low-dimensional subspace or surface. Finding that region gives coordinates that compress the data, remove noise, make distances more meaningful and let us look at the data in two dimensions. This note develops the linear answer, principal component analysis, from two directions that meet in one eigenvalue problem, gives it a probabilistic reading as a latent-variable model, and then turns to nonlinear methods that keep neighbourhoods rather than variance.

Notation. $x_1, \dots, x_n \in \mathbb R^d$ are the data points, centred so that $\sum_i x_i = 0$, and $X \in \mathbb R^{n \times d}$ is the matrix with rows $x_i^\top$. $S = \frac1n X^\top X$ is the sample covariance matrix, with eigenvalues $\lambda_1 \ge \dots \ge \lambda_d \ge 0$ and orthonormal eigenvectors $v_1, \dots, v_d$. $k < d$ is the target dimension, $W \in \mathbb R^{d \times k}$ a matrix with orthonormal columns ($W^\top W = I_k$), and $z_i = W^\top x_i \in \mathbb R^k$ the new coordinates of $x_i$. $X = U \Sigma V^\top$ is the singular value decomposition, with singular values $\sigma_1 \ge \sigma_2 \ge \dots$.

Why reduce the dimension#

Four reasons, often at once:

• Looking. Plots have two axes. A map to two dimensions is the only way to see all 64 features of a dataset in one scatter plot.
• Compression. Storing $k$ coordinates instead of $d$ saves memory and makes everything downstream faster.
• Denoising. If the signal lives in a few directions and the noise is spread over all of them, discarding the weak directions removes mostly noise.
• Better distances. In high dimension, distances between points concentrate around a common value (the shell of the Gaussian), which hurts nearest neighbours and the clustering methods of this section. Distances in a well-chosen low-dimensional space are more informative.

PCA as variance maximization#

Project the data onto a single direction, a unit vector $w$. The projections $w^\top x_i$ have mean zero, because the data are centred, and variance

$$\frac1n \sum_{i=1}^n (w^\top x_i)^2 = w^\top \Big(\frac1n \sum_{i=1}^n x_i x_i^\top\Big) w = w^\top S w .$$

A direction that keeps the points spread out keeps the differences between them, so ask for the direction of largest variance: maximize $w^\top S w$ subject to $w^\top w = 1$. The Lagrangian $w^\top S w - \lambda (w^\top w - 1)$ has gradient $2Sw - 2\lambda w$, which vanishes when

$$S w = \lambda w .$$

So $w$ must be an eigenvector of $S$, and then the variance is $w^\top S w = \lambda w^\top w = \lambda$: the best direction is the eigenvector with the largest eigenvalue, $v_1$, the first principal component. The second component is the direction of largest variance orthogonal to the first. Writing any unit $w$ in the eigenbasis, $w = \sum_j c_j v_j$ with $\sum_j c_j^2 = 1$, gives $w^\top S w = \sum_j \lambda_j c_j^2$; with $c_1 = 0$ forced by orthogonality, the maximum is $\lambda_2$, at $w = v_2$. Continuing, the $k$ leading eigenvectors are the $k$ principal components, and the variance of the data along $v_j$ is exactly $\lambda_j$.

The multivariate Gaussian note met the same eigenvectors as the axes of the covariance ellipse. PCA is that geometry estimated from data: the principal components are the axes of the ellipse that best describes the cloud, ordered by length.

PCA as reconstruction-error minimization#

Now ask a different question. Represent each point by $k$ numbers $z_i = W^\top x_i$ and reconstruct it as $\hat x_i = W z_i = W W^\top x_i$. Which $W$ makes the reconstructions closest to the data? Because $W W^\top$ is an orthogonal projection, the residual $x_i - \hat x_i$ is orthogonal to $\hat x_i$, and by Pythagoras

$$\lVert x_i - W W^\top x_i \rVert^2 = \lVert x_i \rVert^2 - \lVert W^\top x_i \rVert^2 .$$

Averaging over the points,

$$\frac1n \sum_{i=1}^n \lVert x_i - \hat x_i \rVert^2 = \operatorname{tr} S - \operatorname{tr}(W^\top S W) .$$

The first term is the total variance, fixed by the data; the second is the variance kept by the projection. Minimizing the reconstruction error is therefore the same problem as maximizing the kept variance, and its solution is the same: $W = (v_1, \dots, v_k)$, with minimal error

$$\frac1n \sum_{i=1}^n \lVert x_i - \hat x_i \rVert^2 = \sum_{j > k} \lambda_j ,$$

the sum of the discarded eigenvalues. (That the $k$ leading eigenvectors maximize $\operatorname{tr}(W^\top S W)$ over all orthonormal $W$ follows from the same eigenbasis argument applied to each column, or from Ky Fan’s theorem.) The two derivations are the two historical origins of PCA: Pearson (1901) looked for the line and plane of closest fit, Hotelling (1933) for components of maximal variance.

Figure

Two sides of one problem: variance and reconstruction error

Worked example. For the covariance $\Sigma = \begin{pmatrix} 2 & 1.2 \\ 1.2 & 1 \end{pmatrix}$ of the probability course, the eigenvalues are $2.8$ and $0.2$ and the first eigenvector is $v_1 = (0.832, 0.555)$. The first component keeps $2.8 / 3 = 93.3\%$ of the variance, and reconstructing from it alone leaves a mean squared error of $0.2$. The point $x = (1, 1)$ gets the coordinate $z = v_1^\top x = 1.387$ and the reconstruction $\hat x = z v_1 = (1.154, 0.769)$; the residual $(-0.154, 0.231)$ has squared length $0.077$, and indeed $\lVert x \rVert^2 = 2 = 1.387^2 + 0.077$.

Through the SVD#

Every matrix has a singular value decomposition $X = U \Sigma V^\top$, with $U \in \mathbb R^{n \times r}$ and $V \in \mathbb R^{d \times r}$ having orthonormal columns and $\Sigma$ diagonal with $\sigma_1 \ge \dots \ge \sigma_r > 0$ ($r$ the rank). Then

$$S = \frac1n X^\top X = \frac1n V \Sigma U^\top U \Sigma V^\top = V \frac{\Sigma^2}{n} V^\top ,$$

an eigendecomposition: the principal components are the right singular vectors $V$, the eigenvalues are $\lambda_j = \sigma_j^2 / n$, and the coordinates of all points are $XV = U\Sigma$. Two reasons to compute PCA this way. Forming $X^\top X$ squares the condition number, so small eigenvalues lose precision, while the SVD works on $X$ directly. And for large data a truncated or randomized SVD computes only the $k$ leading components, without the full $d \times d$ matrix.

The SVD also says that PCA is optimal in a stronger sense. By the Eckart–Young theorem (Eckart and Young, 1936), the best approximation of $X$ by any matrix of rank $k$ in the Frobenius norm is the truncated SVD $X_k = U_k \Sigma_k V_k^\top$, with error $\lVert X - X_k \rVert_F^2 = \sum_{j > k} \sigma_j^2$. The PCA reconstruction $X W W^\top$ with $W = V_k$ is exactly $X_k$: no rank-$k$ linear summary of the data does better. Truncated SVD of a word–document count matrix is latent semantic analysis, the baseline of the GloVe note.

Explained variance and choosing k#

The share of variance kept by the first $k$ components is the explained variance ratio

$$\frac{\lambda_1 + \dots + \lambda_k}{\lambda_1 + \dots + \lambda_d} ,$$

and one minus it is the relative reconstruction error. On the 1797 digit images of scikit-learn, 64 pixels each, the first two components explain 14.9% and 13.6% of the variance; 5 components explain half, 21 explain 90% and 29 explain 95%. With 21 components the mean squared reconstruction error per image is 116.3, the sum of the 43 discarded eigenvalues, as the derivation above predicts.

There is no single rule for $k$. Common choices:

• A share of variance, such as 90% or 95%. Simple, but the threshold is arbitrary.
• The elbow of the scree plot, the eigenvalues in decreasing order: keep the components before the curve flattens into a floor of small, similar eigenvalues, which is what isotropic noise produces.
• The downstream task: choose $k$ by the validation error of the model that uses the components.
• The likelihood of probabilistic PCA (below), on held-out data or with a penalty like BIC.

Figure

How many components to keep

The figure shows the first two rules disagreeing. The data have three directions of signal and noise in all ten; the scree plot shows the three, but they explain only 89.9% of the variance, so a 90% rule keeps a fourth component that is pure noise.

Standardizing and whitening#

PCA depends on the units of the features. A feature measured in grams has a thousand times the standard deviation of the same feature in kilograms and dominates the first component. When features have unrelated units, standardize each to unit variance first; PCA then works with the correlation matrix instead of the covariance. When they share a unit, as pixels do, the raw covariance is usually right.

Whitening goes one step further and rescales each component to unit variance:

$$z_i = \Lambda_k^{-1/2} W^\top x_i, \qquad \Lambda_k = \operatorname{diag}(\lambda_1, \dots, \lambda_k) .$$

The new coordinates are uncorrelated with unit variance, $\frac1n \sum_i z_i z_i^\top = I_k$, the sample version of the whitening map of the multivariate Gaussian note. Whitening helps methods that assume comparable, uncorrelated inputs, such as independent component analysis or distance-based methods. Its risk is the division by $\sqrt{\lambda_j}$: directions with tiny eigenvalues, often pure noise, are amplified to the same scale as the signal, so whitening is combined with truncation to $k$ components or with a small constant added to each $\lambda_j$.

Probabilistic PCA#

PCA as described is a projection, not a model: it assigns no probability to a point. Probabilistic PCA (Tipping and Bishop, 1999) is the continuous analogue of the Gaussian mixture. A latent vector $z \sim \mathcal N(0, I_k)$ is mapped linearly into the data space and blurred by isotropic noise, $x = Wz + \mu + \varepsilon$ with $\varepsilon \sim \mathcal N(0, \sigma^2 I_d)$; by the affine-map rule, $x \sim \mathcal N(\mu, W W^\top + \sigma^2 I_d)$, a Gaussian whose covariance is “rank $k$ plus noise”. Its maximum likelihood solution is PCA in disguise: the columns of $W$ span the $k$ leading eigenvectors, $W = V_k (\Lambda_k - \sigma^2 I)^{1/2} R$ for any rotation $R$, and the noise variance is the average of the discarded eigenvalues, $\sigma^2 = \frac{1}{d - k}\sum_{j > k} \lambda_j$ (Bishop, 2006, Section 12.2.1). In the scree figure that average is $0.250$, the true noise variance; for the digits with $k = 10$ it is $5.83$. The posterior mean of the latent vector is $\mathbb E[z \mid x] = (W^\top W + \sigma^2 I_k)^{-1} W^\top (x - \mu)$. It is not the vector of principal coordinates $V_k^\top (x - \mu)$, because $W$ holds scaled and rotated loadings, not the orthonormal eigenvectors: as $\sigma^2 \to 0$ it tends to $R^\top \Lambda_k^{-1/2} V_k^\top (x - \mu)$, the principal coordinates divided by their standard deviations and rotated by $R$. What becomes the PCA projection is the reconstruction in data space, $W \, \mathbb E[z \mid x] + \mu \to V_k V_k^\top (x - \mu) + \mu$. The model can also be fitted by EM, which avoids the $d \times d$ covariance and handles missing values, and because it is a density it gives a likelihood for choosing $k$ and a score for unusual points, the starting point of the next note.

Kernel PCA#

PCA finds linear structure only. Two concentric circles have no direction of large variance that separates them: every projection to a line mixes the rings. Kernel PCA (Schölkopf, Smola and Müller, 1998) runs PCA in a feature space $\varphi(x)$ reached through a kernel $k(x, x') = \varphi(x)^\top \varphi(x')$, as in the kernels note. The principal directions in feature space are combinations of the mapped points, $\sum_i \alpha_i \varphi(x_i)$, and the eigenvalue problem becomes one for the $n \times n$ kernel matrix $K_{ij} = k(x_i, x_j)$, centred in feature space,

$$\tilde K = K - \mathbf 1_n K - K \mathbf 1_n + \mathbf 1_n K \mathbf 1_n , \qquad \tilde K \alpha = n \lambda \alpha ,$$

where $\mathbf 1_n$ is the $n \times n$ matrix with all entries $1/n$ and $\alpha$ is scaled so that the direction has unit length, $n\lambda\, \alpha^\top \alpha = 1$. The coordinate of a point on a component is $\sum_i \alpha_i \tilde k(x_i, x)$, computed with kernels only. With $k(x, x') = (x^\top x')^2$ the feature space contains $x_1^2$, $x_2^2$ and $x_1 x_2$, so the squared radius $x_1^2 + x_2^2$, which separates the circles, is a linear function there. In the Code section, a logistic regression on two ordinary principal components of the circles reaches a training accuracy of 0.502, chance, and on three kernel components 1.0. Both are fitted and scored on the same 400 points, so they show what each representation can express linearly, not a test score. Only the third kernel component separates the rings by itself: kernel PCA ranks directions by variance in feature space, not by usefulness. The price is an $n \times n$ eigenproblem, and the components have no simple interpretation in the input space.

t-SNE and UMAP: neighbourhoods instead of variance#

For visualization, preserving large distances is often the wrong goal. A 2-D linear projection of the digits keeps the directions of largest spread, and the ten digit classes overlap heavily. t-SNE (van der Maaten and Hinton, 2008) asks instead that points close in the data stay close in the plane.

1. Turn distances in the data into neighbour probabilities. For each point $x_i$, $p_{j \mid i} \propto \exp\big(-\lVert x_i - x_j \rVert^2 / 2\sigma_i^2\big)$, normalized over $j \ne i$, where each bandwidth $\sigma_i$ is chosen so that the distribution has a set perplexity, $2^{H(p_{\cdot \mid i})}$, roughly the number of neighbours that count; typical values are 5 to 50. Symmetrize: $p_{ij} = (p_{j \mid i} + p_{i \mid j}) / 2n$.
2. Define the same kind of probabilities for points $y_i$ in the plane, with a heavy-tailed Student-$t$ kernel: $q_{ij} \propto (1 + \lVert y_i - y_j \rVert^2)^{-1}$, normalized over all pairs.
3. Move the $y_i$ by gradient descent to minimize $D_{\mathrm{KL}}(P \,\|\, Q) = \sum_{i \ne j} p_{ij} \log (p_{ij} / q_{ij})$.

The asymmetry of the KL divergence decides what is preserved. A pair with large $p_{ij}$, true neighbours, placed far apart (small $q_{ij}$) costs a lot; a pair with small $p_{ij}$ placed close costs little. t-SNE keeps neighbourhoods and is free with everything else. The heavy tail answers the crowding problem: a point in 64 dimensions can have many equidistant neighbours, a point in the plane cannot, and the $t$ kernel lets moderately distant points be pushed far away instead of piling up in the middle. In the Code section, 58.4% of the 10 nearest neighbours of each digit image in 64 dimensions are among its 10 nearest neighbours in a t-SNE map, against 11.8% in the 2-D PCA projection.

UMAP (McInnes, Healy and Melville, 2018) follows the same plan with different choices: a weighted $k$-nearest-neighbour graph in the data, a similar graph in the plane, and a cross-entropy between their edge weights optimized by stochastic gradient descent. It is much faster, scales to millions of points and can embed new points, and it is often said to keep more global structure; the reading rules below apply to it as well.

How to read such maps (Wattenberg, Viégas and Johnson, 2016):

• Hyperparameters really matter. The perplexity (for UMAP, the number of neighbours) changes the picture; look at several values.
• Cluster sizes mean nothing. t-SNE expands dense clusters and contracts sparse ones, so the area of a cluster says nothing about its spread in the data.
• Distances between clusters may mean nothing. Two clusters far apart in the map may be close in the data, and the other way round.
• Random noise can look clustered, especially at low perplexity.
• Shapes and topology need more than one plot. A ring or a line in the map may be an artefact of one run.

What a map shows reliably is which points are neighbours. It is for looking, not for measuring: its axes have no meaning, there is no map of new points in plain t-SNE, and its coordinates are poor features for a model. The autoencoder note reads t-SNE plots of learned codes with exactly these cautions, and the word embeddings note compares PCA and t-SNE projections of word vectors.

Autoencoders: nonlinear PCA#

PCA is the best linear encoder–decoder pair: $z = W^\top x$, $\hat x = Wz$, trained to minimize squared reconstruction error. Training a linear network with a $k$-unit bottleneck on the same loss by gradient descent reaches the same subspace: its minimum spans the principal subspace, and it has no other local minima, only saddle points (Baldi and Hornik, 1989). Replacing the two linear maps by neural networks, an encoder $z = f_\varphi(x)$ and a decoder $\hat x = g_\theta(z)$, gives an autoencoder, which can follow a curved surface instead of a flat subspace. The autoencoder note of the generative models course builds convolutional autoencoders for images on this idea, and the probabilistic version, with the encoder as an approximate posterior, is the variational autoencoder, probabilistic PCA made nonlinear.

What changed#

	PCA	Probabilistic PCA	Kernel PCA	t-SNE / UMAP	Autoencoder
Map	linear, $z = W^\top x$	linear, with a density	linear in feature space	none: positions optimized	neural encoder and decoder
Preserves	variance, large distances	the same, plus a likelihood	variance in feature space	neighbourhoods	what reconstruction needs
New points	project	project	kernels with the training set	not in plain t-SNE	encode
Solution	SVD, exact	closed form or EM	$n \times n$ eigenproblem	non-convex, random	non-convex, gradient descent
Use	compression, denoising, features	the same, density, missing data	nonlinear features	visualization only	representations, generation

Limitations#

• Linear and variance-driven. PCA keeps the directions of largest variance, which need not be the directions that matter for a task; a small-variance direction can carry all the class information.
• Sensitive to scale and outliers. Units decide the components, and a single far point can capture one; robust PCA variants exist for that.
• Components are hard to interpret. Each is a mixture of all features, with signs fixed only up to $\pm 1$.
• Nonlinear maps are for looking. t-SNE and UMAP distort sizes and distances, depend on hyperparameters and random seeds, and plain t-SNE gives no function for new points. UMAP has one (transform in the umap-learn library), but the new positions inherit the same distortions; a pipeline that feeds an embedding to a model must fit it on the training split only.
• The intrinsic dimension is a modelling choice. Every rule for $k$ encodes an assumption about what is signal and what is noise.

What you should remember#

• PCA projects centred data onto the leading eigenvectors of the covariance $S$; the variance along $v_j$ is $\lambda_j$.
• Maximizing kept variance and minimizing squared reconstruction error are the same problem, since $\frac1n\sum \lVert x_i - \hat x_i\rVert^2 = \operatorname{tr} S - \operatorname{tr}(W^\top S W)$; the minimal error is $\sum_{j > k} \lambda_j$.
• Compute PCA by the SVD $X = U\Sigma V^\top$: components $V$, eigenvalues $\sigma_j^2/n$, coordinates $U\Sigma$; by Eckart–Young no rank-$k$ approximation is better.
• Choose $k$ by explained variance, the scree elbow, the downstream task or the likelihood; the rules can disagree.
• Standardize features with different units; whitening rescales components to unit variance and amplifies noise directions.
• Probabilistic PCA, $x \sim \mathcal N(\mu, WW^\top + \sigma^2 I)$, has PCA as its maximum likelihood solution and $\sigma^2$ the mean of the discarded eigenvalues.
• Kernel PCA is PCA of the centred kernel matrix: nonlinear components at the cost of an $n \times n$ eigenproblem.
• t-SNE and UMAP preserve neighbourhoods; cluster sizes and gaps between clusters in their maps carry little information.
• A linear autoencoder learns the PCA subspace; a nonlinear one generalizes it.

Code#

PCA by the SVD, explained variance and reconstruction error on the digits, probabilistic PCA, whitening, kernel PCA, and neighbourhoods in PCA and t-SNE maps:

import numpy as np
from sklearn.datasets import load_digits, make_circles
from sklearn.decomposition import PCA, KernelPCA
from sklearn.manifold import TSNE
from sklearn.neighbors import NearestNeighbors
from sklearn.linear_model import LogisticRegression

# the population: eigenvalues of the covariance and the share of the first component
S = np.array([[2.0, 1.2], [1.2, 1.0]])
lam = np.linalg.eigvalsh(S)[::-1]
print(lam.round(3), round(lam[0] / lam.sum(), 3))                  # [2.8 0.2] 0.933

# PCA of a sample through the SVD of the centred data matrix: the same as the eigenvectors of the covariance
rng = np.random.default_rng(0)
X = rng.multivariate_normal([0, 0], S, size=2000)
Xc = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
ev = np.linalg.eigvalsh(Xc.T @ Xc / len(X))[::-1]
print((s ** 2 / len(X)).round(3), ev.round(3), np.abs(Vt[0]).round(3))   # [2.812 0.198] [2.812 0.198] [0.826 0.563]

# digits: 1797 images of 8 x 8 pixels; share of the first two components, components needed for 50%, 90%, 95%
D = load_digits().data
p = PCA().fit(D)
cum = np.cumsum(p.explained_variance_ratio_)
print(p.explained_variance_ratio_[:2].round(3), [int(np.searchsorted(cum, q) + 1) for q in (0.5, 0.9, 0.95)])
# [0.149 0.136] [5, 21, 29]

# reconstruction error with k components = the sum of the discarded eigenvalues (per image, with the 1/n covariance)
for k in (2, 10, 21):
    pk = PCA(k).fit(D)
    err = ((D - pk.inverse_transform(pk.transform(D))) ** 2).sum(axis=1).mean()
    lam_n = p.explained_variance_ * (len(D) - 1) / len(D)
    print(k, round(err, 1), round(lam_n[k:].sum(), 1), round(cum[k - 1], 3))
# 2 858.9 858.9 0.285 | 10 314.5 314.5 0.738 | 21 116.3 116.3 0.903

# probabilistic PCA: the noise variance is the mean of the discarded eigenvalues
p10 = PCA(10).fit(D)
print(round(p10.noise_variance_, 3), round(p.explained_variance_[10:].mean(), 3))   # 5.828 5.828
# whitening: the rescaled coordinates have identity covariance
Z = PCA(10, whiten=True).fit_transform(D)
print(np.allclose(np.cov(Z.T), np.eye(10)))                         # True

# kernel PCA on two concentric circles: training accuracy of a linear classifier on the components
Xr, yr = make_circles(n_samples=400, factor=0.4, noise=0.05, random_state=0)
acc = lambda z: round(LogisticRegression().fit(z, yr).score(z, yr), 3)
z_lin = PCA(2).fit_transform(Xr)
z_ker = KernelPCA(3, kernel='poly', degree=2, gamma=1, coef0=0).fit_transform(Xr)   # k(x, y) = (x . y)^2
print(acc(z_lin), acc(z_ker), [acc(z_ker[:, [j]]) for j in range(3)])   # 0.502 1.0 [0.53, 0.535, 1.0]

# neighbourhoods: share of each image's 10 nearest neighbours in 64-D that stay among its 10 nearest in 2-D
def kept(A, B, k=10):
    ia = NearestNeighbors(n_neighbors=k + 1).fit(A).kneighbors(A, return_distance=False)[:, 1:]
    ib = NearestNeighbors(n_neighbors=k + 1).fit(B).kneighbors(B, return_distance=False)[:, 1:]
    return np.mean([len(set(a) & set(b)) / k for a, b in zip(ia, ib)])
for name, E in (('pca', PCA(2).fit_transform(D)), ('t-sne', TSNE(2, perplexity=30, random_state=0).fit_transform(D))):
    print(name, round(kept(D, E), 3))                               # pca 0.118 | t-sne 0.584

Next#

Probabilistic PCA and the Gaussian mixture are densities: they say not only where the data are, but how probable each point is. A point that the model finds very improbable, or that lies far from the subspace where the data live, is unusual. The next note turns that observation into methods for finding anomalies: density scores, distances to neighbours, random trees that isolate odd points quickly, and one-class support vector machines, and asks how to evaluate a detector when anomalies are rare.

References#

• K. Pearson. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559–572, 1901.
• H. Hotelling. Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24(6), 417–441, 1933.
• C. Eckart, G. Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3), 211–218, 1936.
• M. E. Tipping, C. M. Bishop. Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society B, 61(3), 611–622, 1999.
• B. Schölkopf, A. Smola, K.-R. Müller. Nonlinear Component Analysis as a Kernel Eigenvalue Problem. Neural Computation, 10(5), 1299–1319, 1998.
• L. van der Maaten, G. Hinton. Visualizing Data using t-SNE. JMLR, 9, 2579–2605, 2008.
• L. McInnes, J. Healy, J. Melville. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv:1802.03426, 2018.
• M. Wattenberg, F. Viégas, I. Johnson. How to Use t-SNE Effectively. Distill, 2016.
• P. Baldi, K. Hornik. Neural networks and principal component analysis: Learning from examples without local minima. Neural Networks, 2(1), 53–58, 1989.
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 12.1 (PCA, both formulations), Section 12.2 (probabilistic PCA; 12.2.1, the maximum likelihood solution), Section 12.3 (kernel PCA).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Chapter 10 (dimensionality reduction with PCA), Section 4.5 (the SVD), Section 4.6 (matrix approximation).
• E. Sokolov et al. Машинное обучение (HSE): lecture notes, lecture 12: unsupervised learning (2021); problems on PCA (2025).