Step-by-step tutorial ·

PCA with small pixel images

New to SVD? Start with SVD, explained with a small picture. It introduces the idea without matrix equations before you use SVD in step 5.

An image is an array of numbers. I use small images so that we can follow every step from the pixel values to the principal components and the reconstructed image.

Allow 45–60 minutes. You need basic Python and array indexing. I explain the new quantities as we use them. Start with the interactive example, then run the code cells in order.

Download the notebook · Download the Python script · Read the PCA introduction

In Jupyter or Colab, open the downloaded notebook and run its cells from top to bottom. For a local Python installation, install the packages with python -m pip install numpy matplotlib, then run python pixel-pca.py. No external images are required.

Reconstruct a pixel image

Each image has 6 × 6 pixels. The 24 images vary in left/right and top/bottom brightness contrast. They are synthetic images, with dimensionless brightness from 0 to 1.

Loading the pixel data…

Original and reconstruction: black = 0, white = 1, using the same scale. Residual: blue is negative, white is zero, red is positive, on a fixed scale from −0.4 to +0.4. The residual is original brightness minus reconstructed brightness. Pixel indices run from 0 to 5.

Show the mean image and principal-component maps

The mean image uses the brightness scale 0–1. Component maps use blue for negative weights and red for positive weights, on the fixed scale −1/6 to +1/6. A component and its negative describe the same direction; reversing its scores too leaves reconstruction unchanged.

The first component explains 88.7% of the variance and the second 11.3%. Select zero components to see the mean image; select two to recover the original. This exact recovery follows from the way I created these images. It is not a general result for microscope images.

1. Create one small image

A grayscale image stores one brightness value per pixel. I use values from 0 (black) to 1 (white), with no physical units. The first index selects a row, the second a column. Rows increase downward in this display. These are invented patterns, not microscope measurements.

import numpy as np
import matplotlib.pyplot as plt

image = np.array([
    [0.0, 0.0, 1.0, 1.0],
    [0.0, 0.0, 1.0, 1.0],
    [0.0, 0.0, 1.0, 1.0],
    [0.0, 0.0, 1.0, 1.0],
])
print("Image shape:", image.shape)
print("Pixel at row 0, column 2:", image[0, 2])
plt.imshow(image, cmap="gray", vmin=0, vmax=1, interpolation="nearest")
plt.xlabel("Column index")
plt.ylabel("Row index")
plt.colorbar(label="Brightness (dimensionless)")
plt.show()
assert image.shape == (4, 4)
assert image[0, 2] == 1.0

Check: The shape is (4, 4), and pixel [0, 2] has brightness 1.0.

2. Create a collection of images

PCA needs observations to compare. Here, one complete image is one observation, and each pixel position is one feature. I create 24 images, each with 6 rows and 6 columns. Their brightness changes along two independent patterns: a left/right contrast and a top/bottom contrast. The amplitudes are dimensionless.

size = 6
sample_count = 24
left_right = np.ones((size, size))
left_right[:, :3] = -1

top_bottom = np.ones((size, size))
top_bottom[:3, :] = -1

images = []
for sample in range(sample_count):
    angle = 2 * np.pi * sample / sample_count
    a = 0.28 * np.cos(angle)
    b = 0.10 * np.sin(angle)
    new_image = 0.5 + a * left_right + b * top_bottom
    images.append(new_image)

images = np.array(images)
print("Collection shape:", images.shape)
assert images.shape == (24, 6, 6)
assert images.min() >= 0 and images.max() <= 1

Check: The collection has shape (24, 6, 6). Every brightness value lies between 0 and 1.

3. Put each image in one data row

Flattening reads the pixels row by row. It does not average them or throw them away. The data matrix X has 24 rows (images) and 36 columns (pixel positions). The same column must refer to the same pixel position in every image. Misaligned images would change the meaning of PCA.

X = images.reshape(sample_count, size * size)
print("Data matrix shape:", X.shape)
restored = X[0].reshape(size, size)
assert X.shape == (24, 36)
assert np.array_equal(restored, images[0])

Check: Reshaping the first row back into 6 × 6 recovers the first image exactly.

4. Subtract the mean image

The mean of a feature is its average over all 24 images. I calculate a separate mean for each of the 36 pixel positions. Subtracting these means centers the data. The centered entries can be negative: they represent brightness below the mean, not an invalid negative brightness. I do not standardize the pixels here. All pixels use the same brightness scale; standardization would give a different analysis.

