10 PCA — Exercises¶

Test your understanding of covariance, eigendecomposition, variance maximization, reconstruction, and PCA failure modes.

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

In [1]:

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: Covariance Matrix and First PC¶

Given 3 data points in 2D (already centered):

$$ X_c = \begin{pmatrix} -1 & 2 \\ 0 & -1 \\ 1 & -1 \end{pmatrix} $$

Tasks (by hand, then verify with NumPy):

1. Verify that $X_c$ is mean-centered (column means are zero).
2. Compute the covariance matrix $S = \frac{1}{n} X_c^\top X_c$.
3. Find the eigenvalues $\lambda_1 \ge \lambda_2$ of $S$ by solving $\det(S - \lambda I) = 0$.
4. Find the eigenvector $v_1$ corresponding to $\lambda_1$ (this is the first principal component direction). Normalize to unit length.
5. What fraction of total variance does PC1 capture?

Hand-computation guide:

$$ S = \frac{1}{3} X_c^\top X_c = \frac{1}{3} \begin{pmatrix} (-1)^2 + 0^2 + 1^2 & (-1)(2) + (0)(-1) + (1)(-1) \\ (-1)(2) + (0)(-1) + (1)(-1) & 2^2 + (-1)^2 + (-1)^2 \end{pmatrix} = \frac{1}{3} \begin{pmatrix} 2 & -3 \\ -3 & 6 \end{pmatrix} $$

Characteristic equation: $\lambda^2 - \text{tr}(S)\lambda + \det(S) = 0$

Expected results:

Quantity	Value
$S$	$\begin{pmatrix} 2/3 & -1 \\\\ -1 & 2 \end{pmatrix}$
$\text{tr}(S)$	$8/3 \approx 2.6667$
$\det(S)$	$4/3 - 1 = 1/3 \approx 0.3333$
$\lambda_1$	$\frac{4 + \sqrt{7}}{3} \approx 2.2153$
$\lambda_2$	$\frac{4 - \sqrt{7}}{3} \approx 0.4514$
EVR of PC1	$\lambda_1 / (\lambda_1 + \lambda_2) \approx 83.07\%$

In [2]:

# Verify your hand calculations here
Xc_ex1 = np.array([[-1, 2],
                    [ 0, -1],
                    [ 1, -1]], dtype=float)
n_ex1 = Xc_ex1.shape[0]

# TODO: Step 1 — verify mean-centered
# col_means = ...
# assert np.allclose(col_means, 0, atol=1e-12)

# TODO: Step 2 — compute covariance matrix
# S = ...
# print('S ='); print(S)

# TODO: Step 3 — eigenvalues via characteristic equation
# tr_S = ...
# det_S = ...
# lam1 = ...
# lam2 = ...
# print(f'lambda_1 = {lam1:.4f}, lambda_2 = {lam2:.4f}')

# TODO: Step 4 — eigenvector for lam1, then normalize
# Solve (S - lam1 * I) v = 0
# v1 = ...
# v1 = v1 / np.linalg.norm(v1)

# TODO: Step 5 — explained variance ratio
# evr1 = lam1 / (lam1 + lam2)
# print(f'PC1 captures {evr1:.2%} of total variance')

# Cross-check with np.linalg.eigh
# eigvals, eigvecs = np.linalg.eigh(S)
# print(f'\nnp.linalg.eigh eigenvalues: {eigvals[::-1]}')

# Verify your hand calculations here Xc_ex1 = np.array([[-1, 2], [ 0, -1], [ 1, -1]], dtype=float) n_ex1 = Xc_ex1.shape[0] # TODO: Step 1 — verify mean-centered # col_means = ... # assert np.allclose(col_means, 0, atol=1e-12) # TODO: Step 2 — compute covariance matrix # S = ... # print('S ='); print(S) # TODO: Step 3 — eigenvalues via characteristic equation # tr_S = ... # det_S = ... # lam1 = ... # lam2 = ... # print(f'lambda_1 = {lam1:.4f}, lambda_2 = {lam2:.4f}') # TODO: Step 4 — eigenvector for lam1, then normalize # Solve (S - lam1 * I) v = 0 # v1 = ... # v1 = v1 / np.linalg.norm(v1) # TODO: Step 5 — explained variance ratio # evr1 = lam1 / (lam1 + lam2) # print(f'PC1 captures {evr1:.2%} of total variance') # Cross-check with np.linalg.eigh # eigvals, eigvecs = np.linalg.eigh(S) # print(f'\nnp.linalg.eigh eigenvalues: {eigvals[::-1]}')

