Seeing high-dimensional data in two dimensions

In post 15 we built the digits loading and inspection pipeline, and we ended with a matrix of 60,000 images, each one a 28 by 28 grid of pixels. Flatten that grid and every image becomes a point in 784 dimensions. This post takes those points and asks a single question: what does the data look like when we force it onto a flat sheet of paper? We will fit eleven different methods, from plain Principal Component Analysis to UMAP, and we will score every one of them the same way. The headline lands early: the best map, t-SNE, keeps 42.5 percent of the ten nearest neighbour links that exist in the full 784 dimensions, against 3.9 percent for a random 2D projection and 0.67 percent by chance. That gap is the whole lesson.

The through-line for this post is the shadow. A 2D map of 784-dimensional data is a shadow cast on a wall, and like any shadow it keeps the outline of the object while losing the depth. Some of our methods cast the shadow with a straight light source, some bend the light around the object, and some trace the surface of the object itself. The question we keep returning to is simple: which shadow still tells us who is standing there?

The data and what we found

To test which shadow works, we need the object casting it: the MNIST digits, a standard set of handwritten digit images, 60,000 for training and 10,000 for testing, each one 784 pixels of grayscale ink. Before fitting anything we ran the exploratory pass and wrote down what it told us, because those findings fix every choice that follows.

The raw pixels span the full 0 to 255 range and there are no missing values in either split, so nothing gets imputed and every method sees the same matrix. That matters because when we compare eleven methods, we want the only thing changing to be the method.

Two findings shaped the rest of the post. First, the class balance is near uniform, with the most common digit, the 1, holding 6,742 of 60,000 images, so the majority class base rate, the accuracy of always guessing the most common digit, sits at 0.1124. Any accuracy we report in two dimensions has to be read against that line, not against zero. Second, pixel space is ambiguous on its own. We took a 2,000-image probe set and asked how often an image's nearest neighbour in 784 dimensions carries a different label.

167 of those 2,000 images, 8.3 percent, have a nearest pixel neighbour with a different digit. That is a ceiling. No 2D map can recover neighbour links that were never there in the first place, so the top of our leaderboard is a ceiling rather than a failure. We also hashed every image, computing a short fingerprint that lets us spot exact duplicates, found none to drop, and confirmed that 80.9 percent of pixels are exactly zero, which means the useful variance lives in a minority of pixels and the mean image is a strong baseline to beat.

Linear and kernel methods

The data gave us a baseline; the first way to cast the shadow is linear. A linear map writes each 784-dimensional image as a weighted sum of a few directions. Principal Component Analysis picks the directions with the largest variance. Explained variance tells us how much of the total spread those directions keep. We fit PCA on a fixed pool of 6,000 images, scaled to 0 to 1, and read the curve.

The first component keeps 9.9 percent of the variance and the first two together keep 16.8 percent. To reach 90 percent of the variance we need 86 components. That number is the honest dimensionality of this data, and it is far below 784, which is exactly why MNIST is easy to compress and hard to draw in two dimensions. The first two components carry a small slice of the spread, so the 2D PCA picture is a shadow of a shadow.

The reconstruction snippet prints the MSE, the mean squared error, at each component count. The mean image gives a baseline MSE of 0.06724, and PCA has to beat it.

# Baseline first: the mean image is the best constant reconstruction, so PCA has to beat that MSE.
mean_mse = float(np.mean((X_pool - X_pool.mean(axis=0)) ** 2))
Z_pool = pca.transform(X_pool)
for k in (2, 10, 50, 200):
    rec = pca.mean_ + Z_pool[:, :k] @ pca.components_[:k]
    print(k, float(np.mean((X_pool - rec) ** 2)))