Chapter 11 Experiment 13: Principal Component Analysis of the Iris Dataset (Iris PCA)

Where This Experiment Fits

This experiment is a real-data exercise for the theory of principal component analysis (PCA) in §11.3. Starting from the Eckart–Young theorem (Theorem 8), §11.3 derived the two-step logic of PCA—every projection conserves the total variance (Theorem 10), and then Eckart–Young picks the best projection (Theorem 11). This experiment applies the method to the classic Iris dataset, so that you can see with your own eyes how PCA squeezes 4-dimensional features into a single two-dimensional scatter plot and brings out the cluster structure of the three species.

The experiment also reuses the centering and total-variance tools of Experiment 5 (covariance matrices): the total variance $V_{\text{total}} = \operatorname{tr}(\mathbf{S}) = \frac{1}{m-1}\|\tilde{\mathbf{X}}\|_F^2$ (Exercise 5).

Reused / prerequisite concept	Source	Where it appears in this experiment
Centered data matrix $\tilde{\mathbf{X}}$ , total variance $\operatorname{tr}(\mathbf{S})$	Experiment 5 + §11.3.1	Blocks B, C
The two-step logic of conservation + best projection	Theorem 10, Theorem 11	Block C
Applying the SVD directly to $\tilde{\mathbf{X}}$ (rather than computing $\mathbf{S}$ first)	Remark 5	Block B
The row/column duality “things cluster by kind, people form groups”	Remark 4	Block D

Background and Motivation¶

The Iris dataset (Fisher, 1936) contains 150 iris flowers. For each flower, 4 features are measured—sepal length, sepal width, petal length, and petal width (all in centimeters)—and the species is recorded (setosa, versicolor, virginica, 50 flowers each).

The core question: the 4 features cannot all be drawn in a plane at once. At most we can pick two features for a scatter plot, but then we lose the information in the other two features. Can we find the “most telling” two-dimensional plane—not a coordinate plane spanned by any two of the original features, but a combination of directions in 4-dimensional space that retains the most variance—and pack as much of the information of all 4 features into it as possible?

This is exactly the question that PCA answers. In order, this experiment will:

observe the correlation structure among the original features, and understand why plotting an arbitrary pair of features is not ideal (Block A);
apply the SVD to the centered data and read off the principal directions and their loadings (Block B);
use a scree plot to see how the variance is distributed among the principal components, and verify the conservation of total variance (Block C);
project onto the best plane, spanned by the first 2 principal components, watch the three species separate into clusters, and use a biplot to interpret the biological meaning of the principal components (Block D);
close with the reconstruction error rank by rank, bringing Eckart–Young’s “error = sum of the squares of the tail singular values” down to real data (Block E).

Notational Conventions¶

We arrange the data matrix with samples as rows (the mainstream convention in machine learning). Vectors, matrices, and summation indices all start at 0.

Symbol	Meaning	Dimension
$m$	Number of samples (flowers)	$m = 150$
$n$	Number of features	$n = 4$
$\mathbf{X}$	Data matrix; row $k$ contains the features of flower $k$	$\mathbb{R}^{m\times n}$
$\bar{\mathbf{x}}$	Sample mean vector	$\mathbb{R}^{n}$
$\tilde{\mathbf{X}}$	Centered data matrix $\mathbf{X} - \mathbf{1}_m\bar{\mathbf{x}}^{\top}$	$\mathbb{R}^{m\times n}$
$\mathbf{S}$	Sample covariance matrix $\frac{1}{m-1}\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}}$	$\mathbb{R}^{n\times n}$ , symmetric positive semidefinite
$\sigma_i$	The $i$ -th singular value of $\tilde{\mathbf{X}}$	$\sigma_0 \ge \sigma_1 \ge \cdots$
$\mathbf{v}_i$	The $i$ -th right singular vector (the $i$ -th principal direction)	$\mathbb{R}^{n}$
$\mathbf{z}^{(k)}$	Principal component scores of flower $k$	$\mathbb{R}^{n}$
$\rho_i$	Variance share of the $i$ -th principal component, $\sigma_i^2/\sum_j\sigma_j^2$	$[0,1]$

Structure of the Experiment¶

Block	Topic	Core concepts
A	Looking at the data	The 4 original features, species labels, pairwise scatter plots and correlations
B	Centering and the SVD	The SVD of $\tilde{\mathbf{X}}$ , principal directions and loadings, correspondence with the eigenvalues of $\mathbf{S}$
C	Variance decomposition	Scree plot, cumulative explained variance, verifying the conservation of total variance
D	Best projection and clustering	Projection onto PC0–PC1, biplot, the duality “things cluster by kind, people form groups”
E	Reconstruction error	Truncated SVD reconstruction rank by rank, error = sum of the squares of the tail singular values

How to use: run the cells in order from top to bottom. Passages marked ◆ are advanced extensions and may be skipped on a first reading.

