12 Dimensionality Reduction (LDA + t-SNE) — First Principles¶

Goal. Build LDA and t-SNE from scratch using only NumPy, compare them with PCA, and verify against sklearn. Explore when each method succeeds or fails.

Prerequisites. PCA (topic 10), eigenvalues/SVD, basic probability (KL divergence).

Theory. See theory.md for derivations of Fisher's criterion, scatter matrices, t-SNE objective, and gradient intuition.

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)

plt.rcParams["figure.dpi"] = 90
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)

plt.rcParams["figure.dpi"] = 90

1. Problem Setup — WHY¶

PCA finds directions of maximum variance, ignoring class labels. When classes overlap along high-variance directions but separate along low-variance directions, PCA fails to reveal class structure.

In [2]:

Copied!





# Generate 2D data where PCA direction != class-discriminative direction
n_per_class = 100

# Class 0: elongated along direction (1, 0.3), centered at (-1, -0.5)
# Class 1: elongated along direction (1, 0.3), centered at (1, 0.5)
class0 = rng.normal(size=(n_per_class, 2)) * [3.0, 0.4] + [-1.0, -1.5]
class1 = rng.normal(size=(n_per_class, 2)) * [3.0, 0.4] + [1.0, 1.5]
X_demo = np.vstack([class0, class1])
y_demo = np.array([0] * n_per_class + [1] * n_per_class)

# PCA direction (max variance)
Xc_demo = X_demo - X_demo.mean(axis=0)
_, _, Vt_demo = np.linalg.svd(Xc_demo, full_matrices=False)
pc1 = Vt_demo[0]  # first PC

# LDA direction (class separation) — for 2 classes: S_W^{-1} (mu1 - mu0)
mu0 = class0.mean(axis=0)
mu1 = class1.mean(axis=0)
diff = mu1 - mu0
lda_dir = diff / np.linalg.norm(diff)  # simplified direction

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# Left: scatter
for c, label, color in [(0, 'Class 0', 'steelblue'), (1, 'Class 1', 'coral')]:
    mask = y_demo == c
    axes[0].scatter(X_demo[mask, 0], X_demo[mask, 1], alpha=0.4, s=15,
                    color=color, label=label, edgecolor='k', linewidth=0.2)
axes[0].set_xlabel('x₁'); axes[0].set_ylabel('x₂')
axes[0].set_title('Original 2D data'); axes[0].legend(fontsize=8)
axes[0].set_aspect('equal')

# Middle: PCA projection
proj_pca = X_demo @ pc1
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_demo == c
    axes[1].hist(proj_pca[mask], bins=20, alpha=0.5, color=color, label=f'Class {c}')
axes[1].set_xlabel('PC1 projection'); axes[1].set_ylabel('count')
axes[1].set_title('PCA projection — classes overlap'); axes[1].legend(fontsize=8)

# Right: LDA projection
proj_lda = X_demo @ lda_dir
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_demo == c
    axes[2].hist(proj_lda[mask], bins=20, alpha=0.5, color=color, label=f'Class {c}')
axes[2].set_xlabel('LDA projection'); axes[2].set_ylabel('count')
axes[2].set_title('LDA projection — classes separated'); axes[2].legend(fontsize=8)

plt.suptitle('PCA maximizes variance (ignores labels) vs LDA maximizes class separation')
plt.tight_layout()
plt.show()
# Generate 2D data where PCA direction != class-discriminative direction
n_per_class = 100

# Class 0: elongated along direction (1, 0.3), centered at (-1, -0.5)
# Class 1: elongated along direction (1, 0.3), centered at (1, 0.5)
class0 = rng.normal(size=(n_per_class, 2)) * [3.0, 0.4] + [-1.0, -1.5]
class1 = rng.normal(size=(n_per_class, 2)) * [3.0, 0.4] + [1.0, 1.5]
X_demo = np.vstack([class0, class1])
y_demo = np.array([0] * n_per_class + [1] * n_per_class)

# PCA direction (max variance)
Xc_demo = X_demo - X_demo.mean(axis=0)
_, _, Vt_demo = np.linalg.svd(Xc_demo, full_matrices=False)
pc1 = Vt_demo[0]  # first PC

# LDA direction (class separation) — for 2 classes: S_W^{-1} (mu1 - mu0)
mu0 = class0.mean(axis=0)
mu1 = class1.mean(axis=0)
diff = mu1 - mu0
lda_dir = diff / np.linalg.norm(diff)  # simplified direction

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# Left: scatter
for c, label, color in [(0, 'Class 0', 'steelblue'), (1, 'Class 1', 'coral')]:
    mask = y_demo == c
    axes[0].scatter(X_demo[mask, 0], X_demo[mask, 1], alpha=0.4, s=15,
                    color=color, label=label, edgecolor='k', linewidth=0.2)
axes[0].set_xlabel('x₁'); axes[0].set_ylabel('x₂')
axes[0].set_title('Original 2D data'); axes[0].legend(fontsize=8)
axes[0].set_aspect('equal')

# Middle: PCA projection
proj_pca = X_demo @ pc1
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_demo == c
    axes[1].hist(proj_pca[mask], bins=20, alpha=0.5, color=color, label=f'Class {c}')
axes[1].set_xlabel('PC1 projection'); axes[1].set_ylabel('count')
axes[1].set_title('PCA projection — classes overlap'); axes[1].legend(fontsize=8)

# Right: LDA projection
proj_lda = X_demo @ lda_dir
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_demo == c
    axes[2].hist(proj_lda[mask], bins=20, alpha=0.5, color=color, label=f'Class {c}')
axes[2].set_xlabel('LDA projection'); axes[2].set_ylabel('count')
axes[2].set_title('LDA projection — classes separated'); axes[2].legend(fontsize=8)

plt.suptitle('PCA maximizes variance (ignores labels) vs LDA maximizes class separation')
plt.tight_layout()
plt.show()

No description has been provided for this image

2. Mathematical Core — WHAT¶

LDA: Scatter Matrices¶

Within-class scatter: $S_W = \sum_k \sum_{i: y_i=k} (x_i - \mu_k)(x_i - \mu_k)^\top$
Between-class scatter: $S_B = \sum_k n_k (\mu_k - \mu)(\mu_k - \mu)^\top$
Fisher criterion: $J(w) = \frac{w^\top S_B w}{w^\top S_W w}$
Solution: eigenvectors of $S_W^{-1} S_B$ (at most $C-1$ non-trivial)

t-SNE: Neighborhood Preservation¶

