Chapter 11 Experiment 12: ◆ Eigenfaces—Using the SVD to “Understand” Faces

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from ipywidgets import widgets, Output, HBox, VBox, Layout
from sklearn.datasets import fetch_olivetti_faces

np.set_printoptions(precision=4, suppress=True, linewidth=100)

# ── Color scheme: book palette ──────────────────────────────────────
C_BG   = "#F8F8F8"
C_GRID = "#D6D6D6"
C_AXIS = "#000000"
C_V1   = "#57068C"    # violet
C_V2   = "#006385"    # deep blue
C_T1   = "#2AD2C9"    # teal
C_WARN = "#FF5D47"    # orange
C_AUX  = "#AB82C5"    # light violet

# Fonts: the English edition needs no CJK font, so nothing is downloaded when _lang == 'en';
# the Chinese editions use this same block to fetch Noto Sans TC/SC where no CJK font is installed
import os, urllib.request
import matplotlib.font_manager as fm
_lang = 'en'
_cjk = ['Microsoft JhengHei', 'PingFang TC', 'Noto Sans CJK TC', 'Noto Sans TC']
_have = {f.name for f in fm.fontManager.ttflist}
if _lang != 'en' and not _have & set(_cjk):
    _font = os.path.join(os.path.expanduser('~'), '.cache', 'fonts', 'NotoSansTC.ttf')
    try:
        if not os.path.exists(_font):
            os.makedirs(os.path.dirname(_font), exist_ok=True)
            urllib.request.urlretrieve('https://github.com/google/fonts/raw/main/ofl/notosanstc/NotoSansTC%5Bwght%5D.ttf', _font + '.part')
            os.replace(_font + '.part', _font)
        fm.fontManager.addfont(_font)
        _have.add('Noto Sans TC')
    except OSError as err:
        print('Could not download the CJK font; Chinese text in figures may not display:', err)