Environment Setup¶

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from sklearn.datasets import load_iris
import warnings
warnings.filterwarnings('ignore')

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

# --- Color scheme (Accent Mix) ---
C_BG    = "#F8F8F8"   # canvas background
C_GRID  = "#D6D6D6"   # grid lines
C_AXIS  = "#000000"   # axes / title text
C_V1    = "#57068C"   # violet
C_V2    = "#006385"   # deep blue
C_T1    = "#2AD2C9"   # teal
C_T2    = "#8900E1"   # Ultra Violet
C_WARN  = "#FF5D47"   # orange
C_AUX   = "#AB82C5"   # light violet

# the three species always use the same three main colors throughout the experiment
CLASS_COLORS = [C_V1, C_T1, C_WARN]   # setosa, versicolor, virginica

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

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

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value:\n{actual}")
    print(f"  Expected value: {expected_desc}")

print_header("Experiment environment initialized")
print("  Color scheme: book palette")
print("  Plotting engine: Plotly")

Block A: Looking at the Data¶

A1 | Loading the Iris Dataset¶

Each row of the Iris data matrix $\mathbf{X}\in\mathbb{R}^{150\times 4}$ is a flower, and its 4 columns are, in order, sepal length, sepal width, petal length, and petal width. The label vector $\mathbf{y}\in\{0,1,2\}^{150}$ records the species (0=setosa, 1=versicolor, 2=virginica).

PCA is an unsupervised method (the derivation of Theorem 11 does not use the labels at all), but we keep the labels anyway—not for training, but to color the scatter plots afterward, so that we can check whether the low-dimensional structure that PCA finds on its own happens to correspond to the true species.

# ============================================================
print_header("A1 | Loading the Iris dataset")
# ============================================================

iris = load_iris()
X = iris.data.astype(float)          # (150, 4)
y = iris.target                      # (150,)  0/1/2
feature_names = ["sepal length", "sepal width", "petal length", "petal width"]
class_names = ["setosa", "versicolor", "virginica"]

m, n = X.shape

print_step(1, f"Data matrix X ∈ ℝ^{{{m}×{n}}} (first 5 rows)")
print(X[:5])
print(f"\n  m = {m} flowers, n = {n} features")
print(f"  Features (columns): {feature_names}")
print(f"  Species (labels): {class_names}, with {list(np.bincount(y))} flowers each")

print_step(2, "The features have different scales (petals vs sepals)")
for j in range(n):
    print(f"  column {j} {feature_names[j]}: "
          f"mean {X[:, j].mean():6.3f}, "
          f"std. dev. {X[:, j].std(ddof=1):6.3f}, "
          f"range [{X[:, j].min():.1f}, {X[:, j].max():.1f}]")

A2 | Pairwise Scatter Plots: Why Plotting an Arbitrary Pair of Features Is Not Ideal¶

The figure below plots the 4 features against one another in pairs. We see that the petal length/width pair (bottom right) already separates setosa cleanly, but versicolor and virginica still overlap on every single pair of features; on the sepal length/width pair (top left) the three classes are even more entangled.

The root of the problem is that each subplot uses only 2 of the 4 features and throws away the other 2. The goal of PCA is not to pick the best of these subplots but to recombine all 4 features and build the two-dimensional plane in four-dimensional space that retains the most variance.

# ============================================================
print_header("A2 | Pairwise scatter plots of the four features")
# ============================================================

pairs = [(0, 1), (0, 2), (1, 3), (2, 3)]   # pick 4 representative pairs of features
fig = make_subplots(
    rows=2, cols=2,
    subplot_titles=[f"{feature_names[a]} × {feature_names[b]}" for a, b in pairs],
    horizontal_spacing=0.12, vertical_spacing=0.16
)

for idx, (a, b) in enumerate(pairs):
    r, c = idx // 2 + 1, idx % 2 + 1
    for cls in range(3):
        mask = y == cls
        fig.add_trace(go.Scatter(
            x=X[mask, a], y=X[mask, b], mode='markers',
            marker=dict(size=6, color=CLASS_COLORS[cls], opacity=0.75,
                        line=dict(width=0.3, color=C_AXIS)),
            name=class_names[cls], legendgroup=class_names[cls],
            showlegend=(idx == 0)
        ), row=r, col=c)
    fig.update_xaxes(title_text=feature_names[a], gridcolor=C_GRID, row=r, col=c)
    fig.update_yaxes(title_text=feature_names[b], gridcolor=C_GRID, row=r, col=c)

fig.update_layout(
    title=dict(text="Pairwise scatter plots of the original features: versicolor and virginica always partly overlap",
               font=dict(size=15, color=C_AXIS)),
    plot_bgcolor=C_BG, paper_bgcolor=C_BG, height=720, width=860,
    font=dict(family="Microsoft JhengHei, SimHei, sans-serif"),
    legend=dict(orientation="h", yanchor="bottom", y=1.06, xanchor="center", x=0.5)
)
fig.show()
print("  Observation: the petal features discriminate best, but no single pair of features fully separates the three classes")

