PCA for Images
STEP 1 · OVERVIEW

Implement PCA Dimensionality Reduction for Images

In this practical codelab you will write every line of Python code needed to load images, apply PCA, visualise eigenfaces, and reconstruct compressed images — step by step.

By the end of this lab you will be able to:
  • Load and preprocess a set of face images into a data matrix
  • Standardise the matrix and apply PCA using scikit-learn
  • Inspect explained variance and choose the right number of components
  • Visualise the mean face and eigenfaces
  • Reconstruct compressed images and measure quality
  • Plot a scree chart and compare compression levels

What you will build

Raw face images (Olivetti / ORL dataset)
↓ flatten each image to a row vector
↓ stack into matrix X  (n_samples × n_pixels)
↓ standardise
↓ apply PCA → choose k components
↓ reconstruct images from k components
↓ compute reconstruction error
↓ visualise eigenfaces and compressed images
🐍 Language

Python 3.8+

📦 Libraries

NumPy, scikit-learn, Matplotlib

📂 Dataset

Olivetti Faces — built into scikit-learn, no download needed

⏱ Estimated time

45 – 60 minutes

Prerequisite

This lab follows the conceptual codelab Dimensionality Reduction (Module 4.2). Make sure you understand what PCA does before writing the code.

STEP 2 · SETUP

Set Up Your Environment

Create a clean Python environment and install the required packages before writing any code.

Option A — local virtual environment

terminal
python -m venv pca_lab
source pca_lab/bin/activate       # Windows: pca_lab\Scripts\activate
pip install numpy scikit-learn matplotlib scikit-image

Option B — Google Colab / Jupyter

notebook cell
!pip install -q scikit-image   # numpy, sklearn, matplotlib are pre-installed

Verify your setup

pca_images.py
import numpy as np
import sklearn
import matplotlib
print("NumPy      :", np.__version__)
print("scikit-learn:", sklearn.__version__)
print("Matplotlib :", matplotlib.__version__)
print("Setup OK ✓")
expected output NumPy : 1.26.x scikit-learn: 1.4.x Matplotlib : 3.8.x Setup OK ✓
STEP 3 · DATA

Load the Olivetti Faces Dataset

The Olivetti Faces dataset ships with scikit-learn — 400 grayscale 64×64 face images of 40 different people, 10 images each.

Write the import and loading code

pca_images.py — Part 1: imports & data
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_olivetti_faces
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# Load dataset
dataset = fetch_olivetti_faces(shuffle=True, random_state=42)
X = dataset.data          # shape: (400, 4096)
y = dataset.target        # shape: (400,) person identity 0-39

n_samples, n_features = X.shape
image_shape = (64, 64)

print(f"Dataset loaded:")
print(f"  Samples  : {n_samples}")
print(f"  Features : {n_features}  (64 x 64 pixels)")
print(f"  Classes  : {len(np.unique(y))} people")
expected output Dataset loaded: Samples : 400 Features : 4096 (64 x 64 pixels) Classes : 40 people
KEY POINT

Each image is flattened into a row vector of 4,096 numbers (64×64). The full dataset matrix X is therefore 400×4,096 — a very wide table. PCA will reduce the 4,096 columns to a much smaller number.

Visualise a sample of raw images

pca_images.py — Part 1 continued
fig, axes = plt.subplots(2, 5, figsize=(12, 5))
fig.suptitle("Sample raw images from Olivetti Faces", fontsize=14, fontweight="bold")
for ax, image, label in zip(axes.flat, X[:10], y[:10]):
    ax.imshow(image.reshape(image_shape), cmap="gray")
    ax.set_title(f"Person {label}")
    ax.axis("off")
plt.tight_layout()
plt.savefig("01_raw_samples.png", dpi=120, bbox_inches="tight")
plt.show()
What to check

You should see 10 small grayscale face images. Note how lighting, expression, and pose vary across people.

STEP 4 · PREPROCESSING

Standardise the Image Matrix

PCA is sensitive to scale. We standardise each pixel feature so that no pixel dominates just because of its brightness range.

Why standardise?

Without standardisation

PCA directions are dominated by features with the largest numerical range, not the most informative ones.