High-d similarity: Gaussian kernel $p_{ij}$ with per-point bandwidth $\sigma_i$
Low-d similarity: Student-t kernel $q_{ij} \propto (1 + \|z_i - z_j\|^2)^{-1}$
Objective: minimize $D_{\text{KL}}(P \| Q)$
Gradient: attraction (neighbors) + repulsion (non-neighbors)

In [3]:

Copied!





# Visualize scatter matrices on the 2D demo data
mu_global = X_demo.mean(axis=0)

# Within-class scatter
S_W_demo = np.zeros((2, 2))
for c in [0, 1]:
    Xc = X_demo[y_demo == c] - X_demo[y_demo == c].mean(axis=0)
    S_W_demo += Xc.T @ Xc

# Between-class scatter
S_B_demo = np.zeros((2, 2))
for c in [0, 1]:
    n_c = np.sum(y_demo == c)
    mu_c = X_demo[y_demo == c].mean(axis=0)
    diff_c = (mu_c - mu_global).reshape(-1, 1)
    S_B_demo += n_c * (diff_c @ diff_c.T)

print('Within-class scatter S_W:')
print(S_W_demo.round(2))
print(f'\nBetween-class scatter S_B:')
print(S_B_demo.round(2))
print(f'\nS_T = S_W + S_B:')
S_T_demo = S_W_demo + S_B_demo
print(S_T_demo.round(2))

# Verify decomposition: S_T should equal total scatter
Xc_total = X_demo - mu_global
S_T_direct = Xc_total.T @ Xc_total
assert np.allclose(S_T_demo, S_T_direct, atol=1e-10), 'S_T decomposition failed'
print('\nS_W + S_B = S_T verified. ✓')
# Visualize scatter matrices on the 2D demo data
mu_global = X_demo.mean(axis=0)

# Within-class scatter
S_W_demo = np.zeros((2, 2))
for c in [0, 1]:
    Xc = X_demo[y_demo == c] - X_demo[y_demo == c].mean(axis=0)
    S_W_demo += Xc.T @ Xc

# Between-class scatter
S_B_demo = np.zeros((2, 2))
for c in [0, 1]:
    n_c = np.sum(y_demo == c)
    mu_c = X_demo[y_demo == c].mean(axis=0)
    diff_c = (mu_c - mu_global).reshape(-1, 1)
    S_B_demo += n_c * (diff_c @ diff_c.T)

print('Within-class scatter S_W:')
print(S_W_demo.round(2))
print(f'\nBetween-class scatter S_B:')
print(S_B_demo.round(2))
print(f'\nS_T = S_W + S_B:')
S_T_demo = S_W_demo + S_B_demo
print(S_T_demo.round(2))

# Verify decomposition: S_T should equal total scatter
Xc_total = X_demo - mu_global
S_T_direct = Xc_total.T @ Xc_total
assert np.allclose(S_T_demo, S_T_direct, atol=1e-10), 'S_T decomposition failed'
print('\nS_W + S_B = S_T verified. ✓')

Within-class scatter S_W:
[[1748.39   12.96]
 [  12.96   26.72]]

Between-class scatter S_B:
[[246.58 334.33]
 [334.33 453.3 ]]

S_T = S_W + S_B:
[[1994.98  347.29]
 [ 347.29  480.02]]

S_W + S_B = S_T verified. ✓

3. Solution Method — HOW¶

LDA Algorithm¶

Compute class means $\mu_k$ and global mean $\mu$.
Build $S_W$ and $S_B$.
Solve the generalized eigenvalue problem $S_W^{-1} S_B w = \lambda w$.
Take the top $k \leq C - 1$ eigenvectors as projection directions.

t-SNE Algorithm¶

Compute pairwise distances in high-d.
For each point, find $\sigma_i$ via binary search to match target perplexity.
Compute symmetrized $p_{ij}$.
Initialize embedding $z_i$ (e.g., small random or PCA).
Gradient descent on $D_{\text{KL}}(P \| Q)$ with momentum.

4. Implementation — BUILD¶

In [4]:

Copied!





class LDAScratch:
    """Fisher's Linear Discriminant Analysis from scratch.

    Attributes after fit:
        scalings_:   (d, n_components) — discriminant directions (columns)
        eigenvalues_: (n_components,) — corresponding eigenvalues of S_W^{-1} S_B
        classes_:    unique class labels
        means_:      (C, d) — class means
        global_mean_: (d,) — overall mean
    """

    def __init__(self, n_components=None):
        self.n_components = n_components

    def fit(self, X, y):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y)
        n, d = X.shape
        self.classes_ = np.unique(y)
        C = len(self.classes_)
        if self.n_components is None:
            self.n_components = min(d, C - 1)

        # Class means and global mean
        self.global_mean_ = X.mean(axis=0)
        self.means_ = np.array([X[y == c].mean(axis=0) for c in self.classes_])

        # Within-class scatter
        S_W = np.zeros((d, d))
        for c in self.classes_:
            Xc = X[y == c] - X[y == c].mean(axis=0)
            S_W += Xc.T @ Xc

        # Between-class scatter
        S_B = np.zeros((d, d))
        for i, c in enumerate(self.classes_):
            n_c = np.sum(y == c)
            diff = (self.means_[i] - self.global_mean_).reshape(-1, 1)
            S_B += n_c * (diff @ diff.T)

        # Solve generalized eigenvalue problem: S_W^{-1} S_B w = lambda w
        # Use scipy-free approach: np.linalg.eig on S_W^{-1} S_B
        # S_W is singular when features are collinear or classes have few samples
        # (e.g. constant pixels in digits), so add a small shrinkage ridge.
        eps = 1e-6 * np.trace(S_W) / d
        A = np.linalg.solve(S_W + eps * np.eye(d), S_B)  # = (S_W + eps I)^{-1} S_B
        eigenvalues, eigenvectors = np.linalg.eig(A)

        # Eigenvalues may have small imaginary parts due to numerics; take real
        eigenvalues = np.real(eigenvalues)
        eigenvectors = np.real(eigenvectors)

        # Sort by descending eigenvalue
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]

        k = self.n_components
        self.eigenvalues_ = eigenvalues[:k]
        self.scalings_ = eigenvectors[:, :k]
        return self

    def transform(self, X):
        X = np.asarray(X, dtype=float)
        return (X - self.global_mean_) @ self.scalings_

    def fit_transform(self, X, y):
        return self.fit(X, y).transform(X)