Exercise 2 — Coding: Reconstruction Error vs Number of Components¶

The mean squared reconstruction error when keeping $k$ of $d$ components is:

$$ \text{MSE}(k) = \frac{1}{n} \|X_c - X_c W_k W_k^\top\|_F^2 = \sum_{j=k+1}^{d} \lambda_j $$

Tasks:

1. Implement reconstruction_mse(X, k) that fits PCA with $k$ components, reconstructs the data, and returns the MSE.
2. Run it on the dataset below for $k = 1, 2, \ldots, 5$ and verify that:
   • Error decreases monotonically.
   • At $k = d = 5$, error is numerically zero.
   • The errors match the eigenvalue formula: $\text{MSE}(k) = \sum_{j=k+1}^{d} \lambda_j$.
3. Plot reconstruction error vs $k$.

Deterministic check (with SEED = 42):

$k$	MSE (reconstruction)
1	≈ 5.30
2	≈ 2.28
3	≈ 0.89
4	≈ 0.10
5	≈ 0.00

In [3]:

# Fixed dataset for deterministic checks
n_ex2, d_ex2 = 200, 5
rng_ex2 = np.random.default_rng(42)
raw_ex2 = rng_ex2.normal(size=(n_ex2, d_ex2)) * np.array([5.0, 3.0, 2.0, 1.0, 0.5])
R_ex2 = np.linalg.qr(rng_ex2.normal(size=(d_ex2, d_ex2)))[0]
X_ex2 = raw_ex2 @ R_ex2.T

# Fixed dataset for deterministic checks n_ex2, d_ex2 = 200, 5 rng_ex2 = np.random.default_rng(42) raw_ex2 = rng_ex2.normal(size=(n_ex2, d_ex2)) * np.array([5.0, 3.0, 2.0, 1.0, 0.5]) R_ex2 = np.linalg.qr(rng_ex2.normal(size=(d_ex2, d_ex2)))[0] X_ex2 = raw_ex2 @ R_ex2.T

In [4]:

def reconstruction_mse(X, k):
    """Compute mean squared reconstruction error for PCA with k components.

    Steps:
      1. Center X.
      2. Compute SVD (or eigh) to get top-k directions W_k.
      3. Project and reconstruct: X_hat = Xc @ W_k @ W_k.T + mean.
      4. Return mean(sum-of-squares per row).
    """
    # TODO: implement
    pass


# TODO: compute MSE for k = 1, ..., 5
# mse_values = [reconstruction_mse(X_ex2, k) for k in range(1, d_ex2 + 1)]
# for k, mse in enumerate(mse_values, 1):
#     print(f'k={k}  MSE={mse:.4f}')

# TODO: verify monotonic decrease
# for i in range(len(mse_values) - 1):
#     assert mse_values[i] >= mse_values[i + 1] - 1e-10

# TODO: verify MSE ≈ 0 at k = d
# assert mse_values[-1] < 1e-20

# TODO: verify eigenvalue formula
# Xc = X_ex2 - X_ex2.mean(axis=0)
# S = Xc.T @ Xc / n_ex2
# eigvals = np.sort(np.linalg.eigvalsh(S))[::-1]
# for k in range(1, d_ex2 + 1):
#     mse_formula = np.sum(eigvals[k:])
#     assert np.isclose(mse_values[k - 1], mse_formula, atol=1e-8), \
#         f'k={k}: MSE={mse_values[k-1]:.6f} != sum(eigvals[{k}:])={mse_formula:.6f}'
# print('All checks passed. ✓')

def reconstruction_mse(X, k): """Compute mean squared reconstruction error for PCA with k components. Steps: 1. Center X. 2. Compute SVD (or eigh) to get top-k directions W_k. 3. Project and reconstruct: X_hat = Xc @ W_k @ W_k.T + mean. 4. Return mean(sum-of-squares per row). """ # TODO: implement pass # TODO: compute MSE for k = 1, ..., 5 # mse_values = [reconstruction_mse(X_ex2, k) for k in range(1, d_ex2 + 1)] # for k, mse in enumerate(mse_values, 1): # print(f'k={k} MSE={mse:.4f}') # TODO: verify monotonic decrease # for i in range(len(mse_values) - 1): # assert mse_values[i] >= mse_values[i + 1] - 1e-10 # TODO: verify MSE ≈ 0 at k = d # assert mse_values[-1] < 1e-20 # TODO: verify eigenvalue formula # Xc = X_ex2 - X_ex2.mean(axis=0) # S = Xc.T @ Xc / n_ex2 # eigvals = np.sort(np.linalg.eigvalsh(S))[::-1] # for k in range(1, d_ex2 + 1): # mse_formula = np.sum(eigvals[k:]) # assert np.isclose(mse_values[k - 1], mse_formula, atol=1e-8), \ # f'k={k}: MSE={mse_values[k-1]:.6f} != sum(eigvals[{k}:])={mse_formula:.6f}' # print('All checks passed. ✓')