Block B: Centering and the SVD¶

B1 | Center, Then Apply the SVD to $\tilde{\mathbf{X}}$ ¶

Following the standard procedure of §11.3.1: first center to obtain $\tilde{\mathbf{X}} = \mathbf{X} - \mathbf{1}_m\bar{\mathbf{x}}^{\top}$ , then apply the singular value decomposition directly to $\tilde{\mathbf{X}}$ :

\tilde{\mathbf{X}} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top} = \sum_{i=0}^{r-1}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^{\top} = \sum_{i=0}^{r-1}\sigma_i\,|\mathbf{u}_i\rangle\langle\mathbf{v}_i|

(1)

The right singular vectors $\mathbf{v}_0,\mathbf{v}_1,\ldots$ are the principal directions 0, 1, …; the corresponding $\sigma_i^2/(m-1)$ is the variance explained by each principal component.

# ============================================================
print_header("B1 | Center, then apply the SVD")
# ============================================================

print_step(1, "Centering: subtract the mean vector so that the centroid of the sample cloud moves to the origin")
x_bar = X.mean(axis=0)
X_tilde = X - x_bar
print(f"  Mean vector x̄ = {x_bar}")
print(f"  Column means after centering = {X_tilde.mean(axis=0)} (should all be 0)")

print_step(2, "Apply the reduced SVD to X̃: X̃ = U Σ V^T")
U, sigma, Vt = np.linalg.svd(X_tilde, full_matrices=False)
V = Vt.T                      # the columns are the principal directions v_0, v_1, ...
r = np.sum(sigma > 1e-10)
print(f"  Singular values σ = {sigma}")
print(f"  Numerical rank r = {r}")

print_step(3, "The right singular vectors are the principal directions (the columns of V)")
for i in range(n):
    print(f"  v_{i} = {V[:, i]}")

B2 | Agreement of the Two Paths + Interpreting the Principal Component Loadings¶

The first part verifies that the SVD and the eigendecomposition of the covariance matrix give the same variances and directions: the eigenvalues of $\mathbf{S}$ should equal $\sigma_i^2/(m-1)$ , and the eigenvectors should be parallel to $\mathbf{v}_i$ (the signs may be opposite; this is an inherent freedom of the SVD and does not affect the subspaces spanned by the principal components).

The second part reads $\mathbf{V}$ as a loading table: the $j$ -th component of $\mathbf{v}_i$ is “the weight of the $j$ -th original feature in the $i$ -th principal component.” Once you can read the loadings, you can read the biological meaning of each principal component.

# ============================================================
print_header("B2 | SVD vs eigendecomposition of the covariance matrix + loading table")
# ============================================================

print_step(1, "Check: the eigenvalues of S = σ²/(m-1), and the eigenvectors are parallel to v_i")
S = (X_tilde.T @ X_tilde) / (m - 1)          # covariance matrix (computed here only for the check; the proper procedure does not form it)
eigval, eigvec = np.linalg.eigh(S)
eigval = eigval[::-1]                          # eigh sorts from smallest to largest; reverse to largest first
eigvec = eigvec[:, ::-1]
compare_print("Eigenvalues vs σ²/(m-1)", np.round(eigval, 6),
              f"{np.round(sigma**2 / (m - 1), 6)}")
align = np.abs(np.sum(eigvec * V, axis=0))     # |inner product| of each pair of directions; close to 1 if parallel
print(f"  |⟨eigenvector, v_i⟩| = {np.round(align, 6)} (close to 1 means the directions are parallel)")

print_step(2, "Loading table of the principal components: each v_i gives weights on the 4 original features")
lam = sigma**2 / (m - 1)
header = "  " + " " * 8 + "".join(f"{nm:>10}" for nm in ["PC0", "PC1", "PC2", "PC3"])
print(header)
for j in range(n):
    row = "".join(f"{V[j, i]:10.3f}" for i in range(n))
    print(f"  {feature_names[j]:<12}{row}")
print(f"\n  Variance of each principal component λ_i = σ²/(m-1) = {np.round(lam, 4)}")

print_step(3, "Interpretation (by the absolute values of the weights, regardless of sign)")
print("  PC0: petal length, petal width, and sepal length have large weights of the same sign, sepal width is near 0 → the \"overall flower size\" axis")
print("  PC1: sepal length and sepal width have the largest weights → a secondary axis reflecting \"sepal size\"")

Block C: Variance Decomposition and the Scree Plot¶

C1 | How Much Variance Does Each Principal Component Explain?¶

We bring the two formulas of §11.3.1 down to Iris. For $\boldsymbol{\Pi}_k = \mathbf{V}_k\mathbf{V}_k^{\top}$ (the projection onto the first $k$ principal components),

V_{\text{approx}}(\boldsymbol{\Pi}_k) = \frac{1}{m-1}\sum_{i=0}^{k-1}\sigma_i^2,\qquad V_{\text{err}}(\boldsymbol{\Pi}_k) = \frac{1}{m-1}\sum_{i=k}^{r-1}\sigma_i^2