# Quick test on demo data
lda = LDAScratch(n_components=1).fit(X_demo, y_demo)
Z_lda = lda.transform(X_demo)
print(f'LDA eigenvalue: {lda.eigenvalues_[0]:.4f}')
print(f'LDA direction:  {lda.scalings_[:, 0].round(4)}')
print(f'Projection shape: {Z_lda.shape}')
class LDAScratch:
    """Fisher's Linear Discriminant Analysis from scratch.

    Attributes after fit:
        scalings_:   (d, n_components) — discriminant directions (columns)
        eigenvalues_: (n_components,) — corresponding eigenvalues of S_W^{-1} S_B
        classes_:    unique class labels
        means_:      (C, d) — class means
        global_mean_: (d,) — overall mean
    """

    def __init__(self, n_components=None):
        self.n_components = n_components

    def fit(self, X, y):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y)
        n, d = X.shape
        self.classes_ = np.unique(y)
        C = len(self.classes_)
        if self.n_components is None:
            self.n_components = min(d, C - 1)

        # Class means and global mean
        self.global_mean_ = X.mean(axis=0)
        self.means_ = np.array([X[y == c].mean(axis=0) for c in self.classes_])

        # Within-class scatter
        S_W = np.zeros((d, d))
        for c in self.classes_:
            Xc = X[y == c] - X[y == c].mean(axis=0)
            S_W += Xc.T @ Xc

        # Between-class scatter
        S_B = np.zeros((d, d))
        for i, c in enumerate(self.classes_):
            n_c = np.sum(y == c)
            diff = (self.means_[i] - self.global_mean_).reshape(-1, 1)
            S_B += n_c * (diff @ diff.T)

        # Solve generalized eigenvalue problem: S_W^{-1} S_B w = lambda w
        # Use scipy-free approach: np.linalg.eig on S_W^{-1} S_B
        # S_W is singular when features are collinear or classes have few samples
        # (e.g. constant pixels in digits), so add a small shrinkage ridge.
        eps = 1e-6 * np.trace(S_W) / d
        A = np.linalg.solve(S_W + eps * np.eye(d), S_B)  # = (S_W + eps I)^{-1} S_B
        eigenvalues, eigenvectors = np.linalg.eig(A)

        # Eigenvalues may have small imaginary parts due to numerics; take real
        eigenvalues = np.real(eigenvalues)
        eigenvectors = np.real(eigenvectors)

        # Sort by descending eigenvalue
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]

        k = self.n_components
        self.eigenvalues_ = eigenvalues[:k]
        self.scalings_ = eigenvectors[:, :k]
        return self

    def transform(self, X):
        X = np.asarray(X, dtype=float)
        return (X - self.global_mean_) @ self.scalings_

    def fit_transform(self, X, y):
        return self.fit(X, y).transform(X)


# Quick test on demo data
lda = LDAScratch(n_components=1).fit(X_demo, y_demo)
Z_lda = lda.transform(X_demo)
print(f'LDA eigenvalue: {lda.eigenvalues_[0]:.4f}')
print(f'LDA direction:  {lda.scalings_[:, 0].round(4)}')
print(f'Projection shape: {Z_lda.shape}')

LDA eigenvalue: 16.9803
LDA direction:  [-0.0039 -1.    ]
Projection shape: (200, 1)

In [5]:

Copied!





class TSNEScratch:
    """t-SNE from scratch using gradient descent.

    Simplified implementation for educational purposes.
    Uses exact (not Barnes-Hut) computation.
    """

    def __init__(self, n_components=2, perplexity=30.0, n_iter=500,
                 learning_rate=200.0, momentum=0.8, seed=42):
        self.n_components = n_components
        self.perplexity = perplexity
        self.n_iter = n_iter
        self.learning_rate = learning_rate
        self.momentum = momentum
        self.seed = seed

    def _pairwise_distances(self, X):
        """Compute squared Euclidean distance matrix."""
        sum_sq = np.sum(X ** 2, axis=1)
        D = sum_sq[:, None] + sum_sq[None, :] - 2.0 * X @ X.T
        np.fill_diagonal(D, 0.0)
        return np.maximum(D, 0.0)

    def _compute_perplexity_and_p(self, D_sq, target_perplexity):
        """Compute P matrix with per-point bandwidth via binary search."""
        n = D_sq.shape[0]
        P = np.zeros((n, n))
        target_entropy = np.log(target_perplexity)

        for i in range(n):
            # Binary search for sigma_i
            lo, hi = 1e-10, 1e4
            for _ in range(50):  # binary search iterations
                sigma = (lo + hi) / 2.0
                # Conditional probabilities p_{j|i}
                dists = D_sq[i].copy()
                dists[i] = np.inf  # exclude self
                logits = -dists / (2.0 * sigma ** 2)
                logits -= np.max(logits)  # numerical stability
                exp_logits = np.exp(logits)
                exp_logits[i] = 0.0
                sum_exp = np.sum(exp_logits)
                if sum_exp == 0:
                    lo = sigma
                    continue
                p_cond = exp_logits / sum_exp

                # Entropy
                p_safe = np.maximum(p_cond, 1e-12)
                entropy = -np.sum(p_cond * np.log(p_safe))

                if entropy > target_entropy:
                    hi = sigma
                else:
                    lo = sigma

            P[i] = p_cond

        # Symmetrize
        P = (P + P.T) / (2.0 * n)
        P = np.maximum(P, 1e-12)
        return P

    def _compute_q(self, Y):
        """Compute Q matrix using Student-t kernel."""
        D_sq = self._pairwise_distances(Y)
        inv = 1.0 / (1.0 + D_sq)
        np.fill_diagonal(inv, 0.0)
        Q = inv / np.sum(inv)
        Q = np.maximum(Q, 1e-12)
        return Q, inv

    def _kl_divergence(self, P, Q):
        """KL(P || Q)."""
        return np.sum(P * np.log(P / Q))

    def fit_transform(self, X):
        X = np.asarray(X, dtype=float)
        n = X.shape[0]
        rng_tsne = np.random.default_rng(self.seed)

        # Step 1: Compute pairwise distances and P matrix
        D_sq = self._pairwise_distances(X)
        P = self._compute_perplexity_and_p(D_sq, self.perplexity)

        # Early exaggeration: multiply P by 4 for the first 100 iterations
        P_exag = P * 4.0

        # Step 2: Initialize embedding
        Y = rng_tsne.normal(scale=1e-4, size=(n, self.n_components))
        velocity = np.zeros_like(Y)

        self.kl_history_ = []

        for it in range(self.n_iter):
            # Use exaggerated P for first 100 iterations
            P_use = P_exag if it < 100 else P

            Q, inv_dist = self._compute_q(Y)
            self.kl_history_.append(self._kl_divergence(P_use, Q))

            # Gradient: dC/dY_i = 4 * sum_j (p_ij - q_ij) * (y_i - y_j) * (1 + ||y_i-y_j||^2)^{-1}
            PQ_diff = P_use - Q
            grad = np.zeros_like(Y)
            for i in range(n):
                diff = Y[i] - Y  # (n, m)
                grad[i] = 4.0 * np.sum(
                    (PQ_diff[i, :, None]) * diff * inv_dist[i, :, None],
                    axis=0
                )

            # Update with momentum
            velocity = self.momentum * velocity - self.learning_rate * grad
            Y += velocity

            # Center embedding
            Y -= Y.mean(axis=0)

        self.embedding_ = Y
        return Y