With standardisation

Every pixel has mean 0 and variance 1, so PCA looks for the most informative directions regardless of raw scale.

pca_images.py — Part 2: standardisation
# Standardise: zero mean, unit variance per pixel
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("Scaled matrix shape:", X_scaled.shape)
print("Mean of first feature (should be ~0):", round(X_scaled[:, 0].mean(), 4))
print("Std  of first feature (should be ~1):", round(X_scaled[:, 0].std(), 4))
expected output Scaled matrix shape: (400, 4096) Mean of first feature (should be ~0): 0.0 Std of first feature (should be ~1): 1.0
Raw X  (400 × 4096, pixel range 0.0 – 1.0)
↓ StandardScaler().fit_transform(X)
X_scaled  (400 × 4096, mean≈0, std≈1)

Why should we standardise before applying PCA?

STEP 5 · PCA

Apply PCA to the Image Matrix

Use scikit-learn’s PCA class to find the principal components of the face dataset, keeping 150 components to start.

pca_images.py — Part 3: apply PCA
N_COMPONENTS = 150   # we will experiment with this later

pca = PCA(n_components=N_COMPONENTS, whiten=True, random_state=42)
X_pca = pca.fit_transform(X_scaled)

print(f"PCA applied:")
print(f"  Original  : {X_scaled.shape}")
print(f"  Compressed: {X_pca.shape}")
print(f"  Ratio     : {X_pca.shape[1]/X_scaled.shape[1]:.2%} of original features")

cum_var = np.cumsum(pca.explained_variance_ratio_)
print(f"\nVariance explained by top {N_COMPONENTS} components: {cum_var[-1]:.2%}")
expected output PCA applied: Original : (400, 4096) Compressed: (400, 150) Ratio : 3.66% of original features Variance explained by top 150 components: 99.xx%

Understand the output

VariableShapeMeaning
X_scaled(400, 4096)Standardised full image matrix
X_pca(400, 150)Compressed representation
pca.components_(150, 4096)Eigenfaces (principal directions)
pca.explained_variance_ratio_(150,)Fraction of variance each PC captures
pca.mean_(4096,)Mean face (per-pixel average)
PARAMETER

whiten=True divides each component by its standard deviation. This gives components with unit variance and often helps downstream classifiers.

STEP 6 · ANALYSIS

Plot Explained Variance (Scree Chart)

The scree plot shows how much information each principal component carries, helping us decide how many components to keep.

pca_images.py — Part 4: scree plot
# Fit full PCA to see ALL components
pca_full = PCA(random_state=42)
pca_full.fit(X_scaled)

ev_ratio = pca_full.explained_variance_ratio_
cum_ev   = np.cumsum(ev_ratio)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
fig.suptitle("Explained Variance -- Olivetti Faces", fontsize=14, fontweight="bold")

ax1.bar(range(1, 51), ev_ratio[:50], color="#2856d8", edgecolor="#101010", linewidth=0.8)
ax1.set_xlabel("Principal Component")
ax1.set_ylabel("Explained Variance Ratio")
ax1.set_title("Individual (first 50 components)")

ax2.plot(range(1, len(cum_ev)+1), cum_ev, color="#2856d8", linewidth=2)
ax2.axhline(y=0.90, color="#ff9bb7", linestyle="--", linewidth=1.5, label="90%")
ax2.axhline(y=0.95, color="#ffd83d", linestyle="--", linewidth=1.5, label="95%")
ax2.axhline(y=0.99, color="#75e6c5", linestyle="--", linewidth=1.5, label="99%")
ax2.set_xlabel("Number of Components")
ax2.set_ylabel("Cumulative Explained Variance")
ax2.set_title("Cumulative")
ax2.legend()
ax2.set_xlim(0, 200)

plt.tight_layout()
plt.savefig("02_scree_plot.png", dpi=120, bbox_inches="tight")
plt.show()

for threshold in [0.80, 0.90, 0.95, 0.99]:
    n = np.argmax(cum_ev >= threshold) + 1
    print(f"  {threshold:.0%} variance -> {n} components  ({n/X_scaled.shape[1]:.2%} of features)")