(2)

and Theorem 10 guarantees that whatever the value of $k$ , the two always add up to the total variance $V_{\text{total}}$ . We define the variance share of the $i$ -th principal component as

\rho_i = \frac{\sigma_i^2}{\sum_{j}\sigma_j^2}

(3)

and draw it as a scree plot: the bars are the $\rho_i$ of the principal components, and the line is the cumulative share $\sum_{j\le i}\rho_j$ . The scree plot of Iris shows an extremely steep “elbow”—the first bar alone takes more than 90% of the variance, and this is the fundamental reason why PCA can compress 4 dimensions into 2 with almost no distortion.

# ============================================================
print_header("C1 | Variance shares and conservation of total variance")
# ============================================================

print_step(1, "Compute the variance share ρ of each principal component and the cumulative share")
var_each = sigma**2 / (m - 1)               # variance of each principal component
rho = sigma**2 / (sigma**2).sum()            # variance share
cum = np.cumsum(rho)                          # cumulative share
V_total = np.trace(S)                         # total variance = tr(S)
for i in range(n):
    print(f"  PC{i}: σ={sigma[i]:7.3f}  variance={var_each[i]:7.4f}  "
          f"ρ={rho[i]:7.2%}  cumulative={cum[i]:7.2%}")
print(f"\n  Total variance V_total = tr(S) = {V_total:.4f}")
print(f"  The first 2 principal components together explain {cum[1]:.2%} of the variance → reducing to 2 dimensions causes almost no distortion")

print_step(2, "Verify conservation: for every k, V_approx(k) + V_err(k) = V_total")
print("   k │  V_approx     V_err       sum      V_total")
print("  ───┼─────────────────────────────────────────────")
for k in range(n + 1):
    V_approx = var_each[:k].sum()
    V_err = var_each[k:].sum()
    print(f"   {k} │ {V_approx:9.4f}  {V_err:9.4f}  {V_approx + V_err:9.4f}   {V_total:9.4f}")
compare_print("Conservation check", f"the \"sum\" column for every k", "always equals V_total (independent of k)")

# ============================================================
print_header("C2 | Scree plot")
# ============================================================

pc_labels = [f"PC{i}" for i in range(n)]
fig = make_subplots(specs=[[{"secondary_y": True}]])

# --- Bars: variance share of each principal component ---
fig.add_trace(go.Bar(
    x=pc_labels, y=rho, name="Individual share ρ_i",
    marker=dict(color=[C_V1, C_V2, C_T1, C_AUX]),
    text=[f"{v:.1%}" for v in rho], textposition='outside'
), secondary_y=False)

# --- Line: cumulative share ---
fig.add_trace(go.Scatter(
    x=pc_labels, y=cum, name="Cumulative share", mode='lines+markers',
    line=dict(color=C_WARN, width=2.5),
    marker=dict(size=9, color=C_WARN)
), secondary_y=True)

# --- 95% reference line ---
fig.add_hline(y=0.95, line=dict(color=C_AXIS, width=1, dash="dot"),
              annotation_text="95%", annotation_position="bottom right",
              secondary_y=True)

fig.update_layout(
    title=dict(text="Iris scree plot: PC0 alone takes 92.5% of the variance, the first 2 principal components 97.8%",
               font=dict(size=15, color=C_AXIS)),
    plot_bgcolor=C_BG, paper_bgcolor=C_BG, height=460, width=720,
    font=dict(family="Microsoft JhengHei, SimHei, sans-serif"),
    legend=dict(orientation="h", yanchor="bottom", y=1.02, xanchor="center", x=0.5),
    bargap=0.45
)
fig.update_yaxes(title_text="Individual variance share", range=[0, 1.05], gridcolor=C_GRID, secondary_y=False)
fig.update_yaxes(title_text="Cumulative variance share", range=[0, 1.05], showgrid=False, secondary_y=True)
fig.show()
print("  The elbow is extremely steep: the first bar dominates, which is the fundamental reason PCA can reduce the dimension drastically")

Block D: Projecting onto the Best Plane and Watching the Clusters Appear¶

D1 | Principal Component Scores: Projecting the 4-Dimensional Flowers onto PC0–PC1¶

The principal component scores of a flower are simply its coordinates in the new coordinate system (the principal directions). Computing them for the whole batch takes a single matrix multiplication:

\mathbf{Z} = \tilde{\mathbf{X}}\mathbf{V}\in\mathbb{R}^{m\times n},\qquad \mathbf{z}^{(k)} = \mathbf{V}^{\top}\tilde{\mathbf{x}}^{(k)}

(4)

Column 0 of $\mathbf{Z}$ holds the PC0 scores of all the flowers, and column 1 the PC1 scores. Plotting the first 2 columns as a scatter plot gives the “best two-dimensional projection” of §11.3—among all two-dimensional planes it retains the most variance (Theorem 11).

