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
# I store the brightness values as rows of pixels.
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.0Check: The shape is (4, 4), and pixel [0, 2] has brightness 1.0.
Read the code, line by line
np.array collects the four lists of pixel values into a two-dimensional array. Each inner list is one row. Python counts from zero, so image[0, 2] means the first row and the third column. image.shape reports the number of rows and columns.
imshow displays the numbers as an image. cmap="gray" selects a grayscale display. I fix vmin=0 and vmax=1 so that black and white mean the same thing in every image. interpolation="nearest" keeps the individual pixels visible instead of smoothing their edges. The color bar labels the brightness scale.
An assert is a small check. If its condition is false, Python stops with an AssertionError. I use it to check the result before moving on.
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))
# I change the first three columns to make a left/right contrast.
left_right[:, :3] = -1
top_bottom = np.ones((size, size))
top_bottom[:3, :] = -1
# I build the collection one image at a time.
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() <= 1Check: The collection has shape (24, 6, 6). Every brightness value lies between 0 and 1.
Read the code, line by line
np.ones((size, size)) makes a 6 × 6 array filled with 1. The slice [:, :3] selects every row and the first three columns. Replacing those entries with −1 creates a left/right pattern. The slice [:3, :] selects the first three rows and every column, giving a top/bottom pattern.
The loop runs 24 times. sample takes the values 0 through 23. The angle is measured in radians; the loop samples a complete cycle at equally spaced angles. The two amplitudes, a and b, change the strength of the two patterns. I add both patterns to a mean brightness of 0.5.
images.append(new_image) adds one image to the list. After the loop, np.array(images) collects them into one array. The shape (24, 6, 6) means 24 images, six rows per image and six columns per image. The bounds check confirms that all constructed brightness values remain 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.
# I keep all pixel values and arrange each image in one row.
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.
Read the code, line by line
PCA expects a table with observations in rows and features in columns. reshape(sample_count, size * size) changes the arrangement of the numbers without changing the numbers themselves. The default order reads each image row by row.
There are 24 observations and 36 features. I then reshape the first data row back into a 6 × 6 image. np.array_equal checks that it is identical to the first original image. This confirms that the conversion is reversible.
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.
# I calculate one mean for each pixel position.
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.
Read the code, line by line
X.mean(axis=0) averages down the image rows of the data table. It produces one mean for each pixel position, so mean_row contains 36 values. This is different from taking one average over the entire table.
When I subtract mean_row from X, NumPy subtracts the same 36 means from each of the 24 observations. mean_image arranges those means into an image so I can inspect them.
np.allclose allows small floating-point rounding differences. I use it rather than exact equality for calculations involving averages and decimal numbers.
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
# I project each centered image onto the component directions.
scores = X_centered @ components.T
# I calculate the sample variance along each principal direction.
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.
Read the code, line by line
np.linalg.svd returns three arrays. Here, U has shape (24, 24), singular_values has length 24 and Vt has shape (24, 36). The option full_matrices=False avoids storing additional directions that are unnecessary for this calculation.
Each row of components is a direction in the 36-feature pixel space. .T exchanges rows and columns. The operator @ performs matrix multiplication. Thus X_centered @ components.T calculates a score for each observation along each direction.
The operator **2 squares each singular value. I divide by sample_count - 1, which is 23, to obtain sample variances. Dividing each variance by their sum gives the variance fractions. [:2] selects the first two entries.
The final check uses np.eye(sample_count), an identity matrix with ones on its diagonal and zeros elsewhere. It checks that the component directions are normalized and mutually perpendicular. It does not show that these directions identify physical mechanisms.
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
# I add the retained patterns to the mean image.
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()
# I check that both constructed patterns recover the complete collection.
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.
Read the code, line by line
sample = 3 selects the fourth image because indexing begins at zero. k = 1 keeps only the first component. scores[sample, :k] selects that image's first score, and components[:k] selects the corresponding component direction. Their matrix product gives the centered approximation. Adding mean_row restores the mean brightness.
The residual is the original image minus its reconstruction, pixel by pixel. residual**2 squares each difference, np.mean averages the 36 squared differences and np.sqrt takes the square root. This gives the RMSE.
plt.subplots(1, 3) creates a figure with one row of three plotting areas. axes[0], axes[1] and axes[2] refer to these areas. I use the same 0–1 scale for the original and reconstructed images. The residual uses a separate signed scale from −0.4 to +0.4.
zip pairs each plotting area with its title. The loop labels all three areas. tight_layout adjusts the spacing, and show displays the figure. The final checks confirm that two components recover all 24 images and that the selected one-component example has the expected error.
Try these exercises
- Keep zero components and change the image index. Explain why the reconstruction stays the same.
- Keep one component. Find an image with a large residual and one with a nearly zero residual.
- Keep two components. Compare the residual and RMSE with the one-component result.
- 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.
Next: Generate colored microstructures and play the PCA reconstruction game.