Reading the chart
  • The first few bars are tall — those components explain a lot.
  • The curve flattens quickly — a small number of components carries most information.
  • This is the power of PCA: enormous compression with minimal information loss.

How do you use the cumulative explained variance curve to choose the number of components?

STEP 7 · VISUALISE

Visualise the Mean Face and Eigenfaces

The mean face is the average of all training images. Eigenfaces are the principal directions — they look like ghostly faces capturing the main patterns of variation.

Plot the mean face

pca_images.py — Part 5: mean face
mean_face = pca.mean_.reshape(image_shape)

plt.figure(figsize=(3, 3))
plt.imshow(mean_face, cmap="gray")
plt.title("Mean Face", fontsize=13, fontweight="bold")
plt.axis("off")
plt.tight_layout()
plt.savefig("03_mean_face.png", dpi=120, bbox_inches="tight")
plt.show()
What is the mean face?

It is the pixel-wise average of all 400 training images. PCA then finds the directions of variation around this mean.

Plot the top eigenfaces

pca_images.py — Part 5 continued: eigenfaces
fig, axes = plt.subplots(3, 5, figsize=(14, 9))
fig.suptitle("Top Eigenfaces (principal components)", fontsize=14, fontweight="bold")

for idx, ax in enumerate(axes.flat):
    eigenface = pca.components_[idx].reshape(image_shape)
    # Normalise each eigenface to [0,1] for display
    eigenface_display = (eigenface - eigenface.min()) / (eigenface.max() - eigenface.min())
    ax.imshow(eigenface_display, cmap="gray")
    ax.set_title(f"PC {idx+1}\n{pca.explained_variance_ratio_[idx]:.2%}", fontsize=8)
    ax.axis("off")

plt.tight_layout()
plt.savefig("04_eigenfaces.png", dpi=120, bbox_inches="tight")
plt.show()
What do eigenfaces look like?

They look like ghostly, abstract face patterns. PC1 captures the most common variation (e.g. overall lighting). Later PCs capture finer details like eye shape or hair texture.

INSIGHT

Each original face image can be expressed as a weighted sum of eigenfaces plus the mean face. PCA compression means keeping only the weights for the top-k eigenfaces.

STEP 8 · RECONSTRUCTION

Reconstruct Images from Compressed Representations

After compressing with PCA, we can reconstruct approximate versions of the original images. Quality improves as we keep more components.

pca_images.py — Part 6: reconstruction
def reconstruct(X_scaled, n_components, image_shape=(64, 64)):
    """Compress X with n_components PCA dimensions and reconstruct."""
    pca_k   = PCA(n_components=n_components, whiten=True, random_state=42)
    X_proj  = pca_k.fit_transform(X_scaled)
    X_recon = pca_k.inverse_transform(X_proj)
    X_recon = scaler.inverse_transform(X_recon)
    X_recon = np.clip(X_recon, 0, 1)
    return X_recon

levels = [1, 5, 10, 20, 50, 100, 150, 200]
sample_idx = 0
original_face = X[sample_idx].reshape(image_shape)

fig, axes = plt.subplots(2, 4, figsize=(14, 7))
fig.suptitle("Image Reconstruction at Different Compression Levels", fontsize=13, fontweight="bold")

for ax, k in zip(axes.flat, levels):
    recon      = reconstruct(X_scaled, n_components=k)
    recon_face = recon[sample_idx].reshape(image_shape)
    mse        = np.mean((original_face - recon_face) ** 2)
    ax.imshow(recon_face, cmap="gray", vmin=0, vmax=1)
    ax.set_title(f"k={k}\nMSE={mse:.4f}", fontsize=9)
    ax.axis("off")

plt.tight_layout()
plt.savefig("05_reconstruction.png", dpi=120, bbox_inches="tight")
plt.show()
what to observe k=1 -> very blurry, only overall brightness captured k=5 -> rough face shape visible k=20 -> recognisable face, some smoothing k=100 -> nearly identical to original k=200 -> very close to lossless
MSE (Mean Squared Error)

MSE measures the average squared difference between original and reconstructed pixel values. Lower MSE = better reconstruction. As k increases, MSE decreases.

STEP 9 · ANALYSIS