The key point: PCA does not use the species labels at all; it blindly looks for the plane of maximal variance. Yet as soon as the figure below is colored by the true species, three groups emerge clearly—setosa is flung far away, while versicolor and virginica line up as two segments along PC0 whose boundaries barely touch. This shows that the species differences in Iris are essentially a one-dimensional “flower size” gradient (PC0), which is exactly the direction of maximal variance that PCA picks up automatically.

# ============================================================
print_header("D1 | Projecting onto the best plane PC0–PC1")
# ============================================================

print_step(1, "One matrix multiplication gives all the principal component scores Z = X̃ V")
Z = X_tilde @ V
print(f"  Score matrix Z ∈ ℝ^{{{Z.shape[0]}×{Z.shape[1]}}} (first 3 rows)")
print(Z[:3])

print_step(2, "Draw the PC0–PC1 scatter plot (colored by the true species, although PCA does not use the labels)")
fig = go.Figure()
for cls in range(3):
    mask = y == cls
    fig.add_trace(go.Scatter(
        x=Z[mask, 0], y=Z[mask, 1], mode='markers',
        marker=dict(size=8, color=CLASS_COLORS[cls], opacity=0.8,
                    line=dict(width=0.4, color=C_AXIS)),
        name=class_names[cls]
    ))
fig.add_hline(y=0, line=dict(color=C_GRID, width=1, dash="dot"))
fig.add_vline(x=0, line=dict(color=C_GRID, width=1, dash="dot"))
fig.update_layout(
    title=dict(text=f"Projection of Iris onto the best two-dimensional plane (PC0–PC1, retaining {cum[1]:.1%} of the variance)",
               font=dict(size=15, color=C_AXIS)),
    xaxis_title=f"PC0 ({rho[0]:.1%}, flower size axis)",
    yaxis_title=f"PC1 ({rho[1]:.1%}, sepal size axis)",
    plot_bgcolor=C_BG, paper_bgcolor=C_BG, height=520, width=720,
    font=dict(family="Microsoft JhengHei, SimHei, sans-serif"),
    legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01)
)
fig.update_xaxes(gridcolor=C_GRID)
fig.update_yaxes(gridcolor=C_GRID)
fig.show()
print("  setosa lies far from the other two groups along PC0; versicolor and virginica line up as two segments along PC0 that barely touch")

D2 | Biplot: Overlaying Sample Scores and Feature Loadings in One Figure¶

A biplot draws two kinds of things at once in the same PC0–PC1 plane, corresponding to the row/column duality of Remark 4:

Points = sample scores (the first two columns of $\mathbf{Z}$ ), that is, “people form groups”—flowers of the same species gather together.
Arrows = the original feature axes projected onto the PC plane, that is, “things cluster by kind”—letting us see which features point in the same direction.

Where do the arrows come from? The $j$ -th original feature axis is the standard basis vector $\mathbf{e}_j\in\mathbb{R}^4$ ; rewriting it in the new coordinates given by the principal components, its coordinates are $(\mathbf{v}_0^{\top}\mathbf{e}_j,\ \mathbf{v}_1^{\top}\mathbf{e}_j) = (V_{j0},\ V_{j1})$ , that is, the first two components of row $j$ of $\mathbf{V}$ . This is precisely the orthogonal projection of $\mathbf{e}_j$ onto the PC0–PC1 plane. To make them visible against the cloud of points, all the arrows are then multiplied by one and the same magnification factor (which changes only their lengths, not their directions). How to read the plot:

Arrow length (before projection $\|\mathbf{e}_j\|=1$ ): the closer to 1, the more nearly the whole feature lies in the plane and is faithfully represented; the shorter, the more of it “sticks out of the plane,” hidden in the discarded PC2/PC3.
Angle with the PC0 axis: whether the feature leans toward PC0 or PC1.
Angle between two arrows: a small angle $\Rightarrow$ the two features are positively correlated in this plane.

# ============================================================
print_header("D2 | Biplot: scores + loadings overlaid")
# ============================================================

print_step(1, "Arrows = the original feature axes e_j projected onto the PC plane (first two components of row j of V)")
proj = V[:, :2]                                   # row j = the projection of e_j onto PC0-PC1
arrow_len = np.hypot(proj[:, 0], proj[:, 1])      # length of each arrow (≤ 1)
# common magnification factor: stretch the longest arrow to ~85% of the range of the points (changes only lengths, not directions)
scale = 0.85 * np.max(np.abs(Z[:, :2])) / arrow_len.max()
arrows = proj * scale
for j in range(n):
    print(f"  {feature_names[j]:<12} projection=({proj[j,0]:+.3f}, {proj[j,1]:+.3f})  "
          f"length={arrow_len[j]:.3f}  ({'faithful' if arrow_len[j] > 0.7 else 'partly sticks out of the plane'})")

