12 Dimensionality Reduction — Exercises¶

Test your understanding of scatter matrices, Fisher's LDA, t-SNE mechanics, and the tradeoffs between PCA, LDA, and t-SNE.

Prerequisites. Read theory.md and work through first_principles.ipynb before attempting these.

In [1]:

Copied!





import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)
import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)

Exercise 1 — Hand Calculation: Scatter Matrices for 2 Classes in 2D¶

Given two classes with the following data points:

Class 0: $x_1 = (1, 2)$, $x_2 = (3, 4)$, $x_3 = (2, 3)$

Class 1: $x_4 = (5, 5)$, $x_5 = (7, 7)$, $x_6 = (6, 8)$

Tasks (by hand, then verify with NumPy):

Compute the class means $\mu_0$, $\mu_1$ and the global mean $\mu$.
Compute the within-class scatter matrix $S_W$.
Compute the between-class scatter matrix $S_B$.
Verify that $S_T = S_W + S_B$ where $S_T = \sum_i (x_i - \mu)(x_i - \mu)^\top$.
For the 2-class case, compute the LDA direction $w^\ast = S_W^{-1}(\mu_1 - \mu_0)$ and normalize to unit length.

Hand-computation guide:

$$\mu_0 = \frac{1}{3}\bigl((1,2) + (3,4) + (2,3)\bigr) = (2, 3)$$

$$\mu_1 = \frac{1}{3}\bigl((5,5) + (7,7) + (6,8)\bigr) = (6, 20/3)$$

$$\mu = \frac{1}{6}\sum_{i=1}^{6} x_i = (4, 29/6)$$

Within-class scatter for class 0: $$S_W^{(0)} = \sum_{i=1}^{3}(x_i - \mu_0)(x_i - \mu_0)^\top$$

Centering class 0: $(−1,−1)$, $(1,1)$, $(0,0)$

$$S_W^{(0)} = \begin{pmatrix}1 & 1\\1 & 1\end{pmatrix} + \begin{pmatrix}1 & 1\\1 & 1\end{pmatrix} + \begin{pmatrix}0 & 0\\0 & 0\end{pmatrix} = \begin{pmatrix}2 & 2\\2 & 2\end{pmatrix}$$

Expected results:

Quantity	Value
$\mu_0$	$(2, 3)$
$\mu_1$	$(6, 20/3) \approx (6, 6.667)$
$\mu$	$(4, 29/6) \approx (4, 4.833)$
$\mu_1 - \mu_0$	$(4, 11/3) \approx (4, 3.667)$

In [2]:

Copied!





# Verify your hand calculations
X_c0 = np.array([[1, 2], [3, 4], [2, 3]], dtype=float)
X_c1 = np.array([[5, 5], [7, 7], [6, 8]], dtype=float)
X_all = np.vstack([X_c0, X_c1])
y_all = np.array([0, 0, 0, 1, 1, 1])

# TODO: Step 1 — class means and global mean
# mu0 = ...
# mu1 = ...
# mu = ...
# print(f'mu0 = {mu0}')
# print(f'mu1 = {mu1}')
# print(f'mu  = {mu}')

# TODO: Step 2 — within-class scatter
# S_W = np.zeros((2, 2))
# for c, Xc in [(0, X_c0), (1, X_c1)]:
#     centered = Xc - Xc.mean(axis=0)
#     S_W += centered.T @ centered
# print(f'\nS_W =\n{S_W}')

# TODO: Step 3 — between-class scatter
# S_B = np.zeros((2, 2))
# for mc, nc in [(mu0, 3), (mu1, 3)]:
#     d = (mc - mu).reshape(-1, 1)
#     S_B += nc * (d @ d.T)
# print(f'\nS_B =\n{S_B}')

# TODO: Step 4 — verify S_T = S_W + S_B
# S_T = (X_all - mu).T @ (X_all - mu)
# assert np.allclose(S_W + S_B, S_T, atol=1e-10)
# print('\nS_W + S_B = S_T verified. ✓')