Plot Reconstruction Error vs Number of Components

Quantify how reconstruction quality improves as we add more principal components, and find the sweet spot between compression and accuracy.

pca_images.py — Part 7: MSE curve
k_values  = list(range(1, 21)) + list(range(25, 201, 5))
mse_values = []
X_original = np.clip(scaler.inverse_transform(X_scaled), 0, 1)

for k in k_values:
    X_recon_k = reconstruct(X_scaled, n_components=k)
    mse = np.mean((X_original - X_recon_k) ** 2)
    mse_values.append(mse)

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(k_values, mse_values, color="#2856d8", linewidth=2, marker=".", markersize=4)
ax.fill_between(k_values, mse_values, alpha=0.15, color="#2856d8")
ax.set_xlabel("Number of PCA Components (k)", fontsize=12)
ax.set_ylabel("Mean Squared Error", fontsize=12)
ax.set_title("Reconstruction Error vs. Number of Components", fontsize=13, fontweight="bold")
ax.axvline(x=50,  color="#ffd83d", linestyle="--", linewidth=1.5, label="k=50")
ax.axvline(x=100, color="#75e6c5", linestyle="--", linewidth=1.5, label="k=100")
ax.legend()
ax.set_xlim(0, max(k_values))
plt.tight_layout()
plt.savefig("06_mse_curve.png", dpi=120, bbox_inches="tight")
plt.show()
Interpretation
  • The curve drops steeply at first — the first few components reduce error dramatically.
  • After ~50–100 components, gains become much smaller (diminishing returns).
  • This matches the scree plot from Step 6.

What happens to reconstruction MSE as we increase k (number of PCA components)?

STEP 10 · VISUALISE

Original vs Reconstructed — Gallery

Create a polished comparison showing original images next to their reconstructions at k = 50 and k = 150 components.

pca_images.py — Part 8: comparison gallery
X_recon_50  = reconstruct(X_scaled, n_components=50)
X_recon_150 = reconstruct(X_scaled, n_components=150)
X_original_unscaled = np.clip(scaler.inverse_transform(X_scaled), 0, 1)

n_display = 6
fig, axes = plt.subplots(3, n_display, figsize=(16, 9))
fig.suptitle("Original  |  k=50 components  |  k=150 components",
             fontsize=13, fontweight="bold")

for col in range(n_display):
    idx = col * 5
    axes[0, col].imshow(X_original_unscaled[idx].reshape(image_shape), cmap="gray", vmin=0, vmax=1)
    axes[0, col].set_title(f"Original\n(Person {y[idx]})", fontsize=8)
    axes[0, col].axis("off")

    axes[1, col].imshow(X_recon_50[idx].reshape(image_shape), cmap="gray", vmin=0, vmax=1)
    mse_50 = np.mean((X_original_unscaled[idx] - X_recon_50[idx]) ** 2)
    axes[1, col].set_title(f"k=50\nMSE={mse_50:.4f}", fontsize=8)
    axes[1, col].axis("off")

    axes[2, col].imshow(X_recon_150[idx].reshape(image_shape), cmap="gray", vmin=0, vmax=1)
    mse_150 = np.mean((X_original_unscaled[idx] - X_recon_150[idx]) ** 2)
    axes[2, col].set_title(f"k=150\nMSE={mse_150:.4f}", fontsize=8)
    axes[2, col].axis("off")

plt.tight_layout()
plt.savefig("07_comparison_gallery.png", dpi=120, bbox_inches="tight")
plt.show()
What you will see

Row 1: originals. Row 2 (k=50): blurry but recognisable. Row 3 (k=150): nearly identical to row 1. Yet k=150 uses only 150/4096 ≈ 3.7% of the original storage per image.

STEP 11 · ANALYSIS

Measure Compression Ratio

Quantify how much space is saved by PCA compression compared to storing the full images.

pca_images.py — Part 9: compression analysis
original_size = n_samples * n_features  # 400 x 4096 = 1,638,400

print("=" * 60)
print(f"{'Compression Summary':^60}")
print("=" * 60)
print(f"{'Components':>12}  {'Size (floats)':>14}  {'Saved':>7}  {'Var. explained':>16}")
print("-" * 60)