In [5]:

# TODO: bar plot of MSE vs k
# fig, ax = plt.subplots(figsize=(6, 4))
# ax.bar(range(1, d_ex2 + 1), mse_values, color='steelblue', alpha=0.7)
# ax.set_xlabel('number of components k')
# ax.set_ylabel('mean squared reconstruction error')
# ax.set_title('Reconstruction error vs k')
# plt.tight_layout()
# plt.show()

# TODO: bar plot of MSE vs k # fig, ax = plt.subplots(figsize=(6, 4)) # ax.bar(range(1, d_ex2 + 1), mse_values, color='steelblue', alpha=0.7) # ax.set_xlabel('number of components k') # ax.set_ylabel('mean squared reconstruction error') # ax.set_title('Reconstruction error vs k') # plt.tight_layout() # plt.show()

Exercise 3 — Conceptual: When Does PCA Fail?¶

PCA makes specific assumptions. Describe why PCA fails in each scenario and suggest an alternative method:

Questions:

1. Non-linear manifolds. Your data lies on a Swiss roll (a 2D surface curled up in 3D). PCA projects onto a 2D plane. Why does this lose the intrinsic neighborhood structure? What happens to nearby points on the manifold that are far in Euclidean distance? (Suggest: kernel PCA, t-SNE, UMAP, Isomap.)

2. Categorical features. You have a dataset with columns like color ∈ {red, green, blue} encoded as integers {0, 1, 2}. Why is computing a covariance matrix on these values meaningless? What assumption of PCA does this violate? (Suggest: MCA — Multiple Correspondence Analysis, or encode then use a method that respects discrete structure.)

3. Different scales without standardization. Feature $x_1$ is measured in meters (range 0–1) and feature $x_2$ in millimeters (range 0–1000). Without standardization, which feature dominates PC1? Explain using the formula $\text{Var}_v = v^\top S v$ and the diagonal entries of $S$.

4. High noise dimensions. You have 2 signal dimensions with variance 10 each, and 100 noise dimensions with variance 1 each. What fraction of total variance do the signal dimensions carry? Would a 95% explained-variance threshold retain the signal, or mostly noise? Compute the numbers.

In [6]:

# Scratch space for Q3 and Q4 calculations

# Q3: compute covariance for scaled vs unscaled features
# x1_meters = rng.uniform(0, 1, size=100)
# x2_mm = rng.uniform(0, 1000, size=100)
# X_unscaled = np.column_stack([x1_meters, x2_mm])
# S_unscaled = np.cov(X_unscaled, rowvar=False, ddof=0)
# print('S (unscaled):\n', S_unscaled)
# print('PC1 is dominated by x2 because S[1,1] >> S[0,0]')

# Q4: signal vs noise variance
# signal_var = 2 * 10   # 2 dimensions × variance 10
# noise_var = 100 * 1   # 100 dimensions × variance 1
# total_var = signal_var + noise_var
# print(f'Signal fraction: {signal_var / total_var:.1%}')
# print(f'To reach 95% EVR, need to retain noise dims too')

# Scratch space for Q3 and Q4 calculations # Q3: compute covariance for scaled vs unscaled features # x1_meters = rng.uniform(0, 1, size=100) # x2_mm = rng.uniform(0, 1000, size=100) # X_unscaled = np.column_stack([x1_meters, x2_mm]) # S_unscaled = np.cov(X_unscaled, rowvar=False, ddof=0) # print('S (unscaled):\n', S_unscaled) # print('PC1 is dominated by x2 because S[1,1] >> S[0,0]') # Q4: signal vs noise variance # signal_var = 2 * 10 # 2 dimensions × variance 10 # noise_var = 100 * 1 # 100 dimensions × variance 1 # total_var = signal_var + noise_var # print(f'Signal fraction: {signal_var / total_var:.1%}') # print(f'To reach 95% EVR, need to retain noise dims too')