print('TSNEScratch class defined.')
class TSNEScratch:
    """t-SNE from scratch using gradient descent.

    Simplified implementation for educational purposes.
    Uses exact (not Barnes-Hut) computation.
    """

    def __init__(self, n_components=2, perplexity=30.0, n_iter=500,
                 learning_rate=200.0, momentum=0.8, seed=42):
        self.n_components = n_components
        self.perplexity = perplexity
        self.n_iter = n_iter
        self.learning_rate = learning_rate
        self.momentum = momentum
        self.seed = seed

    def _pairwise_distances(self, X):
        """Compute squared Euclidean distance matrix."""
        sum_sq = np.sum(X ** 2, axis=1)
        D = sum_sq[:, None] + sum_sq[None, :] - 2.0 * X @ X.T
        np.fill_diagonal(D, 0.0)
        return np.maximum(D, 0.0)

    def _compute_perplexity_and_p(self, D_sq, target_perplexity):
        """Compute P matrix with per-point bandwidth via binary search."""
        n = D_sq.shape[0]
        P = np.zeros((n, n))
        target_entropy = np.log(target_perplexity)

        for i in range(n):
            # Binary search for sigma_i
            lo, hi = 1e-10, 1e4
            for _ in range(50):  # binary search iterations
                sigma = (lo + hi) / 2.0
                # Conditional probabilities p_{j|i}
                dists = D_sq[i].copy()
                dists[i] = np.inf  # exclude self
                logits = -dists / (2.0 * sigma ** 2)
                logits -= np.max(logits)  # numerical stability
                exp_logits = np.exp(logits)
                exp_logits[i] = 0.0
                sum_exp = np.sum(exp_logits)
                if sum_exp == 0:
                    lo = sigma
                    continue
                p_cond = exp_logits / sum_exp

                # Entropy
                p_safe = np.maximum(p_cond, 1e-12)
                entropy = -np.sum(p_cond * np.log(p_safe))

                if entropy > target_entropy:
                    hi = sigma
                else:
                    lo = sigma

            P[i] = p_cond

        # Symmetrize
        P = (P + P.T) / (2.0 * n)
        P = np.maximum(P, 1e-12)
        return P

    def _compute_q(self, Y):
        """Compute Q matrix using Student-t kernel."""
        D_sq = self._pairwise_distances(Y)
        inv = 1.0 / (1.0 + D_sq)
        np.fill_diagonal(inv, 0.0)
        Q = inv / np.sum(inv)
        Q = np.maximum(Q, 1e-12)
        return Q, inv

    def _kl_divergence(self, P, Q):
        """KL(P || Q)."""
        return np.sum(P * np.log(P / Q))

    def fit_transform(self, X):
        X = np.asarray(X, dtype=float)
        n = X.shape[0]
        rng_tsne = np.random.default_rng(self.seed)

        # Step 1: Compute pairwise distances and P matrix
        D_sq = self._pairwise_distances(X)
        P = self._compute_perplexity_and_p(D_sq, self.perplexity)

        # Early exaggeration: multiply P by 4 for the first 100 iterations
        P_exag = P * 4.0

        # Step 2: Initialize embedding
        Y = rng_tsne.normal(scale=1e-4, size=(n, self.n_components))
        velocity = np.zeros_like(Y)

        self.kl_history_ = []

        for it in range(self.n_iter):
            # Use exaggerated P for first 100 iterations
            P_use = P_exag if it < 100 else P

            Q, inv_dist = self._compute_q(Y)
            self.kl_history_.append(self._kl_divergence(P_use, Q))

            # Gradient: dC/dY_i = 4 * sum_j (p_ij - q_ij) * (y_i - y_j) * (1 + ||y_i-y_j||^2)^{-1}
            PQ_diff = P_use - Q
            grad = np.zeros_like(Y)
            for i in range(n):
                diff = Y[i] - Y  # (n, m)
                grad[i] = 4.0 * np.sum(
                    (PQ_diff[i, :, None]) * diff * inv_dist[i, :, None],
                    axis=0
                )

            # Update with momentum
            velocity = self.momentum * velocity - self.learning_rate * grad
            Y += velocity

            # Center embedding
            Y -= Y.mean(axis=0)

        self.embedding_ = Y
        return Y


print('TSNEScratch class defined.')

TSNEScratch class defined.

In [6]:

Copied!





# Load Iris dataset for comparison experiments
try:
    from sklearn.datasets import load_iris, load_digits
    iris = load_iris()
    X_iris, y_iris = iris.data, iris.target
    print(f'Iris: {X_iris.shape[0]} samples, {X_iris.shape[1]} features, '
          f'{len(np.unique(y_iris))} classes')

    # Subsample digits for speed (t-SNE is O(n^2))
    digits = load_digits()
    X_digits_full, y_digits_full = digits.data, digits.target
    # Take 300 samples for tractable t-SNE
    idx_sub = rng.choice(len(X_digits_full), size=300, replace=False)
    X_digits = X_digits_full[idx_sub]
    y_digits = y_digits_full[idx_sub]
    print(f'Digits (subsampled): {X_digits.shape[0]} samples, '
          f'{X_digits.shape[1]} features, {len(np.unique(y_digits))} classes')
    HAS_SKLEARN_DATA = True
except ModuleNotFoundError:
    print('scikit-learn not installed; using synthetic data.')
    HAS_SKLEARN_DATA = False
    # Synthetic 3-class data in 4D
    n_synth = 50
    X_iris = np.vstack([
        rng.normal(loc=[0, 0, 0, 0], scale=0.5, size=(n_synth, 4)),
        rng.normal(loc=[3, 3, 0, 0], scale=0.5, size=(n_synth, 4)),
        rng.normal(loc=[0, 3, 3, 0], scale=0.5, size=(n_synth, 4)),
    ])
    y_iris = np.array([0]*n_synth + [1]*n_synth + [2]*n_synth)
    X_digits = X_iris.copy()
    y_digits = y_iris.copy()