# TODO: Step 5 — LDA direction for 2 classes
# w_star = np.linalg.solve(S_W, mu1 - mu0)
# w_star = w_star / np.linalg.norm(w_star)
# print(f'\nLDA direction w* = {w_star.round(4)}')
# Verify your hand calculations
X_c0 = np.array([[1, 2], [3, 4], [2, 3]], dtype=float)
X_c1 = np.array([[5, 5], [7, 7], [6, 8]], dtype=float)
X_all = np.vstack([X_c0, X_c1])
y_all = np.array([0, 0, 0, 1, 1, 1])

# TODO: Step 1 — class means and global mean
# mu0 = ...
# mu1 = ...
# mu = ...
# print(f'mu0 = {mu0}')
# print(f'mu1 = {mu1}')
# print(f'mu  = {mu}')

# TODO: Step 2 — within-class scatter
# S_W = np.zeros((2, 2))
# for c, Xc in [(0, X_c0), (1, X_c1)]:
#     centered = Xc - Xc.mean(axis=0)
#     S_W += centered.T @ centered
# print(f'\nS_W =\n{S_W}')

# TODO: Step 3 — between-class scatter
# S_B = np.zeros((2, 2))
# for mc, nc in [(mu0, 3), (mu1, 3)]:
#     d = (mc - mu).reshape(-1, 1)
#     S_B += nc * (d @ d.T)
# print(f'\nS_B =\n{S_B}')

# TODO: Step 4 — verify S_T = S_W + S_B
# S_T = (X_all - mu).T @ (X_all - mu)
# assert np.allclose(S_W + S_B, S_T, atol=1e-10)
# print('\nS_W + S_B = S_T verified. ✓')

# TODO: Step 5 — LDA direction for 2 classes
# w_star = np.linalg.solve(S_W, mu1 - mu0)
# w_star = w_star / np.linalg.norm(w_star)
# print(f'\nLDA direction w* = {w_star.round(4)}')

Exercise 2 — Coding: Implement LDA Projection with Deterministic Check¶

Implement a function lda_project(X, y, n_components) that:

Computes within-class scatter $S_W$ and between-class scatter $S_B$.
Solves $S_W^{-1} S_B w = \lambda w$ for the top eigenvectors.
Projects the data onto the top n_components discriminant directions.

Test dataset: 3 classes in 3D, 50 points each (using rng with seed 42).

Deterministic checks:

The projected data should have shape (150, 2) for n_components=2.
The within-class scatter of the projected data should be smaller than that of a random projection (measured by trace of projected $S_W$).
The Fisher criterion $J = \text{tr}(S_W^{-1} S_B)$ on the original data should equal the sum of the top eigenvalues.

Verification values (with SEED = 42):

Quantity	Expected
Projection shape	(150, 2)
Top 2 eigenvalues sum	> 0
LDA separation > random	True

In [3]:

Copied!





# Fixed 3-class 3D dataset
rng_ex2 = np.random.default_rng(42)
n_class = 50
X_ex2 = np.vstack([
    rng_ex2.normal(loc=[0, 0, 0], scale=1.0, size=(n_class, 3)),
    rng_ex2.normal(loc=[4, 0, 0], scale=1.0, size=(n_class, 3)),
    rng_ex2.normal(loc=[2, 4, 0], scale=1.0, size=(n_class, 3)),
])
y_ex2 = np.array([0]*n_class + [1]*n_class + [2]*n_class)
print(f'Dataset: {X_ex2.shape[0]} samples, {X_ex2.shape[1]} features, '
      f'{len(np.unique(y_ex2))} classes')
# Fixed 3-class 3D dataset
rng_ex2 = np.random.default_rng(42)
n_class = 50
X_ex2 = np.vstack([
    rng_ex2.normal(loc=[0, 0, 0], scale=1.0, size=(n_class, 3)),
    rng_ex2.normal(loc=[4, 0, 0], scale=1.0, size=(n_class, 3)),
    rng_ex2.normal(loc=[2, 4, 0], scale=1.0, size=(n_class, 3)),
])
y_ex2 = np.array([0]*n_class + [1]*n_class + [2]*n_class)
print(f'Dataset: {X_ex2.shape[0]} samples, {X_ex2.shape[1]} features, '
      f'{len(np.unique(y_ex2))} classes')