for k in [5, 10, 20, 50, 100, 150, 200]:
    compressed_size = n_samples * k + k * n_features
    saving = 1 - compressed_size / original_size
    pca_k  = PCA(n_components=k, random_state=42).fit(X_scaled)
    var_exp = np.sum(pca_k.explained_variance_ratio_)
    print(f"{k:>12}  {compressed_size:>14,}  {saving:>6.1%}  {var_exp:>15.2%}")

print("=" * 60)
print(f"{'Original':>12}  {original_size:>14,}  {'--':>7}  {'100.00%':>16}")
Trade-off

k = 50 saves ~87% of storage while preserving most recognisable image quality. This is the fundamental power and trade-off of PCA compression.

STEP 12 · COMPLETE CODE

The Complete Script

Below is the full, self-contained Python script containing all parts from Steps 2–11. Copy it into pca_images.py and run it end-to-end.

pca_images.py — COMPLETE SCRIPT
"""
pca_images.py
Implement Dimensionality Reduction using PCA for Images
Practical Codelab -- Computer Vision
"""

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_olivetti_faces
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# 1. Load dataset
dataset = fetch_olivetti_faces(shuffle=True, random_state=42)
X = dataset.data          # (400, 4096)
y = dataset.target        # (400,)
n_samples, n_features = X.shape
image_shape = (64, 64)
print(f"Loaded: {n_samples} images, {n_features} features")

# 2. Visualise samples
fig, axes = plt.subplots(2, 5, figsize=(12, 5))
fig.suptitle("Sample Raw Images", fontsize=13, fontweight="bold")
for ax, img, lbl in zip(axes.flat, X[:10], y[:10]):
    ax.imshow(img.reshape(image_shape), cmap="gray"); ax.set_title(f"P{lbl}"); ax.axis("off")
plt.tight_layout(); plt.savefig("01_raw_samples.png", dpi=120, bbox_inches="tight"); plt.show()

# 3. Standardise
scaler   = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 4. Apply PCA (150 components)
N_COMPONENTS = 150
pca   = PCA(n_components=N_COMPONENTS, whiten=True, random_state=42)
X_pca = pca.fit_transform(X_scaled)
print(f"Compressed: {X_scaled.shape} -> {X_pca.shape}")
print(f"Variance explained: {np.sum(pca.explained_variance_ratio_):.2%}")

# 5. Scree plot
pca_full = PCA(random_state=42).fit(X_scaled)
ev       = pca_full.explained_variance_ratio_
cum_ev   = np.cumsum(ev)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
fig.suptitle("Explained Variance", fontsize=13, fontweight="bold")
ax1.bar(range(1, 51), ev[:50], color="#2856d8", edgecolor="#101010", linewidth=0.8)
ax1.set(xlabel="Component", ylabel="Variance Ratio", title="Individual (first 50)")
ax2.plot(range(1, len(cum_ev)+1), cum_ev, color="#2856d8", linewidth=2)
for t, c in [(0.90, "#ff9bb7"), (0.95, "#ffd83d"), (0.99, "#75e6c5")]:
    ax2.axhline(t, color=c, linestyle="--", linewidth=1.5, label=f"{t:.0%}")
ax2.set(xlabel="Components", ylabel="Cumulative Variance", title="Cumulative", xlim=(0,200))
ax2.legend()
plt.tight_layout(); plt.savefig("02_scree_plot.png", dpi=120, bbox_inches="tight"); plt.show()
for t in [0.80, 0.90, 0.95, 0.99]:
    print(f"  {t:.0%} variance -> {np.argmax(cum_ev >= t)+1} components")

# 6. Mean face
plt.figure(figsize=(3, 3))
plt.imshow(pca.mean_.reshape(image_shape), cmap="gray")
plt.title("Mean Face", fontweight="bold"); plt.axis("off")
plt.tight_layout(); plt.savefig("03_mean_face.png", dpi=120, bbox_inches="tight"); plt.show()