# Load Iris dataset for comparison experiments
try:
    from sklearn.datasets import load_iris, load_digits
    iris = load_iris()
    X_iris, y_iris = iris.data, iris.target
    print(f'Iris: {X_iris.shape[0]} samples, {X_iris.shape[1]} features, '
          f'{len(np.unique(y_iris))} classes')

    # Subsample digits for speed (t-SNE is O(n^2))
    digits = load_digits()
    X_digits_full, y_digits_full = digits.data, digits.target
    # Take 300 samples for tractable t-SNE
    idx_sub = rng.choice(len(X_digits_full), size=300, replace=False)
    X_digits = X_digits_full[idx_sub]
    y_digits = y_digits_full[idx_sub]
    print(f'Digits (subsampled): {X_digits.shape[0]} samples, '
          f'{X_digits.shape[1]} features, {len(np.unique(y_digits))} classes')
    HAS_SKLEARN_DATA = True
except ModuleNotFoundError:
    print('scikit-learn not installed; using synthetic data.')
    HAS_SKLEARN_DATA = False
    # Synthetic 3-class data in 4D
    n_synth = 50
    X_iris = np.vstack([
        rng.normal(loc=[0, 0, 0, 0], scale=0.5, size=(n_synth, 4)),
        rng.normal(loc=[3, 3, 0, 0], scale=0.5, size=(n_synth, 4)),
        rng.normal(loc=[0, 3, 3, 0], scale=0.5, size=(n_synth, 4)),
    ])
    y_iris = np.array([0]*n_synth + [1]*n_synth + [2]*n_synth)
    X_digits = X_iris.copy()
    y_digits = y_iris.copy()

Iris: 150 samples, 4 features, 3 classes
Digits (subsampled): 300 samples, 64 features, 10 classes

In [7]:

Copied!





# PCA vs LDA on Iris

# PCA (from scratch, 2 components)
X_iris_c = X_iris - X_iris.mean(axis=0)
_, _, Vt_iris = np.linalg.svd(X_iris_c, full_matrices=False)
Z_pca_iris = X_iris_c @ Vt_iris[:2].T

# LDA (from scratch, 2 components — Iris has 3 classes)
lda_iris = LDAScratch(n_components=2).fit(X_iris, y_iris)
Z_lda_iris = lda_iris.transform(X_iris)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
colors = ['steelblue', 'coral', 'seagreen']
class_names = ['Class 0', 'Class 1', 'Class 2']

for ax, Z, title in [(axes[0], Z_pca_iris, 'PCA (unsupervised)'),
                      (axes[1], Z_lda_iris, 'LDA (supervised)')]:
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
                   color=colors[c], label=class_names[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('Iris dataset: PCA vs LDA projection')
plt.tight_layout()
plt.show()

print(f'LDA eigenvalues: {lda_iris.eigenvalues_.round(4)}')
# PCA vs LDA on Iris

# PCA (from scratch, 2 components)
X_iris_c = X_iris - X_iris.mean(axis=0)
_, _, Vt_iris = np.linalg.svd(X_iris_c, full_matrices=False)
Z_pca_iris = X_iris_c @ Vt_iris[:2].T

# LDA (from scratch, 2 components — Iris has 3 classes)
lda_iris = LDAScratch(n_components=2).fit(X_iris, y_iris)
Z_lda_iris = lda_iris.transform(X_iris)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
colors = ['steelblue', 'coral', 'seagreen']
class_names = ['Class 0', 'Class 1', 'Class 2']

for ax, Z, title in [(axes[0], Z_pca_iris, 'PCA (unsupervised)'),
                      (axes[1], Z_lda_iris, 'LDA (supervised)')]:
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
                   color=colors[c], label=class_names[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('Iris dataset: PCA vs LDA projection')
plt.tight_layout()
plt.show()

print(f'LDA eigenvalues: {lda_iris.eigenvalues_.round(4)}')

LDA eigenvalues: [32.1919  0.2854]

In [8]:

Copied!





# t-SNE on Iris
tsne_iris = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500,
                        learning_rate=200.0, seed=42)
Z_tsne_iris = tsne_iris.fit_transform(X_iris)

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))
for ax, Z, title in [(axes[0], Z_pca_iris, 'PCA'),
                      (axes[1], Z_lda_iris, 'LDA'),
                      (axes[2], Z_tsne_iris, 't-SNE')]:
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
                   color=colors[c], label=class_names[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('Iris: PCA vs LDA vs t-SNE')
plt.tight_layout()
plt.show()

print(f'Final t-SNE KL divergence: {tsne_iris.kl_history_[-1]:.4f}')
# t-SNE on Iris
tsne_iris = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500,
                        learning_rate=200.0, seed=42)
Z_tsne_iris = tsne_iris.fit_transform(X_iris)

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))
for ax, Z, title in [(axes[0], Z_pca_iris, 'PCA'),
                      (axes[1], Z_lda_iris, 'LDA'),
                      (axes[2], Z_tsne_iris, 't-SNE')]:
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
                   color=colors[c], label=class_names[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('Iris: PCA vs LDA vs t-SNE')
plt.tight_layout()
plt.show()

print(f'Final t-SNE KL divergence: {tsne_iris.kl_history_[-1]:.4f}')

Final t-SNE KL divergence: 0.1267

In [9]:

Copied!





# KL divergence convergence plot
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(tsne_iris.kl_history_, color='steelblue', linewidth=1)
ax.axvline(100, color='crimson', ls=':', label='End early exaggeration')
ax.set_xlabel('Iteration')
ax.set_ylabel('KL divergence')
ax.set_title('t-SNE convergence on Iris')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
# KL divergence convergence plot
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(tsne_iris.kl_history_, color='steelblue', linewidth=1)
ax.axvline(100, color='crimson', ls=':', label='End early exaggeration')
ax.set_xlabel('Iteration')
ax.set_ylabel('KL divergence')
ax.set_title('t-SNE convergence on Iris')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

In [10]:

Copied!





# PCA vs LDA vs t-SNE on Digits (subsampled)

# PCA
X_dig_c = X_digits - X_digits.mean(axis=0)
_, _, Vt_dig = np.linalg.svd(X_dig_c, full_matrices=False)
Z_pca_dig = X_dig_c @ Vt_dig[:2].T

# LDA (10 classes → max 9 components, take 2)
lda_dig = LDAScratch(n_components=2).fit(X_digits, y_digits)
Z_lda_dig = lda_dig.transform(X_digits)

# t-SNE
tsne_dig = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500,
                       learning_rate=200.0, seed=42)
Z_tsne_dig = tsne_dig.fit_transform(X_digits)

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))
cmap = plt.cm.tab10
for ax, Z, title in [(axes[0], Z_pca_dig, 'PCA'),
                      (axes[1], Z_lda_dig, 'LDA'),
                      (axes[2], Z_tsne_dig, 't-SNE')]:
    for c in range(10):
        mask = y_digits == c
        ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=15,
                   color=cmap(c / 10), label=str(c),
                   edgecolor='k', linewidth=0.2)
    ax.set_xlabel('Component 1'); ax.set_ylabel('Component 2')
    ax.set_title(title)