Dataset: 150 samples, 3 features, 3 classes

In [4]:

Copied!





def lda_project(X, y, n_components):
    """Project data onto top LDA discriminant directions.

    Args:
        X: (n, d) data matrix
        y: (n,) class labels
        n_components: number of discriminant directions

    Returns:
        Z: (n, n_components) projected data
        eigenvalues: (n_components,) eigenvalues of S_W^{-1} S_B
        W: (d, n_components) projection matrix
    """
    # TODO: implement
    # 1. Compute class means and global mean
    # 2. Build S_W and S_B
    # 3. Solve eigenvalue problem on S_W^{-1} S_B
    # 4. Select top n_components eigenvectors
    # 5. Project: Z = (X - global_mean) @ W
    pass


# TODO: run and verify
# Z_ex2, eigvals_ex2, W_ex2 = lda_project(X_ex2, y_ex2, n_components=2)

# Check 1: shape
# assert Z_ex2.shape == (150, 2), f'Wrong shape: {Z_ex2.shape}'
# print(f'Projection shape: {Z_ex2.shape} ✓')

# Check 2: LDA separation better than random projection
# rng_check = np.random.default_rng(42)
# W_rand = rng_check.normal(size=(3, 2))
# W_rand, _ = np.linalg.qr(W_rand)  # orthonormalize
# Z_rand = (X_ex2 - X_ex2.mean(axis=0)) @ W_rand
#
# def fisher_ratio(Z, y):
#     """Compute tr(S_W^{-1} S_B) in projected space."""
#     mu = Z.mean(axis=0)
#     S_W = np.zeros((Z.shape[1], Z.shape[1]))
#     S_B = np.zeros_like(S_W)
#     for c in np.unique(y):
#         Zc = Z[y == c]
#         mc = Zc.mean(axis=0)
#         Zc_c = Zc - mc
#         S_W += Zc_c.T @ Zc_c
#         d = (mc - mu).reshape(-1, 1)
#         S_B += len(Zc) * d @ d.T
#     return np.trace(np.linalg.solve(S_W, S_B))
#
# j_lda = fisher_ratio(Z_ex2, y_ex2)
# j_rand = fisher_ratio(Z_rand, y_ex2)
# print(f'Fisher ratio (LDA):    {j_lda:.4f}')
# print(f'Fisher ratio (random): {j_rand:.4f}')
# assert j_lda > j_rand, 'LDA should separate better than random!'
# print('LDA separation > random separation. ✓')

# Check 3: eigenvalue sum
# print(f'Top 2 eigenvalues: {eigvals_ex2.round(4)}')
# assert np.all(eigvals_ex2 > 0), 'Eigenvalues should be positive'
# print('All checks passed. ✓')
def lda_project(X, y, n_components):
    """Project data onto top LDA discriminant directions.

    Args:
        X: (n, d) data matrix
        y: (n,) class labels
        n_components: number of discriminant directions

    Returns:
        Z: (n, n_components) projected data
        eigenvalues: (n_components,) eigenvalues of S_W^{-1} S_B
        W: (d, n_components) projection matrix
    """
    # TODO: implement
    # 1. Compute class means and global mean
    # 2. Build S_W and S_B
    # 3. Solve eigenvalue problem on S_W^{-1} S_B
    # 4. Select top n_components eigenvectors
    # 5. Project: Z = (X - global_mean) @ W
    pass


# TODO: run and verify
# Z_ex2, eigvals_ex2, W_ex2 = lda_project(X_ex2, y_ex2, n_components=2)

# Check 1: shape
# assert Z_ex2.shape == (150, 2), f'Wrong shape: {Z_ex2.shape}'
# print(f'Projection shape: {Z_ex2.shape} ✓')