print_step(2, "Draw the biplot (equal aspect ratio, so that the angles can be trusted)")
fig = go.Figure()
# --- Sample scores (people form groups) ---
for cls in range(3):
    mask = y == cls
    fig.add_trace(go.Scatter(
        x=Z[mask, 0], y=Z[mask, 1], mode='markers',
        marker=dict(size=7, color=CLASS_COLORS[cls], opacity=0.55,
                    line=dict(width=0.3, color=C_AXIS)),
        name=class_names[cls]
    ))
# --- Feature loading arrows (things cluster by kind) ---
for j in range(n):
    fig.add_annotation(
        x=arrows[j, 0], y=arrows[j, 1], ax=0, ay=0,
        xref="x", yref="y", axref="x", ayref="y",
        showarrow=True, arrowhead=3, arrowsize=1.4, arrowwidth=2, arrowcolor=C_AXIS
    )
    fig.add_trace(go.Scatter(
        x=[arrows[j, 0] * 1.10], y=[arrows[j, 1] * 1.10],
        mode='text', text=[feature_names[j]],
        textfont=dict(size=12, color=C_AXIS), showlegend=False
    ))
fig.add_hline(y=0, line=dict(color=C_GRID, width=1, dash="dot"))
fig.add_vline(x=0, line=dict(color=C_GRID, width=1, dash="dot"))
fig.update_layout(
    title=dict(text="Iris biplot: sample scores (people form groups) + projected feature axes (things cluster by kind)",
               font=dict(size=15, color=C_AXIS)),
    xaxis_title="PC0", yaxis_title="PC1",
    plot_bgcolor=C_BG, paper_bgcolor=C_BG, height=560, width=740,
    font=dict(family="Microsoft JhengHei, SimHei, sans-serif"),
    legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01)
)
fig.update_xaxes(gridcolor=C_GRID)
fig.update_yaxes(gridcolor=C_GRID, scaleanchor="x", scaleratio=1)
fig.show()
print("  The petal length arrow is the longest and closest to horizontal → it is the main support of the PC0 \"size axis\"")
print("  The sepal width arrow points almost straight down → PC1 is mainly distinguished by it (sepal shape)")
print("  The petal width arrow is the shortest (0.37) → most of it lies in the discarded PC2/PC3, so it is represented most weakly in this plane")

D3 | Interactive: Choose Any Two Principal Components as Coordinate Axes¶

When PCA compresses 4 dimensions into 2, it takes PC0–PC1 by default (retaining the most variance). Sometimes, however, a minor principal component hides another kind of structure (the Swiss roll of §11.3.2 is exactly such a counterexample). The two drop-down menus below let you choose freely which two principal components to use for the $x$ - and $y$ -axes.

Here we use plotly’s native drop-downs (updatemenus). They are controlled purely on the front end and do not depend on a running kernel, so they update instantly both in JupyterLab and in the statically built HTML book. Try comparing PC0–PC1 (the default, with the best separation) with PC2–PC3 (residual directions, where the three classes overlap almost completely), and see that on Iris “retaining the most variance” and “discriminating best” happen to coincide—but this coincidence is not a matter of course; it happens because the species differences happen to line up with the direction of maximal variance.

# ============================================================
print_header("D3 | Interactive: choose the principal component axes (plotly native drop-downs)")
# ============================================================

masks = [y == cls for cls in range(3)]

fig = go.Figure()
for cls in range(3):
    fig.add_trace(go.Scatter(
        x=Z[masks[cls], 0], y=Z[masks[cls], 1], mode='markers',
        marker=dict(size=8, color=CLASS_COLORS[cls], opacity=0.8,
                    line=dict(width=0.4, color=C_AXIS)),
        name=class_names[cls]
    ))

# Each x-button changes only the x of the 3 traces (and the x-axis title); each y-button changes only y.
# The two menus restyle independently without overwriting each other, so any (PCi, PCj) can be combined freely.
def x_button(i):
    return dict(method="update", label=f"PC{i}",
                args=[{"x": [Z[masks[cls], i] for cls in range(3)]},
                      {"xaxis.title.text": f"PC{i} ({rho[i]:.1%})"}])

def y_button(j):
    return dict(method="update", label=f"PC{j}",
                args=[{"y": [Z[masks[cls], j] for cls in range(3)]},
                      {"yaxis.title.text": f"PC{j} ({rho[j]:.1%})"}])