axes[2].legend(fontsize=6, ncol=2, loc='best', title='digit')
plt.suptitle('Digits dataset: PCA vs LDA vs t-SNE')
plt.tight_layout()
plt.show()
# PCA vs LDA vs t-SNE on Digits (subsampled)

# PCA
X_dig_c = X_digits - X_digits.mean(axis=0)
_, _, Vt_dig = np.linalg.svd(X_dig_c, full_matrices=False)
Z_pca_dig = X_dig_c @ Vt_dig[:2].T

# LDA (10 classes → max 9 components, take 2)
lda_dig = LDAScratch(n_components=2).fit(X_digits, y_digits)
Z_lda_dig = lda_dig.transform(X_digits)

# t-SNE
tsne_dig = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500,
                       learning_rate=200.0, seed=42)
Z_tsne_dig = tsne_dig.fit_transform(X_digits)

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))
cmap = plt.cm.tab10
for ax, Z, title in [(axes[0], Z_pca_dig, 'PCA'),
                      (axes[1], Z_lda_dig, 'LDA'),
                      (axes[2], Z_tsne_dig, 't-SNE')]:
    for c in range(10):
        mask = y_digits == c
        ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=15,
                   color=cmap(c / 10), label=str(c),
                   edgecolor='k', linewidth=0.2)
    ax.set_xlabel('Component 1'); ax.set_ylabel('Component 2')
    ax.set_title(title)

axes[2].legend(fontsize=6, ncol=2, loc='best', title='digit')
plt.suptitle('Digits dataset: PCA vs LDA vs t-SNE')
plt.tight_layout()
plt.show()

5. Library Comparison¶

In [11]:

Copied!





try:
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as SkLDA
except ModuleNotFoundError:
    print('scikit-learn not installed; skipping LDA comparison.')
else:
    # Compare LDA projections on Iris
    sk_lda = SkLDA(n_components=2).fit(X_iris, y_iris)
    Z_sk_lda = sk_lda.transform(X_iris)
    Z_our_lda = lda_iris.transform(X_iris)

    # Directions may differ by sign; check alignment per component
    for j in range(2):
        # Correlate projections (sign-invariant)
        corr = abs(np.corrcoef(Z_our_lda[:, j], Z_sk_lda[:, j])[0, 1])
        print(f'LDA component {j+1}: |correlation| with sklearn = {corr:.6f}')
        assert corr > 0.99, f'LDA component {j+1} does not match sklearn'

    print('LDA matches sklearn. ✓')
try:
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as SkLDA
except ModuleNotFoundError:
    print('scikit-learn not installed; skipping LDA comparison.')
else:
    # Compare LDA projections on Iris
    sk_lda = SkLDA(n_components=2).fit(X_iris, y_iris)
    Z_sk_lda = sk_lda.transform(X_iris)
    Z_our_lda = lda_iris.transform(X_iris)

    # Directions may differ by sign; check alignment per component
    for j in range(2):
        # Correlate projections (sign-invariant)
        corr = abs(np.corrcoef(Z_our_lda[:, j], Z_sk_lda[:, j])[0, 1])
        print(f'LDA component {j+1}: |correlation| with sklearn = {corr:.6f}')
        assert corr > 0.99, f'LDA component {j+1} does not match sklearn'

    print('LDA matches sklearn. ✓')

LDA component 1: |correlation| with sklearn = 1.000000
LDA component 2: |correlation| with sklearn = 1.000000
LDA matches sklearn. ✓

In [12]:

Copied!





try:
    from sklearn.manifold import TSNE as SkTSNE
except ModuleNotFoundError:
    print('scikit-learn not installed; skipping t-SNE comparison.')