Exercise 4 — Coding: Standardized vs Unstandardized PCA¶

When features have different scales, PCA on raw data vs standardized data gives very different principal directions.

Tasks:

1. Create a 2D dataset where $x_1 \sim \mathcal{N}(0, 1)$ and $x_2 \sim \mathcal{N}(0, 100)$ (10× larger standard deviation), with correlation $\rho = 0.6$ between them.
2. Fit PCA (1 component) on the raw data and on the standardized data (z-scored: subtract mean, divide by std).
3. Compare the PC1 directions. The raw-data PC1 should be nearly aligned with the $x_2$ axis. The standardized PC1 should reflect the correlation structure.
4. Verify that the explained variance ratios differ.

Deterministic check (with SEED = 42):

• Raw PC1: angle to $x_2$ axis $< 5°$
• Standardized PC1: angle to $x_2$ axis $> 30°$
• Raw EVR₁ $>$ Standardized EVR₁

In [7]:

# Generate correlated 2D data with different scales
n_ex4 = 300
rng_ex4 = np.random.default_rng(42)
rho = 0.6
# Covariance matrix: var(x1)=1, var(x2)=10000, cov=0.6*1*100=60
cov_ex4 = np.array([[1.0, rho * 100],
                     [rho * 100, 10000.0]])
X_ex4 = rng_ex4.multivariate_normal([0, 0], cov_ex4, size=n_ex4)


def pca_1_direction(X):
    """Return the first principal component direction (unit vector)."""
    # TODO: center, SVD, return first right singular vector
    pass


def standardize(X):
    """Z-score standardization: (X - mean) / std."""
    # TODO: implement
    pass


# TODO: compute PC1 on raw data
# pc1_raw = pca_1_direction(X_ex4)

# TODO: compute PC1 on standardized data
# X_std = standardize(X_ex4)
# pc1_std = pca_1_direction(X_std)

# Angle to x2 axis (= [0, 1])
# x2_axis = np.array([0.0, 1.0])
# angle_raw = np.degrees(np.arccos(np.clip(abs(pc1_raw @ x2_axis), 0, 1)))
# angle_std = np.degrees(np.arccos(np.clip(abs(pc1_std @ x2_axis), 0, 1)))
# print(f'Raw PC1 angle to x2 axis:          {angle_raw:.1f}°')
# print(f'Standardized PC1 angle to x2 axis:  {angle_std:.1f}°')

# Deterministic checks
# assert angle_raw < 5.0, f'Raw PC1 should be near x2 axis, got {angle_raw:.1f}°'
# assert angle_std > 30.0, f'Standardized PC1 should differ, got {angle_std:.1f}°'
# print('Direction checks passed. ✓')

# Generate correlated 2D data with different scales n_ex4 = 300 rng_ex4 = np.random.default_rng(42) rho = 0.6 # Covariance matrix: var(x1)=1, var(x2)=10000, cov=0.6*1*100=60 cov_ex4 = np.array([[1.0, rho * 100], [rho * 100, 10000.0]]) X_ex4 = rng_ex4.multivariate_normal([0, 0], cov_ex4, size=n_ex4) def pca_1_direction(X): """Return the first principal component direction (unit vector).""" # TODO: center, SVD, return first right singular vector pass def standardize(X): """Z-score standardization: (X - mean) / std.""" # TODO: implement pass # TODO: compute PC1 on raw data # pc1_raw = pca_1_direction(X_ex4) # TODO: compute PC1 on standardized data # X_std = standardize(X_ex4) # pc1_std = pca_1_direction(X_std) # Angle to x2 axis (= [0, 1]) # x2_axis = np.array([0.0, 1.0]) # angle_raw = np.degrees(np.arccos(np.clip(abs(pc1_raw @ x2_axis), 0, 1))) # angle_std = np.degrees(np.arccos(np.clip(abs(pc1_std @ x2_axis), 0, 1))) # print(f'Raw PC1 angle to x2 axis: {angle_raw:.1f}°') # print(f'Standardized PC1 angle to x2 axis: {angle_std:.1f}°') # Deterministic checks # assert angle_raw < 5.0, f'Raw PC1 should be near x2 axis, got {angle_raw:.1f}°' # assert angle_std > 30.0, f'Standardized PC1 should differ, got {angle_std:.1f}°' # print('Direction checks passed. ✓')

In [8]:

# TODO: side-by-side scatter plots showing PC1 direction on raw vs standardized data
# fig, axes = plt.subplots(1, 2, figsize=(12, 5))
#
# # Left: raw data + PC1 arrow
# axes[0].scatter(X_ex4[:, 0], X_ex4[:, 1], alpha=0.3, s=10)
# m = X_ex4.mean(axis=0)
# axes[0].annotate('', xy=m + pc1_raw * 150, xytext=m - pc1_raw * 150,
#                  arrowprops=dict(arrowstyle='->', color='crimson', lw=2))
# axes[0].set_xlabel('x₁ (σ=1)'); axes[0].set_ylabel('x₂ (σ=100)')
# axes[0].set_title('Raw PCA — PC1 aligned with large-scale feature')
#
# # Right: standardized data + PC1 arrow
# axes[1].scatter(X_std[:, 0], X_std[:, 1], alpha=0.3, s=10)
# m_s = X_std.mean(axis=0)
# axes[1].annotate('', xy=m_s + pc1_std * 3, xytext=m_s - pc1_std * 3,
#                  arrowprops=dict(arrowstyle='->', color='crimson', lw=2))
# axes[1].set_xlabel('z₁'); axes[1].set_ylabel('z₂')
# axes[1].set_title('Standardized PCA — PC1 reflects correlation')
# axes[1].set_aspect('equal')
#
# plt.suptitle('Scale dominance: PCA without standardization is misled')
# plt.tight_layout()
# plt.show()

# TODO: side-by-side scatter plots showing PC1 direction on raw vs standardized data # fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # # # Left: raw data + PC1 arrow # axes[0].scatter(X_ex4[:, 0], X_ex4[:, 1], alpha=0.3, s=10) # m = X_ex4.mean(axis=0) # axes[0].annotate('', xy=m + pc1_raw * 150, xytext=m - pc1_raw * 150, # arrowprops=dict(arrowstyle='->', color='crimson', lw=2)) # axes[0].set_xlabel('x₁ (σ=1)'); axes[0].set_ylabel('x₂ (σ=100)') # axes[0].set_title('Raw PCA — PC1 aligned with large-scale feature') # # # Right: standardized data + PC1 arrow # axes[1].scatter(X_std[:, 0], X_std[:, 1], alpha=0.3, s=10) # m_s = X_std.mean(axis=0) # axes[1].annotate('', xy=m_s + pc1_std * 3, xytext=m_s - pc1_std * 3, # arrowprops=dict(arrowstyle='->', color='crimson', lw=2)) # axes[1].set_xlabel('z₁'); axes[1].set_ylabel('z₂') # axes[1].set_title('Standardized PCA — PC1 reflects correlation') # axes[1].set_aspect('equal') # # plt.suptitle('Scale dominance: PCA without standardization is misled') # plt.tight_layout() # plt.show()

Exercise 5 — Conceptual: SVD and Eigendecomposition for PCA¶

The theory shows two routes to PCA:

• Route A (Eigendecomposition): Form $S = \frac{1}{n} X_c^\top X_c$, then eigendecompose: $S = V \Lambda V^\top$.
• Route B (SVD): Compute thin SVD $X_c = U \Sigma V^\top$, then $\lambda_j = \sigma_j^2 / n$ and principal directions are columns of $V$.

Questions:

1. Algebraic equivalence. Starting from $X_c = U \Sigma V^\top$, show that $S = V (\Sigma^2 / n) V^\top$. Why does this prove that the right singular vectors of $X_c$ are the eigenvectors of $S$?

2. Numerical stability. The condition number of $S = X_c^\top X_c$ is the square of the condition number of $X_c$. If $X_c$ has condition number $10^6$, what is the condition number of $S$? Why does this make the eigendecomposition of $S$ less reliable than the SVD of $X_c$ directly?

3. Computational cost. Route A requires forming the $d \times d$ matrix $S$, costing $O(n d^2)$, then eigendecomposing it in $O(d^3)$. Route B computes the thin SVD of the $n \times d$ matrix in $O(\min(n d^2, n^2 d))$. In what regime ($n \gg d$ vs $n \ll d$) is each route more efficient?

4. Score computation. The PCA scores (projected coordinates) are $Z = X_c V_k$. Show that this can also be written as $Z = U_k \Sigma_k$ using only the SVD factors. Why is this useful when $d \gg n$?

In [9]:

# Numerical verification of SVD-eigendecomposition equivalence