# 7. Eigenfaces
fig, axes = plt.subplots(3, 5, figsize=(14, 9))
fig.suptitle("Top 15 Eigenfaces", fontsize=13, fontweight="bold")
for idx, ax in enumerate(axes.flat):
    ef = pca.components_[idx].reshape(image_shape)
    ef = (ef - ef.min()) / (ef.max() - ef.min())
    ax.imshow(ef, cmap="gray")
    ax.set_title(f"PC{idx+1}: {pca.explained_variance_ratio_[idx]:.2%}", fontsize=8); ax.axis("off")
plt.tight_layout(); plt.savefig("04_eigenfaces.png", dpi=120, bbox_inches="tight"); plt.show()

# 8. Reconstruction helper
def reconstruct(X_sc, k):
    pca_k = PCA(n_components=k, whiten=True, random_state=42)
    X_r   = pca_k.inverse_transform(pca_k.fit_transform(X_sc))
    return np.clip(scaler.inverse_transform(X_r), 0, 1)

# 9. Reconstruction gallery
levels = [1, 5, 10, 20, 50, 100, 150, 200]
orig   = np.clip(scaler.inverse_transform(X_scaled), 0, 1)
fig, axes = plt.subplots(2, 4, figsize=(14, 7))
fig.suptitle("Reconstruction at Different Compression Levels", fontsize=13, fontweight="bold")
for ax, k in zip(axes.flat, levels):
    rc  = reconstruct(X_scaled, k)[0].reshape(image_shape)
    mse = np.mean((orig[0].reshape(image_shape) - rc)**2)
    ax.imshow(rc, cmap="gray", vmin=0, vmax=1)
    ax.set_title(f"k={k}\nMSE={mse:.4f}", fontsize=9); ax.axis("off")
plt.tight_layout(); plt.savefig("05_reconstruction.png", dpi=120, bbox_inches="tight"); plt.show()

# 10. MSE curve
k_vals  = list(range(1, 21)) + list(range(25, 201, 5))
mse_arr = [np.mean((orig - reconstruct(X_scaled, k))**2) for k in k_vals]
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(k_vals, mse_arr, color="#2856d8", linewidth=2, marker=".", markersize=4)
ax.fill_between(k_vals, mse_arr, alpha=0.15, color="#2856d8")
ax.set(xlabel="k (components)", ylabel="MSE", title="Reconstruction Error vs. Components")
ax.axvline(50,  color="#ffd83d", linestyle="--", linewidth=1.5, label="k=50")
ax.axvline(100, color="#75e6c5", linestyle="--", linewidth=1.5, label="k=100")
ax.legend()
plt.tight_layout(); plt.savefig("06_mse_curve.png", dpi=120, bbox_inches="tight"); plt.show()

# 11. Comparison gallery
r50  = reconstruct(X_scaled, 50)
r150 = reconstruct(X_scaled, 150)
fig, axes = plt.subplots(3, 6, figsize=(16, 9))
fig.suptitle("Original  |  k=50  |  k=150", fontsize=13, fontweight="bold")
for col in range(6):
    idx = col * 5
    for row, (data, label) in enumerate([(orig,"Original"),(r50,"k=50"),(r150,"k=150")]):
        mse = "" if row==0 else f"\nMSE={np.mean((orig[idx]-data[idx])**2):.4f}"
        axes[row,col].imshow(data[idx].reshape(image_shape), cmap="gray", vmin=0, vmax=1)
        axes[row,col].set_title(f"{label}{mse}", fontsize=8); axes[row,col].axis("off")
plt.tight_layout(); plt.savefig("07_comparison_gallery.png", dpi=120, bbox_inches="tight"); plt.show()

# 12. Compression summary
orig_size = n_samples * n_features
print(f"\n{'='*60}\n{'Compression Summary':^60}\n{'='*60}")
print(f"{'k':>10}  {'Floats':>12}  {'Saved':>7}  {'Var%':>8}")
print("-"*45)
for k in [5, 10, 20, 50, 100, 150, 200]:
    sz = n_samples*k + k*n_features
    sv = 1 - sz/orig_size
    ve = np.sum(PCA(n_components=k, random_state=42).fit(X_scaled).explained_variance_ratio_)
    print(f"{k:>10}  {sz:>12,}  {sv:>6.1%}  {ve:>7.2%}")