plt.rcParams['font.family'] = [f for f in _cjk if f in _have] + ['DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def print_header(title):
    print("\n" + "=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

# ============================================================
print_header("Part 1 | Loading the Olivetti Faces dataset")
# ============================================================

# --- Step 1 ---
print_step(1, "Load the dataset")
data   = fetch_olivetti_faces(shuffle=True, random_state=42)
X      = data.data      # shape: (400, 4096)  — data matrix
y      = data.target    # shape: (400,)        — person labels 0–39
imgs   = data.images    # shape: (400, 64, 64) — image form (for display)

print(f"  Shape of the data matrix X: {X.shape}")
print(f"  Interpretation: m = {X.shape[0]} faces, n = {X.shape[1]} pixels (64×64)")
print(f"  Number of people: {len(np.unique(y))}, {X.shape[0]//len(np.unique(y))} photos each")
print(f"  Range of pixel values: [{X.min():.3f}, {X.max():.3f}]")

# --- Step 2 ---
print_step(2, "Compute the mean face and center the data")
mean_vec = X.mean(axis=0)          # mean face vector, shape: (4096,)
X_c      = X - mean_vec            # centered data matrix X̃

print(f"  Mean face ‖f̄‖ = {np.linalg.norm(mean_vec):.3f}")
print(f"  Mean of the centered X̃ (should be close to 0): {X_c.mean():.2e}")

# --- Step 3 ---
print_step(3, "Visualization: mean face + a sample of each person + the effect of centering")
fig = plt.figure(figsize=(14, 5))
gs  = gridspec.GridSpec(2, 9, hspace=0.1, wspace=0.08)

# mean face (two rows high)
ax_mean = fig.add_subplot(gs[:, 0])
ax_mean.imshow(mean_vec.reshape(64, 64), cmap='gray', vmin=0, vmax=1)
ax_mean.set_title('Mean face\n$\\bar{\\mathbf{f}}$', fontsize=9, color=C_V1)
ax_mean.axis('off')

# one photo of each of 8 people (row 0)
person_idx = [np.where(y == p)[0][0] for p in range(8)]
for col, idx in enumerate(person_idx):
    ax = fig.add_subplot(gs[0, col + 1])
    ax.imshow(imgs[idx], cmap='gray', vmin=0, vmax=1)
    ax.set_title(f'Person {y[idx]}', fontsize=7)
    ax.axis('off')

# the corresponding centered images (row 1)
for col, idx in enumerate(person_idx):
    ax = fig.add_subplot(gs[1, col + 1])
    face_c = X_c[idx].reshape(64, 64)
    lim = max(abs(face_c.min()), abs(face_c.max()))
    ax.imshow(face_c, cmap='RdBu', vmin=-lim, vmax=lim)
    ax.set_title('Centered', fontsize=7, color=C_V2)
    ax.axis('off')

fig.suptitle('Olivetti Faces: mean face (left), original faces (top row), after centering (bottom row)',
             fontsize=10, color=C_AXIS, y=1.01)
plt.show()
print("[Note] Centered images: red = brighter than the mean face, blue = darker than the mean face.")
print("       Individual features (eye shape, face shape) stand out more after centering.")

# --- Step 4 ---
print_step(4, "Visualization: all 400 faces")

# display sorted by person number (the order before shuffling)
sorted_idx = np.argsort(y)   # sort by person label

fig, axes = plt.subplots(
    40, 11,                  # 40 people × (1 spacer + 10 photos)
    figsize=(14, 56),
    gridspec_kw={'wspace': 0.04, 'hspace': 0.06}
)

for person in range(40):
    # axes[person, 0]: person label (text shown on an empty axis)
    ax_label = axes[person, 0]
    ax_label.text(0.5, 0.5, f'Person\n{person:02d}',
                  ha='center', va='center',
                  fontsize=7, color=C_V1,
                  transform=ax_label.transAxes)
    ax_label.axis('off')

    # the 10 photos of this person
    person_imgs = sorted_idx[y[sorted_idx] == person]
    for shot in range(10):
        ax = axes[person, shot + 1]
        ax.imshow(imgs[person_imgs[shot]], cmap='gray', vmin=0, vmax=1)
        ax.axis('off')
        if person == 0:
            ax.set_title(f'#{shot}', fontsize=6, pad=2)

# overall title (its second line gives the mean-face norm ‖f̄‖ and the pixel range)
fig.suptitle(
    f'Olivetti Faces overview: 40 people × 10 photos (400 in all, one row per person)\n'
    f'mean face ‖f̄‖ = {np.linalg.norm(mean_vec):.1f}, pixel range [{X.min():.2f}, {X.max():.2f}]',
    fontsize=10, color=C_AXIS, y=0.9995
)
fig.subplots_adjust(top=1.1)
plt.show()

print(f"[Note] Each row shows the 10 photos of one person; you can observe:")
print(f"       • changes in the direction of lighting (differences in brightness between left and right)")
print(f"       • changes in expression (some people smile or differ in glasses)")
print(f"       • slight turns of the head")

# ============================================================
print_header("Part 2 | Extracting eigenfaces with the SVD")
# ============================================================

# --- Step 1 ---
print_step(1, "Apply the SVD to the centered matrix")
U, s, Vt = np.linalg.svd(X_c, full_matrices=False)
# Vt.shape = (400, 4096); Vt[i] = right singular vector i = eigenface i

print(f"  Shape of U  : {U.shape}")
print(f"  Shape of σ  : {s.shape} ({len(s)} singular values in all)")
print(f"  Shape of V^T: {Vt.shape} (each row is a 64×64 eigenface)")
print(f"\n  First 10 singular values: {s[:10].round(1)}")

# --- Step 2 ---
print_step(2, "Energy distribution of the singular values")
total_energy = (s ** 2).sum()
print(f"  {'k':>5}  {'cumulative':>10}  {'energy bar (each block ≈ 5%)'}")
print(f"  {'-'*50}")
for k in [1, 5, 10, 20, 50, 100, 150, 200]:
    energy = (s[:k] ** 2).sum() / total_energy
    bar    = '█' * int(energy * 20)
    print(f"  {k:>5}  {energy:>10.1%}  {bar}")

# --- Step 3 ---
print_step(3, "Visualize the first 10 eigenfaces")
fig, axes = plt.subplots(2, 5, figsize=(13, 5.5),
                          gridspec_kw={'hspace': 0.3, 'wspace': 0.1})

for i in range(10):
    ef      = Vt[i].reshape(64, 64)
    energy  = s[i]**2 / total_energy * 100
    cumul   = (s[:i+1]**2).sum() / total_energy * 100

    # symmetric color scale (eigenfaces have positive and negative values)
    lim = np.percentile(np.abs(ef), 98)
    r, c = divmod(i, 5)
    axes[r, c].imshow(ef, cmap='RdPu', vmin=-lim, vmax=lim)
    axes[r, c].set_title(
        f'$\\mathbf{{v}}_{i}$   $\\sigma_{i}$={s[i]:.1f}\nenergy {energy:.1f}% (cumulative {cumul:.0f}%)',
        fontsize=7.5
    )
    axes[r, c].axis('off')

fig.suptitle('Olivetti Faces: the first 10 eigenfaces (right singular vectors $\mathbf{v}_i$)', 
             fontsize=11, color=C_AXIS)
plt.show()

print("\n[Observation] Eigenfaces do not look like any real face—they are the directions")
print("       in pixel space in which faces are 'most different.' The first few eigenfaces capture")
print("       overall lighting and broad outlines; the later eigenfaces capture finer differences.")

# ============================================================
print_header("Part 3a | Face reconstruction")
# ============================================================

test_idx = 37   # choose a test face (any integer from 0 to 399 will do)

print_step(1, f"Test face: index {test_idx}, person {y[test_idx]}")
test_c = X_c[test_idx]    # centered test face vector
print(f"  ‖x̃‖ = {np.linalg.norm(test_c):.3f}")

print_step(2, "Reconstruction error for different values of k (relative to ‖x̃‖)")
k_list = [1, 2, 5, 10, 20, 50, 100, 150, 200, 300]
print(f"  {'k':>5}  {'rel. error':>10}  {'cumulative':>10}")
print(f"  {'-'*35}")
for k in k_list:
    c_k   = test_c @ Vt[:k].T        # k eigenface coordinates
    recon = c_k @ Vt[:k]             # reconstruction (centered space)
    err   = np.linalg.norm(test_c - recon) / np.linalg.norm(test_c)
    energy = (s[:k]**2).sum() / total_energy
    print(f"  {k:>5}  {err:>10.4f}  {energy:>10.1%}")

print_step(3, "Visualizing the reconstruction")
k_show = [1, 2, 5, 10, 20, 50]
fig, axes = plt.subplots(1, len(k_show) + 2, figsize=(14, 2.8),
                          gridspec_kw={'wspace': 0.12})

axes[0].imshow(imgs[test_idx], cmap='gray', vmin=0, vmax=1)
axes[0].set_title(f'Original\nperson {y[test_idx]}', fontsize=8)
axes[0].axis('off')

axes[1].imshow(mean_vec.reshape(64, 64), cmap='gray', vmin=0, vmax=1)
axes[1].set_title('Mean face\n$\\bar{\\mathbf{f}}$', fontsize=8)
axes[1].axis('off')

for col, k in enumerate(k_show):
    c_k   = test_c @ Vt[:k].T
    recon = (c_k @ Vt[:k] + mean_vec).reshape(64, 64)
    err   = np.linalg.norm(test_c - (c_k @ Vt[:k])) / np.linalg.norm(test_c)
    axes[col + 2].imshow(np.clip(recon, 0, 1), cmap='gray', vmin=0, vmax=1)
    axes[col + 2].set_title(f'$k={k}$\nerr={err:.2f}', fontsize=8)
    axes[col + 2].axis('off')

fig.suptitle(f'Face reconstruction: person {y[test_idx]} (original → mean → superposition of k eigenfaces)',
             fontsize=10, color=C_AXIS)
plt.show()

# ============================================================
print_header("Part 3b | Face recognition (1-NN accuracy)")
# ============================================================

print_step(1, "Precompute the eigenface coordinates of all the faces")
K_MAX = 200
coords_full = X_c @ Vt[:K_MAX].T    # shape: (400, K_MAX)
print(f"  Shape of the coordinate matrix: {coords_full.shape} (each row = the coordinates of one face in the {K_MAX}-dimensional feature space)")

print_step(2, "Compute the 1-NN recognition rate for different values of k")

def recognition_accuracy(coords_k, y):
    """1-nearest-neighbor recognition rate (leave-one-out)"""
    correct = 0
    for i in range(len(y)):
        dists      = np.linalg.norm(coords_k - coords_k[i], axis=1)
        dists[i]   = np.inf          # exclude the face itself
        pred       = y[np.argmin(dists)]
        correct   += int(pred == y[i])
    return correct / len(y)

k_range = [1, 2, 3, 5, 7, 10, 15, 20, 30, 50, 100, 150, 200]
results = []
print(f"  {'k':>5}  {'accuracy':>8}  {'cumulative':>10}")
print(f"  {'-'*35}")
for k in k_range:
    acc    = recognition_accuracy(coords_full[:, :k], y)
    energy = (s[:k]**2).sum() / total_energy
    bar    = '█' * int(acc * 20)
    print(f"  {k:>5}  {acc:>8.1%}  {energy:>10.1%}  {bar}")
    results.append((k, acc, energy))

print_step(3, "Plot the recognition-rate and energy curves")
k_vals, acc_vals, energy_vals = zip(*results)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
for ax in (ax1, ax2):
    ax.set_facecolor(C_BG)
    ax.grid(True, color=C_GRID, linewidth=0.8, alpha=0.7)

ax1.plot(k_vals, [a*100 for a in acc_vals], color=C_V1,
         linewidth=2.2, marker='o', markersize=5, label='1-NN recognition rate')
k_peak = k_vals[acc_vals.index(max(acc_vals))]
ax1.axvline(k_peak, color=C_WARN, linewidth=1.2, linestyle='--',
            label=f'peak k={k_peak}')
ax1.set_xlabel('Number of eigenfaces $k$', fontsize=10)
ax1.set_ylabel('1-NN recognition rate (%)', fontsize=10)
ax1.set_title('Recognition rate vs $k$', fontsize=11)
ax1.legend(fontsize=9)
ax1.set_ylim(0, 105)

ax2.plot(k_vals, [e*100 for e in energy_vals], color=C_V2,
         linewidth=2.2, marker='s', markersize=5, label='Cumulative energy')
ax2.axhline(90, color=C_AUX, linewidth=1.2, linestyle='--', label='90% threshold')
ax2.axhline(99, color=C_T1,  linewidth=1.2, linestyle='--', label='99% threshold')
ax2.set_xlabel('Number of eigenfaces $k$', fontsize=10)
ax2.set_ylabel('Cumulative energy retained (%)', fontsize=10)
ax2.set_title('Cumulative energy vs $k$', fontsize=11)
ax2.legend(fontsize=9)

fig.suptitle('Dependence of the recognition rate and the retained energy on k (Olivetti Faces, 40 people × 10 photos)',
             fontsize=11, color=C_AXIS)
plt.tight_layout()
plt.show()

best_k   = k_vals[acc_vals.index(max(acc_vals))]
best_acc = max(acc_vals)
print(f"\n[Conclusion] Highest recognition rate: {best_acc:.1%}, reached at k = {best_k}.")
print(f"       Cumulative energy at this k = {(s[:best_k]**2).sum()/total_energy:.1%}")
print(f"       Dimension compression ratio = {X.shape[1]} : {best_k} = {X.shape[1]//best_k} : 1")

# ============================================================
# Part 4 | Interactive experiment (ipywidgets + plotly)
# ============================================================

import plotly.graph_objects as go
from plotly.subplots import make_subplots

# precompute the coordinates for 50 values of k (k=1..50) to speed up the interaction
K_INTERACT = 50
coords_interact = X_c @ Vt[:K_INTERACT].T   # (400, 50)

# colors for the 40 people (book palette + tab20 extension)
import matplotlib.cm as cm
cmap20 = cm.get_cmap('tab20', 40)
PERSON_COLORS = [
    f'rgb({int(c[0]*255)},{int(c[1]*255)},{int(c[2]*255)})'
    for c in [cmap20(i) for i in range(40)]
]

# ── Interactive widgets ──────────────────────────────────────
slider_k = widgets.IntSlider(
    value=5, min=1, max=K_INTERACT, step=1,
    description='k (eigenfaces):',
    style={'description_width': '130px'},
    layout=Layout(width='480px'),
    continuous_update=False
)
slider_person = widgets.IntSlider(
    value=0, min=0, max=399, step=1,
    description='Query face index:',
    style={'description_width': '130px'},
    layout=Layout(width='480px'),
    continuous_update=False
)
output_fig = Output()

def update(change):
    k      = slider_k.value
    p_idx  = slider_person.value

    # --- Reconstruction ---
    c_k    = X_c[p_idx] @ Vt[:k].T
    recon  = np.clip(c_k @ Vt[:k] + mean_vec, 0, 1).reshape(64, 64)
    err    = np.linalg.norm(X_c[p_idx] - c_k @ Vt[:k]) / np.linalg.norm(X_c[p_idx])
    energy = (s[:k]**2).sum() / total_energy

    # --- 1-NN recognition ---
    coords_k = coords_interact[:, :k]
    dists    = np.linalg.norm(coords_k - coords_k[p_idx], axis=1)
    dists[p_idx] = np.inf
    nn_idx   = np.argmin(dists)
    correct  = (y[nn_idx] == y[p_idx])

    with output_fig:
        output_fig.clear_output(wait=True)

        fig = make_subplots(
            rows=1, cols=3,
            subplot_titles=(
                f'Original (person {y[p_idx]})',
                f'Reconstruction k={k} (error={err:.3f})',
                f'Feature space (k={k}, energy={energy:.0%})'
            ),
            column_widths=[0.22, 0.22, 0.56],
            horizontal_spacing=0.06
        )

        # original face
        fig.add_trace(
            go.Heatmap(z=imgs[p_idx][::-1], colorscale='gray',
                       zmin=0, zmax=1, showscale=False),
            row=1, col=1
        )
        # reconstructed face
        fig.add_trace(
            go.Heatmap(z=recon[::-1], colorscale='gray',
                       zmin=0, zmax=1, showscale=False),
            row=1, col=2
        )
        # scatter plot (all 40 people)
        for person in range(40):
            mask = y == person
            fig.add_trace(
                go.Scatter(
                    x=coords_interact[mask, 0],
                    y=coords_interact[mask, 1] if k >= 2 else np.zeros(mask.sum()),
                    mode='markers',
                    marker=dict(color=PERSON_COLORS[person], size=5,
                                opacity=0.6, line=dict(width=0.3, color='white')),
                    name=f'Person {person:02d}',
                    showlegend=(person < 10),
                    legendgroup=f'p{person}'
                ),
                row=1, col=3
            )
        # query point ★
        fig.add_trace(
            go.Scatter(
                x=[coords_interact[p_idx, 0]],
                y=[coords_interact[p_idx, 1] if k >= 2 else 0],
                mode='markers+text',
                marker=dict(symbol='star', size=18, color=C_WARN,
                            line=dict(width=1.5, color='black')),
                text=['Query'], textposition='top center',
                name='Query ★', showlegend=True
            ),
            row=1, col=3
        )
        # nearest-neighbor point ◆
        fig.add_trace(
            go.Scatter(
                x=[coords_interact[nn_idx, 0]],
                y=[coords_interact[nn_idx, 1] if k >= 2 else 0],
                mode='markers+text',
                marker=dict(symbol='diamond', size=14,
                            color='#2AD2C9' if correct else C_WARN,
                            line=dict(width=1.5, color='black')),
                text=[f'Nearest\nperson {y[nn_idx]}'],
                textposition='bottom center',
                name=f'Nearest neighbor (person {y[nn_idx]}) {"✓" if correct else "✗"}',
                showlegend=True
            ),
            row=1, col=3
        )

        result_str = f'Recognition {"correct ✓" if correct else "failed ✗"} (nearest neighbor = person {y[nn_idx]})'
        fig.update_layout(
            height=380, width=980,
            plot_bgcolor=C_BG, paper_bgcolor='white',
            title=dict(
                text=f'k={k} | energy={energy:.0%} | {result_str}',
                font=dict(size=12,
                          color='#2AD2C9' if correct else C_WARN)
            ),
            legend=dict(font=dict(size=7), x=1.02, y=1,
                        tracegroupgap=2),
            margin=dict(l=20, r=150, t=65, b=30)
        )
        for col in (1, 2):
            fig.update_xaxes(showticklabels=False, row=1, col=col)
            fig.update_yaxes(showticklabels=False, row=1, col=col)
        fig.update_xaxes(gridcolor=C_GRID, title_text='PC 0', row=1, col=3)
        fig.update_yaxes(gridcolor=C_GRID, title_text='PC 1', row=1, col=3)

        fig.show()

slider_k.observe(update, names='value')
slider_person.observe(update, names='value')
update(None)

display(VBox([
    HBox([slider_k, slider_person]),
    output_fig
]))