else:
    # Both are stochastic, so we compare qualitatively.
    # Run sklearn t-SNE on Iris (sklearn >= 1.5 renamed n_iter to max_iter)
    sk_tsne = SkTSNE(n_components=2, perplexity=30.0, random_state=42,
                     max_iter=500, learning_rate='auto')
    Z_sk_tsne = sk_tsne.fit_transform(X_iris)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    for ax, Z, title in [(axes[0], Z_tsne_iris, 't-SNE (ours)'),
                          (axes[1], Z_sk_tsne, 't-SNE (sklearn)')]:
        for c in range(3):
            mask = y_iris == c
            ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
                       color=colors[c], label=class_names[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('t-SNE comparison: scratch vs sklearn (both preserve local structure)')
    plt.tight_layout()
    plt.show()

    # Quantitative check: both should separate the classes
    # Compute within-class vs between-class distance ratio in embeddings
    def embedding_separation(Z, y):
        """Ratio of mean between-class distance to mean within-class distance."""
        classes = np.unique(y)
        within = []
        between = []
        for c in classes:
            Zc = Z[y == c]
            # Within-class: mean pairwise distance
            if len(Zc) > 1:
                dists = np.sqrt(np.sum((Zc[:, None] - Zc[None, :]) ** 2, axis=2))
                within.append(np.mean(dists[np.triu_indices(len(Zc), k=1)]))
        for i, c1 in enumerate(classes):
            for c2 in classes[i+1:]:
                Z1, Z2 = Z[y == c1], Z[y == c2]
                dists = np.sqrt(np.sum((Z1[:, None] - Z2[None, :]) ** 2, axis=2))
                between.append(np.mean(dists))
        return np.mean(between) / np.mean(within)

    sep_ours = embedding_separation(Z_tsne_iris, y_iris)
    sep_sk = embedding_separation(Z_sk_tsne, y_iris)
    print(f'Separation ratio (ours):    {sep_ours:.2f}')
    print(f'Separation ratio (sklearn): {sep_sk:.2f}')
    print('Both embeddings separate classes (ratio > 1). ✓')
try:
    from sklearn.manifold import TSNE as SkTSNE
except ModuleNotFoundError:
    print('scikit-learn not installed; skipping t-SNE comparison.')
else:
    # Both are stochastic, so we compare qualitatively.
    # Run sklearn t-SNE on Iris (sklearn >= 1.5 renamed n_iter to max_iter)
    sk_tsne = SkTSNE(n_components=2, perplexity=30.0, random_state=42,
                     max_iter=500, learning_rate='auto')
    Z_sk_tsne = sk_tsne.fit_transform(X_iris)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    for ax, Z, title in [(axes[0], Z_tsne_iris, 't-SNE (ours)'),
                          (axes[1], Z_sk_tsne, 't-SNE (sklearn)')]:
        for c in range(3):
            mask = y_iris == c
            ax.scatter(Z[mask, 0], Z[mask, 1], alpha=0.6, s=25,
                       color=colors[c], label=class_names[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('t-SNE comparison: scratch vs sklearn (both preserve local structure)')
    plt.tight_layout()
    plt.show()

    # Quantitative check: both should separate the classes
    # Compute within-class vs between-class distance ratio in embeddings
    def embedding_separation(Z, y):
        """Ratio of mean between-class distance to mean within-class distance."""
        classes = np.unique(y)
        within = []
        between = []
        for c in classes:
            Zc = Z[y == c]
            # Within-class: mean pairwise distance
            if len(Zc) > 1:
                dists = np.sqrt(np.sum((Zc[:, None] - Zc[None, :]) ** 2, axis=2))
                within.append(np.mean(dists[np.triu_indices(len(Zc), k=1)]))
        for i, c1 in enumerate(classes):
            for c2 in classes[i+1:]:
                Z1, Z2 = Z[y == c1], Z[y == c2]
                dists = np.sqrt(np.sum((Z1[:, None] - Z2[None, :]) ** 2, axis=2))
                between.append(np.mean(dists))
        return np.mean(between) / np.mean(within)

    sep_ours = embedding_separation(Z_tsne_iris, y_iris)
    sep_sk = embedding_separation(Z_sk_tsne, y_iris)
    print(f'Separation ratio (ours):    {sep_ours:.2f}')
    print(f'Separation ratio (sklearn): {sep_sk:.2f}')
    print('Both embeddings separate classes (ratio > 1). ✓')

Separation ratio (ours):    7.53
Separation ratio (sklearn): 5.53
Both embeddings separate classes (ratio > 1). ✓

6. Experiments and Failures — VERIFY¶

In [13]:

Copied!





# Verify S_T = S_W + S_B on Iris
mu_iris = X_iris.mean(axis=0)
S_W_iris = np.zeros((X_iris.shape[1], X_iris.shape[1]))
S_B_iris = np.zeros_like(S_W_iris)
for c in np.unique(y_iris):
    Xc = X_iris[y_iris == c]
    mc = Xc.mean(axis=0)
    Xc_centered = Xc - mc
    S_W_iris += Xc_centered.T @ Xc_centered
    diff = (mc - mu_iris).reshape(-1, 1)
    S_B_iris += len(Xc) * diff @ diff.T

S_T_iris = (X_iris - mu_iris).T @ (X_iris - mu_iris)
assert np.allclose(S_W_iris + S_B_iris, S_T_iris, atol=1e-10)
print('S_W + S_B = S_T on Iris dataset. ✓')

# Verify LDA eigenvalue is the Fisher ratio
w = lda_iris.scalings_[:, 0]
J_w = (w @ S_B_iris @ w) / (w @ S_W_iris @ w)
print(f'Fisher criterion J(w) = {J_w:.4f}')
print(f'LDA eigenvalue       = {lda_iris.eigenvalues_[0]:.4f}')
assert np.isclose(J_w, lda_iris.eigenvalues_[0], atol=1e-6)
print('Fisher criterion matches eigenvalue. ✓')
# Verify S_T = S_W + S_B on Iris
mu_iris = X_iris.mean(axis=0)
S_W_iris = np.zeros((X_iris.shape[1], X_iris.shape[1]))
S_B_iris = np.zeros_like(S_W_iris)
for c in np.unique(y_iris):
    Xc = X_iris[y_iris == c]
    mc = Xc.mean(axis=0)
    Xc_centered = Xc - mc
    S_W_iris += Xc_centered.T @ Xc_centered
    diff = (mc - mu_iris).reshape(-1, 1)
    S_B_iris += len(Xc) * diff @ diff.T

S_T_iris = (X_iris - mu_iris).T @ (X_iris - mu_iris)
assert np.allclose(S_W_iris + S_B_iris, S_T_iris, atol=1e-10)
print('S_W + S_B = S_T on Iris dataset. ✓')

# Verify LDA eigenvalue is the Fisher ratio
w = lda_iris.scalings_[:, 0]
J_w = (w @ S_B_iris @ w) / (w @ S_W_iris @ w)
print(f'Fisher criterion J(w) = {J_w:.4f}')
print(f'LDA eigenvalue       = {lda_iris.eigenvalues_[0]:.4f}')
assert np.isclose(J_w, lda_iris.eigenvalues_[0], atol=1e-6)
print('Fisher criterion matches eigenvalue. ✓')

S_W + S_B = S_T on Iris dataset. ✓
Fisher criterion J(w) = 32.1919
LDA eigenvalue       = 32.1919
Fisher criterion matches eigenvalue. ✓

In [14]:

Copied!





# Failure case: t-SNE sensitivity to perplexity
perplexities = [5, 15, 30, 50]
fig, axes = plt.subplots(1, len(perplexities), figsize=(16, 4))

for ax, perp in zip(axes, perplexities):
    tsne_p = TSNEScratch(n_components=2, perplexity=perp, n_iter=500,
                         learning_rate=200.0, seed=42)
    Z_p = tsne_p.fit_transform(X_iris)
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z_p[mask, 0], Z_p[mask, 1], alpha=0.6, s=20,
                   color=colors[c], edgecolor='k', linewidth=0.2)
    ax.set_title(f'perplexity = {perp}')
    ax.set_xlabel('z₁'); ax.set_ylabel('z₂')

plt.suptitle('t-SNE perplexity sensitivity — cluster shape and separation vary')
plt.tight_layout()
plt.show()
print('Low perplexity → fragmented clusters. High perplexity → rounder, more global.')
# Failure case: t-SNE sensitivity to perplexity
perplexities = [5, 15, 30, 50]
fig, axes = plt.subplots(1, len(perplexities), figsize=(16, 4))

for ax, perp in zip(axes, perplexities):
    tsne_p = TSNEScratch(n_components=2, perplexity=perp, n_iter=500,
                         learning_rate=200.0, seed=42)
    Z_p = tsne_p.fit_transform(X_iris)
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z_p[mask, 0], Z_p[mask, 1], alpha=0.6, s=20,
                   color=colors[c], edgecolor='k', linewidth=0.2)
    ax.set_title(f'perplexity = {perp}')
    ax.set_xlabel('z₁'); ax.set_ylabel('z₂')

plt.suptitle('t-SNE perplexity sensitivity — cluster shape and separation vary')
plt.tight_layout()
plt.show()
print('Low perplexity → fragmented clusters. High perplexity → rounder, more global.')

Low perplexity → fragmented clusters. High perplexity → rounder, more global.

In [15]:

Copied!





# Failure case: LDA on non-linearly separable data (concentric circles)
n_circle = 150
theta = rng.uniform(0, 2 * np.pi, n_circle)
r_inner = 1.0 + rng.normal(scale=0.15, size=n_circle)
r_outer = 3.0 + rng.normal(scale=0.3, size=n_circle)

X_inner = np.column_stack([r_inner * np.cos(theta), r_inner * np.sin(theta)])
X_outer = np.column_stack([r_outer * np.cos(theta), r_outer * np.sin(theta)])
X_circles = np.vstack([X_inner, X_outer])
y_circles = np.array([0] * n_circle + [1] * n_circle)

# LDA on circles (linear method → single direction → poor)
lda_circ = LDAScratch(n_components=1).fit(X_circles, y_circles)
Z_lda_circ = lda_circ.transform(X_circles)

# t-SNE on circles
tsne_circ = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500, seed=42)
Z_tsne_circ = tsne_circ.fit_transform(X_circles)

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# Original
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_circles == c
    axes[0].scatter(X_circles[mask, 0], X_circles[mask, 1], alpha=0.5,
                    s=15, color=color, label=f'Class {c}')