fig.update_layout(
    title=dict(text="Interactive principal component projection (switch the axes with the drop-down menus)",
               font=dict(size=15, color=C_AXIS)),
    xaxis=dict(title=f"PC0 ({rho[0]:.1%})", gridcolor=C_GRID),
    yaxis=dict(title=f"PC1 ({rho[1]:.1%})", gridcolor=C_GRID),
    plot_bgcolor=C_BG, paper_bgcolor=C_BG, height=560, width=740,
    font=dict(family="Microsoft JhengHei, SimHei, sans-serif"),
    legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
    margin=dict(t=110),
    updatemenus=[
        dict(buttons=[x_button(i) for i in range(n)], direction="down",
             showactive=True, active=0,
             x=0.12, xanchor="left", y=1.18, yanchor="top",
             bgcolor="#FFFFFF", bordercolor=C_GRID),
        dict(buttons=[y_button(j) for j in range(n)], direction="down",
             showactive=True, active=1,
             x=0.40, xanchor="left", y=1.18, yanchor="top",
             bgcolor="#FFFFFF", bordercolor=C_GRID),
    ],
    annotations=[
        dict(text="x-axis:", x=0.12, xref="paper", y=1.205, yref="paper",
             xanchor="right", showarrow=False, font=dict(size=12, color=C_AXIS)),
        dict(text="y-axis:", x=0.40, xref="paper", y=1.205, yref="paper",
             xanchor="right", showarrow=False, font=dict(size=12, color=C_AXIS)),
    ]
)
fig.show()
print("  The default PC0–PC1 separates best; switch to PC2–PC3 to see the residual directions, where the three classes overlap almost completely")

Block E: Reconstruction Rank by Rank and the Eckart–Young Error Formula¶

E1 | Truncated SVD Reconstruction: The Error Is Exactly the Sum of the Squares of the Tail Singular Values¶

The last piece of the puzzle of the Eckart–Young theorem (Theorem 8): reconstruct the data from the first $k$ principal components,

\tilde{\mathbf{X}}_k = \sum_{i=0}^{k-1}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^{\top} = \tilde{\mathbf{X}}\boldsymbol{\Pi}_k

(5)

Its reconstruction error has a closed form,

\|\tilde{\mathbf{X}} - \tilde{\mathbf{X}}_k\|_F^2 = \sum_{i=k}^{r-1}\sigma_i^2

(6)

that is, “the sum of the squares of the singular values of the discarded principal components.” In this block we reconstruct for $k=0,1,2,3,4$ one at a time and set the measured Frobenius error side by side with the tail sum of squares predicted by the formula—the most direct numerical verification of the whole theory of §11.3 on real data. We also plot the relative reconstruction error $\|\tilde{\mathbf{X}}-\tilde{\mathbf{X}}_k\|_F / \|\tilde{\mathbf{X}}\|_F$ as a curve decreasing in $k$ , confirming by eye that at $k=2$ the error is already negligibly small.

# ============================================================
print_header("E1 | Rank-by-rank reconstruction error vs the Eckart–Young formula")
# ============================================================

print_step(1, "Reconstruct X̃_k for k=0..4 and compare the measured error with the formula")
fro_total = np.linalg.norm(X_tilde, 'fro')
ks, rel_err = [], []
print("   k │ measured ‖X̃-X̃_k‖²_F  formula Σσ²(tail)  rel. error")
print("  ───┼──────────────────────────────────────────────")
for k in range(n + 1):
    Xk = U[:, :k] @ np.diag(sigma[:k]) @ Vt[:k, :]   # truncated SVD reconstruction
    err2_measured = np.linalg.norm(X_tilde - Xk, 'fro') ** 2
    err2_formula = (sigma[k:] ** 2).sum()             # sum of the squares of the tail singular values
    rel = np.sqrt(err2_measured) / fro_total
    ks.append(k); rel_err.append(rel)
    print(f"   {k} │ {err2_measured:14.6f}  {err2_formula:14.6f}  {rel:10.4%}")
compare_print("Eckart–Young check", "the \"measured\" column", "equals the \"formula Σσ²(tail)\" column row by row")

print_step(2, "Plot the relative reconstruction error decreasing with k")
fig = go.Figure()
fig.add_trace(go.Scatter(
    x=ks, y=rel_err, mode='lines+markers',
    line=dict(color=C_V1, width=2.5), marker=dict(size=10, color=C_V1),
    name="Relative reconstruction error"
))
fig.add_hline(y=0.15, line=dict(color=C_WARN, width=1, dash="dot"),
              annotation_text="k=2 (≈15%)", annotation_position="top right")
fig.update_layout(
    title=dict(text="Relative reconstruction error ‖X̃-X̃_k‖_F / ‖X̃‖_F: about 15% at k=2",
               font=dict(size=15, color=C_AXIS)),
    xaxis_title="Number of principal components retained k", yaxis_title="Relative Frobenius error",
    plot_bgcolor=C_BG, paper_bgcolor=C_BG, height=440, width=700,
    font=dict(family="Microsoft JhengHei, SimHei, sans-serif"),
    showlegend=False
)
fig.update_xaxes(gridcolor=C_GRID, dtick=1)
fig.update_yaxes(gridcolor=C_GRID, range=[-0.02, max(rel_err) * 1.05])
fig.show()
print("  At k=2 the relative error is about 15%, and the retained variance share = 1 - 0.15² ≈ 97.8%, consistent with the scree plot")

E2 (Exercise) | How Many Principal Components Are Needed to Reach 99% Cumulative Variance?¶