# Check 2: LDA separation better than random projection
# rng_check = np.random.default_rng(42)
# W_rand = rng_check.normal(size=(3, 2))
# W_rand, _ = np.linalg.qr(W_rand)  # orthonormalize
# Z_rand = (X_ex2 - X_ex2.mean(axis=0)) @ W_rand
#
# def fisher_ratio(Z, y):
#     """Compute tr(S_W^{-1} S_B) in projected space."""
#     mu = Z.mean(axis=0)
#     S_W = np.zeros((Z.shape[1], Z.shape[1]))
#     S_B = np.zeros_like(S_W)
#     for c in np.unique(y):
#         Zc = Z[y == c]
#         mc = Zc.mean(axis=0)
#         Zc_c = Zc - mc
#         S_W += Zc_c.T @ Zc_c
#         d = (mc - mu).reshape(-1, 1)
#         S_B += len(Zc) * d @ d.T
#     return np.trace(np.linalg.solve(S_W, S_B))
#
# j_lda = fisher_ratio(Z_ex2, y_ex2)
# j_rand = fisher_ratio(Z_rand, y_ex2)
# print(f'Fisher ratio (LDA):    {j_lda:.4f}')
# print(f'Fisher ratio (random): {j_rand:.4f}')
# assert j_lda > j_rand, 'LDA should separate better than random!'
# print('LDA separation > random separation. ✓')

# Check 3: eigenvalue sum
# print(f'Top 2 eigenvalues: {eigvals_ex2.round(4)}')
# assert np.all(eigvals_ex2 > 0), 'Eigenvalues should be positive'
# print('All checks passed. ✓')

In [5]:

Copied!





# TODO: scatter plot of LDA projection
# fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# colors = ['steelblue', 'coral', 'seagreen']
#
# for ax, Z, title in [(axes[0], Z_ex2, 'LDA projection'),
#                       (axes[1], Z_rand, 'Random projection')]:
#     for c in range(3):
#         mask = y_ex2 == c
#         ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
#                    color=colors[c], label=f'Class {c}',
#                    edgecolor='k', linewidth=0.3)
#     ax.set_xlabel('Component 1'); ax.set_ylabel('Component 2')
#     ax.set_title(title); ax.legend(fontsize=8)
#
# plt.suptitle('LDA optimally separates classes vs random projection')
# plt.tight_layout()
# plt.show()
# TODO: scatter plot of LDA projection
# fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# colors = ['steelblue', 'coral', 'seagreen']
#
# for ax, Z, title in [(axes[0], Z_ex2, 'LDA projection'),
#                       (axes[1], Z_rand, 'Random projection')]:
#     for c in range(3):
#         mask = y_ex2 == c
#         ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
#                    color=colors[c], label=f'Class {c}',
#                    edgecolor='k', linewidth=0.3)
#     ax.set_xlabel('Component 1'); ax.set_ylabel('Component 2')
#     ax.set_title(title); ax.legend(fontsize=8)
#
# plt.suptitle('LDA optimally separates classes vs random projection')
# plt.tight_layout()
# plt.show()

Exercise 3 — Conceptual: Why Does t-SNE Use Student-t in Low Dimensions?¶

This exercise explores the crowding problem and the design choice behind the Student-t kernel in t-SNE.

Questions:

Volume scaling. Consider a $d$-dimensional hypersphere of radius $r$. Its volume scales as $r^d$. If you have a point at the center and 10 equidistant neighbors at distance $r$ in $d = 50$ dimensions, approximately how many points could fit in the annular shell between $r$ and $2r$? What does this say about the distribution of pairwise distances in high dimensions?
The crowding problem. When mapping from $d = 50$ to $m = 2$, there is far less "room" in 2D. Points that were at moderate distances in high-d must be placed somewhere in 2D. If we use a Gaussian kernel in both high-d and low-d (as in SNE), why does this cause moderately distant points to crush together in the center of the embedding? (Hint: think about the tails of the Gaussian vs Student-t.)
Student-t solution. The Student-t distribution with 1 degree of freedom (Cauchy) has probability density $f(t) \propto (1 + t^2)^{-1}$. Compare its tail behavior to a Gaussian $f(t) \propto e^{-t^2/2}$. Why do the heavy tails of Student-t help resolve the crowding problem? What happens to $q_{ij}$ for moderately distant points compared to a Gaussian kernel?
Practical consequence. Because of the heavy-tailed low-d kernel:
- Are distances between well-separated clusters in a t-SNE plot meaningful?
- Are relative cluster sizes meaningful?
- What is meaningful in a t-SNE plot?