axes[0].set_title('Concentric circles (original)')
axes[0].set_aspect('equal'); axes[0].legend(fontsize=8)

# LDA projection (1D → histogram)
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_circles == c
    axes[1].hist(Z_lda_circ[mask, 0], bins=25, alpha=0.5, color=color,
                 label=f'Class {c}')
axes[1].set_title('LDA projection — classes overlap')
axes[1].set_xlabel('LD1'); axes[1].legend(fontsize=8)

# t-SNE embedding
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_circles == c
    axes[2].scatter(Z_tsne_circ[mask, 0], Z_tsne_circ[mask, 1], alpha=0.5,
                    s=15, color=color, label=f'Class {c}')
axes[2].set_title('t-SNE — classes separated'); axes[2].legend(fontsize=8)

plt.suptitle('LDA fails on non-linear structure; t-SNE handles it')
plt.tight_layout()
plt.show()
# Failure case: LDA on non-linearly separable data (concentric circles)
n_circle = 150
theta = rng.uniform(0, 2 * np.pi, n_circle)
r_inner = 1.0 + rng.normal(scale=0.15, size=n_circle)
r_outer = 3.0 + rng.normal(scale=0.3, size=n_circle)

X_inner = np.column_stack([r_inner * np.cos(theta), r_inner * np.sin(theta)])
X_outer = np.column_stack([r_outer * np.cos(theta), r_outer * np.sin(theta)])
X_circles = np.vstack([X_inner, X_outer])
y_circles = np.array([0] * n_circle + [1] * n_circle)

# LDA on circles (linear method → single direction → poor)
lda_circ = LDAScratch(n_components=1).fit(X_circles, y_circles)
Z_lda_circ = lda_circ.transform(X_circles)

# t-SNE on circles
tsne_circ = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500, seed=42)
Z_tsne_circ = tsne_circ.fit_transform(X_circles)

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# Original
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_circles == c
    axes[0].scatter(X_circles[mask, 0], X_circles[mask, 1], alpha=0.5,
                    s=15, color=color, label=f'Class {c}')
axes[0].set_title('Concentric circles (original)')
axes[0].set_aspect('equal'); axes[0].legend(fontsize=8)

# LDA projection (1D → histogram)
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_circles == c
    axes[1].hist(Z_lda_circ[mask, 0], bins=25, alpha=0.5, color=color,
                 label=f'Class {c}')
axes[1].set_title('LDA projection — classes overlap')
axes[1].set_xlabel('LD1'); axes[1].legend(fontsize=8)

# t-SNE embedding
for c, color in [(0, 'steelblue'), (1, 'coral')]:
    mask = y_circles == c
    axes[2].scatter(Z_tsne_circ[mask, 0], Z_tsne_circ[mask, 1], alpha=0.5,
                    s=15, color=color, label=f'Class {c}')
axes[2].set_title('t-SNE — classes separated'); axes[2].legend(fontsize=8)

plt.suptitle('LDA fails on non-linear structure; t-SNE handles it')
plt.tight_layout()
plt.show()

In [16]:

Copied!





# Failure case: t-SNE stochasticity — different seeds give different layouts
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
for ax, seed in zip(axes, [42, 123, 999]):
    tsne_s = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500,
                         learning_rate=200.0, seed=seed)
    Z_s = tsne_s.fit_transform(X_iris)
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z_s[mask, 0], Z_s[mask, 1], alpha=0.6, s=20,
                   color=colors[c], edgecolor='k', linewidth=0.2)
    ax.set_title(f'seed = {seed}')
    ax.set_xlabel('z₁'); ax.set_ylabel('z₂')

plt.suptitle('t-SNE stochasticity — different seeds, different layouts')
plt.tight_layout()
plt.show()
print('Cluster shapes and relative positions change across runs.')
print('Only structures consistent across multiple runs should be trusted.')
# Failure case: t-SNE stochasticity — different seeds give different layouts
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
for ax, seed in zip(axes, [42, 123, 999]):
    tsne_s = TSNEScratch(n_components=2, perplexity=30.0, n_iter=500,
                         learning_rate=200.0, seed=seed)
    Z_s = tsne_s.fit_transform(X_iris)
    for c in range(3):
        mask = y_iris == c
        ax.scatter(Z_s[mask, 0], Z_s[mask, 1], alpha=0.6, s=20,
                   color=colors[c], edgecolor='k', linewidth=0.2)
    ax.set_title(f'seed = {seed}')
    ax.set_xlabel('z₁'); ax.set_ylabel('z₂')

plt.suptitle('t-SNE stochasticity — different seeds, different layouts')
plt.tight_layout()
plt.show()
print('Cluster shapes and relative positions change across runs.')
print('Only structures consistent across multiple runs should be trusted.')

Cluster shapes and relative positions change across runs.
Only structures consistent across multiple runs should be trusted.

7. Connections¶

Theory: theory.md — scatter matrices, Fisher criterion, t-SNE KL objective
PCA: topic 10 — unsupervised linear baseline
Information Theory: foundations — KL divergence, entropy
Autoencoder: topic 17 — nonlinear parametric alternative
Geometry of ML: synthesis

Takeaway¶

PCA: unsupervised, linear, preserves global variance. Baseline.
LDA: supervised, linear, maximizes class separation. Limited to $C-1$ directions.
t-SNE: nonlinear, preserves local neighborhoods. Great for visualization, but stochastic, non-parametric, and global distances are meaningless.
Always compare methods. No single projection suits all tasks.
t-SNE perplexity matters: try multiple values before drawing conclusions.