# TODO: create a dataset, compute PCA both ways, compare
# n_ex5, d_ex5 = 100, 4
# rng_ex5 = np.random.default_rng(42)
# X_ex5 = rng_ex5.normal(size=(n_ex5, d_ex5)) * [4.0, 2.0, 1.0, 0.3]
# Xc_ex5 = X_ex5 - X_ex5.mean(axis=0)

# Route A: eigendecomposition of covariance
# S_ex5 = Xc_ex5.T @ Xc_ex5 / n_ex5
# eigvals_a, eigvecs_a = np.linalg.eigh(S_ex5)
# idx = eigvals_a.argsort()[::-1]
# eigvals_a = eigvals_a[idx]
# eigvecs_a = eigvecs_a[:, idx]

# Route B: SVD
# U_b, sigma_b, Vt_b = np.linalg.svd(Xc_ex5, full_matrices=False)
# eigvals_b = sigma_b ** 2 / n_ex5
# eigvecs_b = Vt_b.T

# Compare
# assert np.allclose(eigvals_a, eigvals_b, atol=1e-10)
# for j in range(d_ex5):
#     align = abs(float(eigvecs_a[:, j] @ eigvecs_b[:, j]))
#     assert align > 1 - 1e-8, f'PC{j+1} misaligned: {align:.6f}'
# print('Routes A and B produce identical results. ✓')

# Condition number comparison
# cond_Xc = np.linalg.cond(Xc_ex5)
# cond_S = np.linalg.cond(S_ex5)
# print(f'cond(Xc) = {cond_Xc:.2f}')
# print(f'cond(S)  = {cond_S:.2f}')
# print(f'cond(S) ≈ cond(Xc)² = {cond_Xc**2:.2f}')

# Score computation: Z = Xc @ V_k vs Z = U_k @ Sigma_k
# k = 2
# Z_route1 = Xc_ex5 @ eigvecs_b[:, :k]
# Z_route2 = U_b[:, :k] * sigma_b[:k]
# For sign consistency:
# for j in range(k):
#     if Z_route1[0, j] * Z_route2[0, j] < 0:
#         Z_route2[:, j] *= -1
# assert np.allclose(Z_route1, Z_route2, atol=1e-10)
# print(f'Score computation: Xc @ V_k ≡ U_k @ Sigma_k. ✓')

# Numerical verification of SVD-eigendecomposition equivalence # TODO: create a dataset, compute PCA both ways, compare # n_ex5, d_ex5 = 100, 4 # rng_ex5 = np.random.default_rng(42) # X_ex5 = rng_ex5.normal(size=(n_ex5, d_ex5)) * [4.0, 2.0, 1.0, 0.3] # Xc_ex5 = X_ex5 - X_ex5.mean(axis=0) # Route A: eigendecomposition of covariance # S_ex5 = Xc_ex5.T @ Xc_ex5 / n_ex5 # eigvals_a, eigvecs_a = np.linalg.eigh(S_ex5) # idx = eigvals_a.argsort()[::-1] # eigvals_a = eigvals_a[idx] # eigvecs_a = eigvecs_a[:, idx] # Route B: SVD # U_b, sigma_b, Vt_b = np.linalg.svd(Xc_ex5, full_matrices=False) # eigvals_b = sigma_b ** 2 / n_ex5 # eigvecs_b = Vt_b.T # Compare # assert np.allclose(eigvals_a, eigvals_b, atol=1e-10) # for j in range(d_ex5): # align = abs(float(eigvecs_a[:, j] @ eigvecs_b[:, j])) # assert align > 1 - 1e-8, f'PC{j+1} misaligned: {align:.6f}' # print('Routes A and B produce identical results. ✓') # Condition number comparison # cond_Xc = np.linalg.cond(Xc_ex5) # cond_S = np.linalg.cond(S_ex5) # print(f'cond(Xc) = {cond_Xc:.2f}') # print(f'cond(S) = {cond_S:.2f}') # print(f'cond(S) ≈ cond(Xc)² = {cond_Xc**2:.2f}') # Score computation: Z = Xc @ V_k vs Z = U_k @ Sigma_k # k = 2 # Z_route1 = Xc_ex5 @ eigvecs_b[:, :k] # Z_route2 = U_b[:, :k] * sigma_b[:k] # For sign consistency: # for j in range(k): # if Z_route1[0, j] * Z_route2[0, j] < 0: # Z_route2[:, j] *= -1 # assert np.allclose(Z_route1, Z_route2, atol=1e-10) # print(f'Score computation: Xc @ V_k ≡ U_k @ Sigma_k. ✓')