Exercise 1 (The Minimal Number of Retained Dimensions for Iris)

Using the precomputed rho (the variance share of each principal component), write one line of code that finds the smallest $k$ for which the cumulative variance share of the first $k$ principal components first reaches or exceeds 99%. Then use this $k$ to look back at the relative reconstruction error of Block E1 and confirm that the two agree (retaining 99% of the cumulative variance $\Leftrightarrow$ a relative error within about $\sqrt{1-0.99}\approx 10\%$ ).

Solution to Exercise 1

k99 = int(np.searchsorted(np.cumsum(rho), 0.99) + 1)
print(f"Reaching 99% cumulative variance requires k = {k99} principal components")
print(f"Cumulative share at this k = {np.cumsum(rho)[k99-1]:.4%}")
# For comparison: retained variance share = 1 - (relative error)², so relative error = sqrt(1 - cumulative share)
print(f"Corresponding relative reconstruction error ≈ {np.sqrt(1 - np.cumsum(rho)[k99-1]):.4%}")

For Iris, the first 3 principal components together account for about 99.5%, so $k=3$ ; the first 2 account for 97.8% (relative error about 15%), which does not yet reach 99%. This shows that if extremely high fidelity is required, Iris still needs to keep principal component 2 (PC2)—but for the vast majority of visualization and classification purposes, 2 dimensions are more than enough.

Summary of the Experiment¶

This experiment runs the PCA theory of §11.3 in full on the Iris dataset. When you have finished it, you should be able to:

Carry out the standard PCA procedure: centering → SVD of $\tilde{\mathbf{X}}$ → take the right singular vectors as the principal components, and explain why we do not form the covariance matrix $\mathbf{S}$ first (the condition number is squared, Remark 5).
Verify that the SVD and the eigendecomposition of $\mathbf{S}$ give the same variances ( $\lambda_i = \sigma_i^2/(m-1)$ ) and parallel directions.
Interpret the loading table of the principal components: state that PC0 of Iris is the “flower size axis” and PC1 the “sepal size axis.”
Draw and read a scree plot, and point out that the elbow for Iris is extremely steep (PC0 accounts for 92.5%), which is the basis for a drastic reduction in dimension.
Verify numerically the conservation of total variance: for every $k$ , $V_{\text{approx}}(k)+V_{\text{err}}(k)=V_{\text{total}}$ (Theorem 10).
Project and watch PCA separate the three species in the PC0–PC1 plane without labels (Theorem 11), and use a biplot to overlay “things cluster by kind” (feature loadings) and “people form groups” (sample scores) in a single figure (Remark 4).
Verify the Eckart–Young error formula: the rank-by-rank reconstruction error $\|\tilde{\mathbf{X}}-\tilde{\mathbf{X}}_k\|_F^2 = \sum_{i\ge k}\sigma_i^2$ (Theorem 8).

Concept	Key formula	Values for Iris
Centering + SVD	$\tilde{\mathbf{X}} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top}$	$\sigma = (25.10,\ 6.01,\ 3.41,\ 1.89)$
Variance of the principal components	$\lambda_i = \sigma_i^2/(m-1)$	—
Variance share	$\rho_i = \sigma_i^2/\sum_j\sigma_j^2$	$(92.5\%,\ 5.3\%,\ 1.7\%,\ 0.5\%)$
Conservation of total variance	$V_{\text{approx}}(k)+V_{\text{err}}(k)=V_{\text{total}}$	Independent of $k$
Principal component scores	$\mathbf{Z} = \tilde{\mathbf{X}}\mathbf{V}$	The first 2 columns are plotted
Reconstruction error	$\|\tilde{\mathbf{X}}-\tilde{\mathbf{X}}_k\|_F^2 = \sum_{i\ge k}\sigma_i^2$	Relative error $\approx 15\%$ at $k=2$ (retaining $97.8\%$ of the variance)

Concept Map¶

This concept map ties the whole experiment into a single thread: centering → SVD is the engine, and it produces two outputs, the singular values (which rank the variances by priority) and the right singular vectors (the best coordinate axes). The former support the scree plot, the conservation law, and the Eckart–Young error formula; the latter, through the scores $\mathbf{Z}=\tilde{\mathbf{X}}\mathbf{V}$ , project the samples onto the best plane and meet the feature loadings in the biplot.

Looking back: this experiment brings the PCA theory of §11.3 (conservation + best projection) down to real data, reusing the centering and total-variance tools of Experiment 5. Looking ahead: put a different matrix into the same SVD engine and the remaining applications of this chapter grow out of it—replace the matrix with term-document counts and you get latent semantic analysis (LSA); use the singular values as a threshold for regularization and you get the Moore–Penrose pseudoinverse and ridge regression of §11.4. Here Iris is only the gentlest introductory example; the real tension (high variance is not necessarily meaningful) was already foreshadowed by the Swiss roll of §11.3.2.