print(f"{'Original':>10}  {orig_size:>12,}  {'--':>7}  {'100.00%':>8}")
print("\nAll outputs saved. Lab complete!")
Run the script
terminal
python pca_images.py

Seven PNG files will be saved in your working directory as the script runs.

STEP 13 · EXTEND

Extensions and Experiments

Now that the core pipeline works, try these extensions to deepen your understanding.

Extension A — Use your own images

pca_custom.py — custom image loader
import os, numpy as np
from skimage import io, color, transform

IMAGE_DIR   = "my_faces/"      # folder containing .jpg or .png images
TARGET_SIZE = (64, 64)

images = []
for fname in sorted(os.listdir(IMAGE_DIR)):
    if fname.lower().endswith((".jpg", ".png", ".jpeg")):
        img = io.imread(os.path.join(IMAGE_DIR, fname))
        if img.ndim == 3:
            img = color.rgb2gray(img)
        img = transform.resize(img, TARGET_SIZE)
        images.append(img.flatten())

X_custom = np.array(images)
print(f"Custom dataset: {X_custom.shape}")
# Then continue with StandardScaler -> PCA -> reconstruct

Extension B — PCA features for classification

pca_classification.py
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score

X_train, X_test, y_train, y_test = train_test_split(
    X_pca, y, test_size=0.25, random_state=42
)
clf = SVC(kernel="rbf", C=10, gamma="scale")
clf.fit(X_train, y_train)
print(f"Accuracy with {N_COMPONENTS} PCA components: {accuracy_score(y_test, clf.predict(X_test)):.2%}")
What to expect

PCA features often give similar or better accuracy than raw pixels, with much less computation.

🧠 Challenge Questions

  1. With k = 10 components, what percentage of variance is retained? What does the reconstruction look like?
  2. How many components are needed to exceed 95% cumulative variance on this dataset?
  3. Why does whiten=True help downstream classifiers?
  4. What would happen if you applied PCA to colour (RGB) images? How would you modify the pipeline?
  5. If the dataset grows to 5,000 images of 128×128, what changes computationally?
STEP 14 · REVISION

Lab Revision & Wrap-Up

Review the key concepts and code patterns from this lab before moving on.

Key code patterns

1. Load and flatten images into a matrix
X = dataset.data — each row is one flattened image. Shape: (n_samples, n_pixels).
2. Standardise with StandardScaler
X_scaled = StandardScaler().fit_transform(X) — always fit on training data only.
3. Apply PCA
pca = PCA(n_components=k, whiten=True); X_pca = pca.fit_transform(X_scaled) — result shape: (n_samples, k).
4. Inspect explained variance
pca.explained_variance_ratio_ per component; np.cumsum(...) for cumulative.
5. Access eigenfaces
pca.components_[i].reshape(image_shape) — pca.components_ shape: (k, n_pixels).
6. Reconstruct images
X_recon = scaler.inverse_transform(pca.inverse_transform(X_pca)) — then clip to [0,1].

Full pipeline at a glance

fetch_olivetti_faces()
↓ X  (400 × 4096)
↓ StandardScaler().fit_transform(X)
↓ X_scaled  (400 × 4096)
↓ PCA(n_components=k).fit_transform(X_scaled)
↓ X_pca  (400 × k)  ← compressed representation
↓ pca.inverse_transform(X_pca)
↓ scaler.inverse_transform(...)
↓ Reconstructed X  (400 × 4096)

Completion checklist

  • I loaded the Olivetti dataset and checked its shape.
  • I standardised the image matrix using StandardScaler.
  • I applied PCA and confirmed the compressed shape.
  • I plotted and interpreted the scree chart.
  • I visualised the mean face and at least 10 eigenfaces.
  • I reconstructed images at different k values and observed quality differences.
  • I plotted the MSE vs k curve.
  • I computed and interpreted the compression ratio table.

Lab complete! 🎉

You have implemented a full PCA pipeline for image dimensionality reduction — from raw pixels to eigenfaces to reconstruction analysis. This exact workflow is the foundation of classical face recognition systems and image compression algorithms.

← Back to the Conceptual Codelab (Module 4.2)