mean_row = X.mean(axis=0)
X_centered = X - mean_row
mean_image = mean_row.reshape(size, size)
print("Mean brightness:", mean_image.mean())
assert np.allclose(mean_image, 0.5)
assert np.allclose(X_centered.mean(axis=0), 0.0)

Check: The mean image is uniformly 0.5. Every centered pixel column has mean zero within rounding error.

5. Find the principal components

Singular value decomposition (SVD) is a way to calculate PCA for centered data. The rows of Vt give principal directions in pixel space. A direction has 36 entries, so I reshape it into a component map. A component map shows a pattern of change, not an ordinary grayscale image. Its positive and negative entries are dimensionless weights. Scores tell us how strongly each image uses that direction. The array U contains coefficients that become scores when multiplied by the singular values. The nonnegative singular values measure the strengths of the independent variations. Vt stores the principal directions as rows. Squaring a singular value and dividing by 23 gives its sample variance. A variance fraction is that variance divided by the total variance.

U, singular_values, Vt = np.linalg.svd(X_centered, full_matrices=False)
components = Vt
scores = X_centered @ components.T

variances = singular_values**2 / (sample_count - 1)
variance_fraction = variances / variances.sum()
print("First two variance fractions:", variance_fraction[:2])
assert np.allclose(variance_fraction[:2], [0.28**2, 0.10**2] / np.array(0.28**2 + 0.10**2))
assert np.allclose(variance_fraction[:2].sum(), 1.0)
assert np.allclose(components @ components.T, np.eye(sample_count))

Check: The first component explains about 88.7% of the variance; the second explains about 11.3%. All remaining variance is numerical rounding in this constructed example.

6. Reconstruct an image

For k retained components, I multiply the first k scores by the first k component directions, then add the mean image. Here, k is the number of retained components and sample is the image index, from 0 to 23. RMSE is the square root of the mean squared difference over the 36 pixels; it measures the brightness error. With k = 0, every reconstruction is the mean image. One component keeps the dominant left/right variation. Two components reproduce these particular images, because I constructed them using only two independent patterns. Real micrographs generally need more.

sample = 3
k = 1
reconstructed_row = mean_row + scores[sample, :k] @ components[:k]
reconstructed = reconstructed_row.reshape(size, size)
residual = images[sample] - reconstructed
rmse = np.sqrt(np.mean(residual**2))
print("Brightness RMSE:", rmse)

fig, axes = plt.subplots(1, 3, figsize=(10, 3))
axes[0].imshow(images[sample], cmap="gray", vmin=0, vmax=1)
axes[1].imshow(reconstructed, cmap="gray", vmin=0, vmax=1)
axes[2].imshow(residual, cmap="RdBu_r", vmin=-0.4, vmax=0.4)
for ax, title in zip(axes, ["Original", "Reconstruction", "Residual"]):
    ax.set_title(title)
    ax.set_xlabel("Column index")
    ax.set_ylabel("Row index")
plt.tight_layout()
plt.show()

full_reconstruction = mean_row + scores[:, :2] @ components[:2]
assert np.allclose(full_reconstruction, X)
assert np.isclose(rmse, 0.10 / np.sqrt(2))
print("All tutorial checks passed.")

Check: For sample 3 and one component, RMSE is about 0.07071 in brightness units. Two-component reconstruction matches every image within floating-point rounding.

Try these exercises

  1. Keep zero components and change the image index. Explain why the reconstruction stays the same.
  2. Keep one component. Find an image with a large residual and one with a nearly zero residual.
  3. Keep two components. Compare the residual and RMSE with the one-component result.
  4. In the notebook, change the second contrast amplitude from 0.10 to 0.28. Re-run from step 2 and update the fixed checks. What happens to the two variance fractions?
Worked solutions and self-checks

Zero components gives only the mean image, so changing the selected image cannot change it. With one component, image 6 has RMSE 0.10, while image 0 has a residual near zero. With two components, RMSE is near machine precision. Equal contrast amplitudes give two equal variances. There is then no unique first direction in the two-dimensional variation space.

What changes for real images?

Keep image size, alignment and intensity calibration consistent. Decide whether brightness, contrast or spatial arrangement is the variation you want to study. A component with large variance can describe illumination changes rather than a material feature. Color images can be represented by red, green and blue channels, but flattening does not give PCA knowledge of neighboring pixels or interfaces.

If PCA is part of a prediction model, fit the mean and components using the training images only. Apply that fitted transformation to held-out images. Here I analyze the whole synthetic collection to explain PCA; I do not claim a prediction test.