In [6]:

Copied!





# Scratch space: compare Gaussian vs Student-t kernels

# TODO: Plot both kernels to visualize the tail difference
# t = np.linspace(0, 5, 200)
# gaussian = np.exp(-t**2 / 2)
# student_t = 1.0 / (1.0 + t**2)
#
# fig, ax = plt.subplots(figsize=(7, 4))
# ax.plot(t, gaussian, label='Gaussian', color='steelblue', linewidth=2)
# ax.plot(t, student_t, label='Student-t (ν=1)', color='coral', linewidth=2)
# ax.set_xlabel('distance t')
# ax.set_ylabel('kernel value (unnormalized)')
# ax.set_title('Gaussian vs Student-t kernel: tail behavior')
# ax.legend(fontsize=10)
# ax.set_ylim(0, 1.05)
# plt.tight_layout()
# plt.show()
#
# # At distance t=3:
# print(f'Gaussian at t=3:  {np.exp(-9/2):.6f}')
# print(f'Student-t at t=3: {1/(1+9):.6f}')
# print('Student-t assigns ~1000x more probability to distance 3.')
# print('This lets moderately distant points be placed further apart in 2D.')
# Scratch space: compare Gaussian vs Student-t kernels

# TODO: Plot both kernels to visualize the tail difference
# t = np.linspace(0, 5, 200)
# gaussian = np.exp(-t**2 / 2)
# student_t = 1.0 / (1.0 + t**2)
#
# fig, ax = plt.subplots(figsize=(7, 4))
# ax.plot(t, gaussian, label='Gaussian', color='steelblue', linewidth=2)
# ax.plot(t, student_t, label='Student-t (ν=1)', color='coral', linewidth=2)
# ax.set_xlabel('distance t')
# ax.set_ylabel('kernel value (unnormalized)')
# ax.set_title('Gaussian vs Student-t kernel: tail behavior')
# ax.legend(fontsize=10)
# ax.set_ylim(0, 1.05)
# plt.tight_layout()
# plt.show()
#
# # At distance t=3:
# print(f'Gaussian at t=3:  {np.exp(-9/2):.6f}')
# print(f'Student-t at t=3: {1/(1+9):.6f}')
# print('Student-t assigns ~1000x more probability to distance 3.')
# print('This lets moderately distant points be placed further apart in 2D.')

In [7]:

Copied!





# TODO: compute volume ratio for Q1
# For a d-dimensional sphere:
# V(2r) / V(r) = (2r)^d / r^d = 2^d
# Volume of shell [r, 2r] = V(2r) - V(r) = (2^d - 1) * V(r)
#
# d = 50
# ratio = 2**d - 1
# print(f'In d={d} dimensions:')
# print(f'Shell [r, 2r] is {ratio:.2e} times the volume of sphere [0, r]')
# print(f'Almost all volume is in the outer shell.')
# print(f'This means most pairwise distances concentrate near the same value.')
# TODO: compute volume ratio for Q1
# For a d-dimensional sphere:
# V(2r) / V(r) = (2r)^d / r^d = 2^d
# Volume of shell [r, 2r] = V(2r) - V(r) = (2^d - 1) * V(r)
#
# d = 50
# ratio = 2**d - 1
# print(f'In d={d} dimensions:')
# print(f'Shell [r, 2r] is {ratio:.2e} times the volume of sphere [0, r]')
# print(f'Almost all volume is in the outer shell.')
# print(f'This means most pairwise distances concentrate near the same value.')