PCA of the Iris Dataset

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

The iris flower dataset is \(150\) measurements of three species of iris, four features each. This example is loosely modeled on scikit-learn’s own PCA-on-iris and PCA-vs-LDA galleries, and uses PCA purely as a way to look at \(4\)-dimensional data in \(2\) and \(3\) dimensions.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from sklearn.datasets import load_iris

iris = load_iris()
X_raw = iris.data                                     # (150, 4), centimeters
species = iris.target                                 # 0, 1, 2
species_names = iris.target_names
feature_names = [name.replace(" (cm)", "") for name in iris.feature_names]
n, d = X_raw.shape

print(f"{n} flowers, {d} features: {feature_names}")
print(f"species: {list(species_names)}")

150 flowers, 4 features: ['sepal length', 'sepal width', 'petal length', 'petal width']
species: [np.str_('setosa'), np.str_('versicolor'), np.str_('virginica')]

# A muted ink colour and a translucent grid read correctly against both the
# light and the dark version of this site, so the figures need no theme switch.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE, GREEN = "#2a78d6", "#eb6834", "#1baf7a"
SPECIES_COLOR = {0: BLUE, 1: ORANGE, 2: GREEN}


def style(fig, height=420, title=None, showlegend=None):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
        margin=dict(l=70, r=30, t=60 if title else 30, b=55),
        height=height,
        showlegend=showlegend,
        hovermode="closest",
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig


def species_traces(x, y, showlegend):
    """One Scatter trace per species, for a 2D feature or score plot."""
    return [
        go.Scatter(
            x=x[species == s], y=y[species == s], mode="markers",
            name=name, legendgroup=name, showlegend=showlegend,
            marker=dict(color=SPECIES_COLOR[s], size=6, opacity=0.85,
                        line=dict(color="white", width=0.5)),
        )
        for s, name in enumerate(species_names)
    ]

The data

Stacking the \(150\) flowers as rows gives the data matrix \[ X_{\text{raw}}= \begin{pmatrix} - & \mathbf{x}_1^\top & - \\ & \vdots & \\ - & \mathbf{x}_{150}^\top & - \end{pmatrix} \in\mathbb{R}^{150\times4}, \] with columns sepal length, sepal width, petal length, and petal width, all in centimeters. Every pair of features already shows some separation between the three species, with the petal measurements doing most of the work.

pairs = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]
fig = make_subplots(
    rows=2, cols=3, horizontal_spacing=0.07, vertical_spacing=0.16,
    subplot_titles=[f"{feature_names[j]} vs {feature_names[i]}" for i, j in pairs],
)
for idx, (i, j) in enumerate(pairs):
    row, col = idx // 3 + 1, idx % 3 + 1
    for trace in species_traces(X_raw[:, j], X_raw[:, i], showlegend=(idx == 0)):
        fig.add_trace(trace, row=row, col=col)
fig.update_annotations(font=dict(size=12))
style(fig, height=520, title="The 150 x 4 iris data matrix, one panel per pair of features",
     showlegend=True)
fig.update_layout(legend=dict(orientation="h", y=-0.14, x=0.32))
fig.show()

Figure 1: All six pairs of the four raw iris features, colored by species.

PCA as a visualization tool

Four features cannot be plotted directly, and each panel of Figure 1 is only one axis-aligned shadow of the data. PCA instead picks the \(2\)- or \(3\)-dimensional view that keeps the most variance. Here, we project the centered data onto the top \(k\) principal directions and plot the scores \(Z=XV_k\).

from sklearn.decomposition import PCA

pca = PCA(n_components=3).fit(X_raw)       # centers the data internally
Z = pca.transform(X_raw)                   # (150, 3) scores
variance_ratio = pca.explained_variance_ratio_

print("explained variance ratio:", np.round(variance_ratio, 4))
print(f"R_2 = {variance_ratio[:2].sum():.4f},  R_3 = {variance_ratio.sum():.4f}")

explained variance ratio: [0.9246 0.0531 0.0171]
R_2 = 0.9777,  R_3 = 0.9948

PCA in 3D

fig = go.Figure()
for s, name in enumerate(species_names):
    mask = species == s
    fig.add_trace(go.Scatter3d(
        x=Z[mask, 0], y=Z[mask, 1], z=Z[mask, 2], mode="markers", name=name,
        marker=dict(color=SPECIES_COLOR[s], size=4, opacity=0.85,
                    line=dict(color="white", width=0.3)),
    ))
fig.update_scenes(
    xaxis_title=f"z1 (PC1), {variance_ratio[0]:.1%}",
    yaxis_title=f"z2 (PC2), {variance_ratio[1]:.1%}",
    zaxis_title=f"z3 (PC3), {variance_ratio[2]:.1%}",
)
style(fig, height=560, title="Iris scores, top 3 principal components", showlegend=True)
fig.show()

Figure 2: Iris scores in the top-3 principal subspace, colored by species.

PCA in 2D

Dropping to the first two scores keeps most of the total variance (the value of \(R_2\) is printed above), and already separates all three species almost as cleanly as the full 3D plot.

fig = go.Figure()
for trace in species_traces(Z[:, 0], Z[:, 1], showlegend=True):
    fig.add_trace(trace)
fig.update_xaxes(title=f"z1 (PC1), {variance_ratio[0]:.1%} of variance")
fig.update_yaxes(title=f"z2 (PC2), {variance_ratio[1]:.1%} of variance")
style(fig, height=460, title="Iris scores, top 2 principal components")
fig.show()

Figure 3: Iris scores in the top-2 principal subspace, colored by species.

Compared to the raw feature pairs in Figure 1, the 2D principal subspace is the single best 2-dimensional view of the data in the sense of best-fitting subspaces come from the SVD: no other 2-dimensional projection of the centered data captures more of the total variance, or equivalently, leaves a smaller total squared reconstruction error.