import plotly.io as pio; pio.renderers.default = "notebook_connected"

Lecture 03: PCA and Linear Algebra Review – CS 189, Fall 2026

Demonstration: making sense of congressional votes.

We use PCA today before we understand it. By the end of this notebook we will have taken a 441 x 41 table of votes, reduced it to 441 x 2, and recovered a property of Congress that was never supplied to the algorithm. The remainder of the lecture explains why this works.

import numpy as np
import pandas as pd
import yaml
from datetime import datetime
import plotly.express as px
# Uncomment for HTML Export
import plotly.io as pio
pio.renderers.default = "notebook_connected"

Congressional Vote Records

Let's examine how the House of Representatives (of the 116th Congress, 1st session) voted in the month of September 2019.

From the U.S. Senate website:

Roll call votes occur when a representative or senator votes "yea" or "nay," so that the names of members voting on each side are recorded. A voice vote is a vote in which those in favor or against a measure say "yea" or "nay," respectively, without the names or tallies of members voting on each side being recorded.

The data, compiled from ProPublica source, is a "skinny" table of data where each record is a single vote by a member across any roll call in the 116th Congress, 1st session, as downloaded in February 2020. The member of the House, whom we'll call legislator, is denoted by their bioguide alphanumeric ID in http://bioguide.congress.gov/.

votes = pd.read_csv('data/votes.csv')
votes = votes.astype({"roll call": str})
votes

	chamber	session	roll call	member	vote
0	House	1	555	A000374	Not Voting
1	House	1	555	A000370	Yes
2	House	1	555	A000055	No
3	House	1	555	A000371	Yes
4	House	1	555	A000372	No
...	...	...	...	...	...
17823	House	1	515	Y000062	Yes
17824	House	1	515	Y000065	No
17825	House	1	515	Y000033	Yes
17826	House	1	515	Z000017	Yes
17827	House	1	515	P000197	Speaker

17828 rows × 5 columns

votes['vote'].value_counts()

vote
Yes           10373
No             6845
Not Voting      567
Speaker          41
Present           2
Name: count, dtype: int64

This is a "skinny" table, with one row per (member, roll call). To treat each legislator as a datapoint, we pivot so that each row is a legislator and each column is a roll call. We record a 1 for a Yes vote and a 0 otherwise.

def was_yes(s):
    return 1 if s.iloc[0] == "Yes" else 0

vote_pivot = votes.pivot_table(index='member',
                               columns='roll call',
                               values='vote',
                               aggfunc=was_yes,
                               fill_value=0)
print(vote_pivot.shape)
vote_pivot.head()

(441, 41)

roll call	515	516	517	518	519	520	521	522	523	524	...	546	547	548	549	550	551	552	553	554	555
member
A000055	1	0	0	0	1	1	0	1	1	1	...	0	0	1	0	0	1	0	0	1	0
A000367	0	0	0	0	0	0	0	0	0	0	...	0	1	1	1	1	0	1	1	0	1
A000369	1	1	0	0	1	1	0	1	1	1	...	0	0	1	0	0	1	0	0	1	0
A000370	1	1	1	1	1	0	1	0	0	0	...	1	1	1	1	1	0	1	1	1	1
A000371	1	1	1	1	1	0	1	0	0	0	...	1	1	1	1	1	0	1	1	1	1

5 rows × 41 columns

So our data matrix $X$ has 441 rows (legislators) and 41 columns (roll calls).

Each legislator is a point in 41-dimensional space.

X = vote_pivot.to_numpy()
X.shape

(441, 41)

A first attempt: plotting the raw columns

We can only look at two dimensions at a time on a screen. The obvious approach is to select two columns and produce a scatter plot.

px.scatter(vote_pivot, x='555', y='553',
           title='Two roll calls at a time', width=700, height=500)

That is four dots.

Every legislator sits at one of four corners, so 441 points collapse onto 4 visible marks. Even if we jittered them apart, there would be $\binom{41}{2} = 820$ such plots to inspect.

Selecting two of the original columns is not adequate. We require two new columns, constructed from all 41.

PCA in three lines

from sklearn.decomposition import PCA

model = PCA(n_components=2)
Z = model.fit_transform(vote_pivot)
Z.shape

(441, 2)

Each of the 441 legislators is now described by 2 numbers rather than 41. We plot them below.

px.scatter(x=Z[:, 0], y=Z[:, 1],
           labels={'x': 'z1', 'y': 'z2'},
           title='Vote data projected onto 2 dimensions',
           width=800, height=600, opacity=0.7)

Two clusters appear, although nothing in the input indicated that clusters should exist. The algorithm never saw a party label, only zeros and ones.

We now bring in the identity of each legislator, from unitedstates/congress-legislators.

# Static copy of the 2019 roster, so it matches our voting data.
legislators_data = yaml.safe_load(open('data/legislators-2019.yaml'))

def to_date(s):
    return datetime.strptime(s, '%Y-%m-%d')

legs = pd.DataFrame(
    columns=['leg_id', 'first', 'last', 'state', 'chamber', 'party', 'birthday'],
    data=[[x['id']['bioguide'],
           x['name']['first'],
           x['name']['last'],
           x['terms'][-1]['state'],
           x['terms'][-1]['type'],
           x['terms'][-1]['party'],
           to_date(x['bio']['birthday'])] for x in legislators_data])
legs['age'] = 2024 - legs['birthday'].dt.year
legs.head()

	leg_id	first	last	state	chamber	party	birthday	age
0	B000944	Sherrod	Brown	OH	sen	Democrat	1952-11-09	72
1	C000127	Maria	Cantwell	WA	sen	Democrat	1958-10-13	66
2	C000141	Benjamin	Cardin	MD	sen	Democrat	1943-10-05	81
3	C000174	Thomas	Carper	DE	sen	Democrat	1947-01-23	77
4	C001070	Robert	Casey	PA	sen	Democrat	1960-04-13	64

vote_2d = pd.DataFrame(Z, index=vote_pivot.index, columns=['z1', 'z2'])
vote_2d = vote_2d.join(legs.set_index('leg_id'))
vote_2d.head()

	z1	z2	first	last	state	chamber	party	birthday	age
member
A000055	-3.061356	0.364191	Robert	Aderholt	AL	rep	Republican	1965-07-22	59.0
A000367	-0.188870	-2.433565	Justin	Amash	MI	rep	Independent	1980-04-18	44.0
A000369	-2.844370	0.821619	Mark	Amodei	NV	rep	Republican	1958-06-12	66.0
A000370	2.607536	0.127977	Alma	Adams	NC	rep	Democrat	1946-05-27	78.0
A000371	2.607536	0.127977	Pete	Aguilar	CA	rep	Democrat	1979-06-19	45.0

px.scatter(vote_2d, x='z1', y='z2', color='party',
           title='Vote data, colored by party (PCA never saw this column)',
           width=800, height=600, opacity=0.7,
           color_discrete_map={'Democrat': 'blue', 'Republican': 'red', 'Independent': 'green'},
           hover_data=['first', 'last', 'state'],
           render_mode='svg')

The structure that PCA found corresponds to party affiliation.

There is substantial overplotting, since many legislators vote identically and therefore land on exactly the same point. We add jitter to reveal the density.

rng = np.random.default_rng(42)
vote_2d['z1_jittered'] = vote_2d['z1'] + rng.normal(0, 0.1, len(vote_2d))
vote_2d['z2_jittered'] = vote_2d['z2'] + rng.normal(0, 0.1, len(vote_2d))

px.scatter(vote_2d, x='z1_jittered', y='z2_jittered', color='party', size='age',
           title='Vote data (jittered)',
           width=800, height=600, opacity=0.7, size_max=10,
           color_discrete_map={'Democrat': 'blue', 'Republican': 'red', 'Independent': 'green'},
           hover_data=['first', 'last', 'state', 'party'])

How far does this go? If we use only the sign of the first coordinate, how often does it agree with party affiliation?

labeled = vote_2d.dropna(subset=['party'])
guess = np.where(labeled['z1'] > 0, 'Democrat', 'Republican')
(guess == labeled['party']).mean()

np.float64(0.979498861047836)

A single number per legislator, derived from all 41, recovers party affiliation for approximately 98% of the House.

How much information was discarded?

We replaced 41 columns with 2. The attribute explained_variance_ratio_ reports the fraction of the spread in the data accounted for by each new coordinate.

model.explained_variance_ratio_

array([0.80299948, 0.05260076])

model.explained_variance_ratio_.sum()

np.float64(0.8556002417105795)

The first coordinate alone accounts for roughly 80%. The full profile is shown below.

model10 = PCA(n_components=10).fit(vote_pivot)

px.line(y=model10.explained_variance_ratio_, markers=True,
        labels={'x': 'component', 'y': 'fraction of total spread'},
        title='Spread accounted for by each component',
        width=700, height=450)

There is a sharp decrease after the first component, followed by a long flat tail. This shape is what makes the two-dimensional plot trustworthy: no third direction accounts for an appreciable share of the spread.

Note that the data matrix is not actually low rank:

np.linalg.matrix_rank(X - X.mean(axis=0))

np.int64(41)

The rank is 41, which is full. This is therefore not a case of exactly redundant columns that may be deleted. The data is only approximately low dimensional, and that distinction is central to what follows.

What the model consists of

The call to fit estimated something. We inspect it below.

model.components_.shape

(2, 41)

model.components_

array([[ 0.02883302,  0.11337258,  0.18416809,  0.18362783,  0.00522399,
        -0.17220376,  0.1673402 , -0.1538741 , -0.16984021, -0.17817072,
         0.17498465, -0.17397851, -0.17174664, -0.17322967, -0.17549737,
         0.17788976, -0.17150954, -0.00193311,  0.182243  ,  0.18352844,
         0.00717233,  0.18372791,  0.18357099,  0.11563949, -0.1503738 ,
         0.17999325,  0.02464832,  0.18251365,  0.18306103,  0.10361375,
        -0.16405232,  0.18510959,  0.18517762,  0.01248441,  0.18324527,
         0.1827058 , -0.16030314,  0.18154702,  0.17405867,  0.00848885,
         0.17863309],
       [ 0.30363475,  0.24614996,  0.04782263,  0.05221178,  0.21887092,
         0.14505236,  0.05492475,  0.13180032,  0.14510786,  0.15696842,
         0.06032129,  0.15905892,  0.16137791,  0.15834598,  0.16409885,
         0.06411763,  0.1687306 ,  0.24195532,  0.08355767,  0.08140683,
         0.21845543,  0.07368227,  0.07218549,  0.24780437,  0.10844364,
         0.07043742,  0.2972522 ,  0.06551231,  0.0631842 ,  0.21127781,
         0.19274816,  0.06528183,  0.05498272,  0.21267071,  0.04526099,
         0.04844204,  0.20533786,  0.04531407,  0.07800742,  0.25967231,
         0.07677587]])

Two rows of 41 numbers. This is the entire model.

Each row is a set of weights over the 41 roll calls, and a legislator's new coordinates are the dot products of their voting record with these two rows:

w1 = model.components_[0]
Xc = X - X.mean(axis=0)          # PCA centers the data internally

# sklearn's answer for the first legislator, vs. a dot product we compute ourselves
print(Z[0, 0], Xc[0] @ w1)

-3.0613563028703172 -3.0613563028703172

px.bar(x=vote_pivot.columns, y=w1,
       labels={'x': 'roll call', 'y': 'weight in the first component'},
       title='The first row of components_',
       width=900, height=400)

# Center each column: subtract the House-wide Yes rate for that roll call.
vote_pivot_centered = vote_pivot - vote_pivot.mean()

# Per party, the average deviation from the House-wide Yes rate on each roll call.
party_yes_deviation = vote_pivot_centered.join(labeled['party']).groupby('party').mean()

party_yes_deviation_long = (party_yes_deviation
                            .reset_index()
                            .melt(id_vars='party',
                                  var_name='roll call',
                                  value_name='Yes Rate Centered'))

fig = px.bar(party_yes_deviation_long,
             x='roll call', y='Yes Rate Centered',
             facet_row='party', color='party',
             color_discrete_map={'Democrat': 'blue', 'Republican': 'red', 'Independent': 'green'},
             title='Party Yes rate relative to the House average, by roll call',
             width=900, height=800)
fig.for_each_annotation(lambda a: a.update(text=a.text.split('=')[-1]))  # 'party=Democrat' -> 'Democrat'
fig.update_layout(showlegend=False)  # the facet titles already name the party
fig

The method therefore reduces to a single question: how do we find the right $k$ rows of length $d$?

• How are those rows determined?
• In what sense are they the optimal choice?
• Why does projecting onto them preserve the structure of interest?

These are the subject of the remainder of the lecture.

Compressing images: PCA on Fashion-MNIST

The congressional votes were 41-dimensional. We now apply the same three lines to a dataset where each datapoint has 784 dimensions and, unlike a voting record, can be looked at directly. Fashion-MNIST is 60,000 grayscale 28 x 28 photographs of clothing in 10 categories.

# Fetch the Data
import torchvision
data = torchvision.datasets.FashionMNIST(root='data', train=True, download=True)

# Preprocess the data into numpy arrays
images = data.data.numpy().astype(float)
targets = data.targets.numpy() # integer encoding of class labels
class_dict = {i:class_name for i,class_name in enumerate(data.classes)}
labels = np.array([class_dict[t] for t in targets]) # raw class labels
n = len(images)

print("Loaded FashionMNIST dataset with {} samples.".format(n))
print("Classes: {}".format(class_dict))
print("Image shape: {}".format(images[0].shape))
print("Image dtype: {}".format(images[0].dtype))
print("Image 0:\n", images[0])

Loaded FashionMNIST dataset with 60000 samples.
Classes: {0: 'T-shirt/top', 1: 'Trouser', 2: 'Pullover', 3: 'Dress', 4: 'Coat', 5: 'Sandal', 6: 'Shirt', 7: 'Sneaker', 8: 'Bag', 9: 'Ankle boot'}
Image shape: (28, 28)
Image dtype: float64
Image 0:
 [[  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
    0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
    0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
    0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   0.
    0.  13.  73.   0.   0.   1.   4.   0.   0.   0.   0.   1.   1.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   3.   0.
   36. 136. 127.  62.  54.   0.   0.   0.   1.   3.   4.   0.   0.   3.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   6.   0.
  102. 204. 176. 134. 144. 123.  23.   0.   0.   0.   0.  12.  10.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
  155. 236. 207. 178. 107. 156. 161. 109.  64.  23.  77. 130.  72.  15.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   0.  69.
  207. 223. 218. 216. 216. 163. 127. 121. 122. 146. 141.  88. 172.  66.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.   1.   0. 200.
  232. 232. 233. 229. 223. 223. 215. 213. 164. 127. 123. 196. 229.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0. 183.
  225. 216. 223. 228. 235. 227. 224. 222. 224. 221. 223. 245. 173.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0. 193.
  228. 218. 213. 198. 180. 212. 210. 211. 213. 223. 220. 243. 202.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   3.   0.  12. 219.
  220. 212. 218. 192. 169. 227. 208. 218. 224. 212. 226. 197. 209.  52.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   6.   0.  99. 244.
  222. 220. 218. 203. 198. 221. 215. 213. 222. 220. 245. 119. 167.  56.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   4.   0.   0.  55. 236.
  228. 230. 228. 240. 232. 213. 218. 223. 234. 217. 217. 209.  92.   0.]
 [  0.   0.   1.   4.   6.   7.   2.   0.   0.   0.   0.   0. 237. 226.
  217. 223. 222. 219. 222. 221. 216. 223. 229. 215. 218. 255.  77.   0.]
 [  0.   3.   0.   0.   0.   0.   0.   0.   0.  62. 145. 204. 228. 207.
  213. 221. 218. 208. 211. 218. 224. 223. 219. 215. 224. 244. 159.   0.]
 [  0.   0.   0.   0.  18.  44.  82. 107. 189. 228. 220. 222. 217. 226.
  200. 205. 211. 230. 224. 234. 176. 188. 250. 248. 233. 238. 215.   0.]
 [  0.  57. 187. 208. 224. 221. 224. 208. 204. 214. 208. 209. 200. 159.
  245. 193. 206. 223. 255. 255. 221. 234. 221. 211. 220. 232. 246.   0.]
 [  3. 202. 228. 224. 221. 211. 211. 214. 205. 205. 205. 220. 240.  80.
  150. 255. 229. 221. 188. 154. 191. 210. 204. 209. 222. 228. 225.   0.]
 [ 98. 233. 198. 210. 222. 229. 229. 234. 249. 220. 194. 215. 217. 241.
   65.  73. 106. 117. 168. 219. 221. 215. 217. 223. 223. 224. 229.  29.]
 [ 75. 204. 212. 204. 193. 205. 211. 225. 216. 185. 197. 206. 198. 213.
  240. 195. 227. 245. 239. 223. 218. 212. 209. 222. 220. 221. 230.  67.]
 [ 48. 203. 183. 194. 213. 197. 185. 190. 194. 192. 202. 214. 219. 221.
  220. 236. 225. 216. 199. 206. 186. 181. 177. 172. 181. 205. 206. 115.]
 [  0. 122. 219. 193. 179. 171. 183. 196. 204. 210. 213. 207. 211. 210.
  200. 196. 194. 191. 195. 191. 198. 192. 176. 156. 167. 177. 210.  92.]
 [  0.   0.  74. 189. 212. 191. 175. 172. 175. 181. 185. 188. 189. 188.
  193. 198. 204. 209. 210. 210. 211. 188. 188. 194. 192. 216. 170.   0.]
 [  2.   0.   0.   0.  66. 200. 222. 237. 239. 242. 246. 243. 244. 221.
  220. 193. 191. 179. 182. 182. 181. 176. 166. 168.  99.  58.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.  40.  61.  44.  72.  41.  35.   0.
    0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
    0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
    0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.]]

px.imshow(images[0], color_continuous_scale='gray_r') 

def show_images(images, max_images=40, ncols=5, labels = None):
    """Visualize a subset of images from the dataset.
    Args:
        images (np.ndarray): Array of images to visualize [img,row,col].
        max_images (int): Maximum number of images to display.
        ncols (int): Number of columns in the grid.
        labels (np.ndarray, optional): Labels for the images, used for facet titles.
    Returns:
        plotly.graph_objects.Figure: A Plotly figure object containing the images.
    """
    n = min(images.shape[0], max_images) # number of images to show
    px_height = 220 # height of each image in pixels
    fig = px.imshow(images[:n, :, :], color_continuous_scale='gray_r', 
                    facet_col = 0, facet_col_wrap=ncols,
                    height = px_height * int(np.ceil(n/ncols)))
    fig.update_layout(coloraxis_showscale=False)
    if labels is not None:
        # Extract the facet number and replace with the label.
        fig.for_each_annotation(lambda a: a.update(text=labels[int(a.text.split("=")[-1])]))
    return fig

show_images(images, 20, labels=labels)

Each image is a datapoint in 784-dimensional space

The voting data gave us one row per legislator and one column per roll call. We do the same thing here: one row per image, one column per pixel. A 28 x 28 image becomes a single row of $28 \times 28 = 784$ numbers.

X_img = images.reshape(n, -1)   # (60000, 28, 28) -> (60000, 784)
X_img.shape

(60000, 784)

PCA centers the data before it does anything else, so the first thing it computes is the mean of those 60,000 rows. Reshaped back to 28 x 28, the mean is itself an image.

mean_image = X_img.mean(axis=0)

px.imshow(mean_image.reshape(28, 28), color_continuous_scale='gray_r',
          title='The average of all 60,000 images', width=400, height=400)

The principal components are also images

For the votes, components_ was a $k \times 41$ matrix, and each row was a set of weights over the 41 roll calls. Here it is a $k \times 784$ matrix, and each row is a set of weights over the 784 pixels. A row of 784 numbers can be reshaped into a 28 x 28 picture, so we can look directly at the model.

We fit 200 components once, and use the leading $k$ of them below. Because the components are nested, components_[:k] is exactly what PCA(n_components=k) would have produced.

pca_img = PCA(n_components=200).fit(X_img)
pca_img.components_.shape

(200, 784)

n_show = 10
comps = pca_img.components_[:n_show].reshape(n_show, 28, 28)

fig = px.imshow(comps, facet_col=0, facet_col_wrap=5,
                color_continuous_scale='RdBu_r', color_continuous_midpoint=0,
                height=440, title='The first 10 principal components, viewed as images')
fig.for_each_annotation(lambda a: a.update(text=f"PC {int(a.text.split('=')[-1]) + 1}"))
fig.update_layout(coloraxis_showscale=False)
fig

Red is a positive weight and blue is a negative one. These are not garments; they are contrasts. PC 1 separates wide dark regions from narrow ones, which is roughly the distinction between a shirt and a shoe. Later components encode sleeves, straps, and the gap between trouser legs. The overall sign of each component is arbitrary.

How many components do we need?

cumvar = np.cumsum(pca_img.explained_variance_ratio_)

px.line(x=np.arange(1, 201), y=cumvar,
        labels={'x': 'number of components k', 'y': 'cumulative fraction of spread'},
        title='Spread captured by the first k components', width=750, height=450)

The curve rises steeply and then flattens: 50 of the 784 directions account for 86% of the spread, and 200 account for 95%. Compare this to the votes, where a single component captured 80%. Images are low dimensional, but not nearly as aggressively so.

Reconstruction

Compression is only useful if we can get the image back. Keeping $k$ scores and then undoing the projection gives

$$\hat{x} = \bar{x} + \sum_{j=1}^{k} z_j w_j$$

a picture rebuilt as the mean image plus a weighted sum of $k$ component images.

Z_img = pca_img.transform(X_img)   # (60000, 200) scores

def reconstruct(k, rows):
    """Rebuild images from only their first k scores."""
    return Z_img[rows, :k] @ pca_img.components_[:k] + mean_image

ks = [1, 2, 5, 10, 25, 50, 100, 200]
i = 0   # the ankle boot from the top of this section

ladder = np.vstack([X_img[i]] + [reconstruct(k, [i]) for k in ks])

# Reconstructions can fall slightly outside [0, 255], so we clip them for display.
show_images(np.clip(ladder, 0, 255).reshape(-1, 28, 28), ncols=3,
            labels=['original (784)'] + [f'k = {k}' for k in ks])

One number produces a dark blob. Ten produce something identifiable as a boot. By 50 the silhouette and the shading are right, and the remaining 734 dimensions mostly carry texture and noise.

Below we do the same at $k = 50$ for eight random images. The top row is the original, the bottom row is 50 numbers.

rng_img = np.random.default_rng(189)
rows = rng_img.choice(n, 8, replace=False)
k = 50

side_by_side = np.vstack([X_img[rows], reconstruct(k, rows)])
show_images(np.clip(side_by_side, 0, 255).reshape(-1, 28, 28), ncols=8,
            labels=[labels[r] for r in rows] + [f'k = {k}' for _ in rows])

What did that actually save?

To store the whole dataset we need the $n \times k$ table of scores, plus the basis we need in order to decode it: the $k \times 784$ components and the 784-pixel mean. The basis is paid for once, no matter how many images we compress.

d = X_img.shape[1]

for k in [10, 50, 100]:
    stored = n * k + k * d + d
    rmse = np.sqrt(((reconstruct(k, slice(None)) - X_img) ** 2).mean())
    print(f"k = {k:3d} | scores {n*k:>9,} + basis {k*d + d:>7,} = {stored:>9,} numbers "
          f"| {n*d/stored:5.1f}x smaller | RMSE {rmse:5.1f}")

k =  10 | scores   600,000 + basis   8,624 =   608,624 numbers |  77.3x smaller | RMSE  39.8
k =  50 | scores 3,000,000 + basis  39,984 = 3,039,984 numbers |  15.5x smaller | RMSE  27.9

k = 100 | scores 6,000,000 + basis  79,184 = 6,079,184 numbers |   7.7x smaller | RMSE  22.3

At $k = 50$ the dataset is 15x smaller and the images survive. The basis is only 1.3% of the stored bytes, so essentially all of the cost is the 50 numbers per image.

This is lossy compression of the same general kind as JPEG. The difference is that JPEG uses a fixed, universal basis (cosines), while PCA learns a basis from this particular collection of images. That is why the components above look like clothing contrasts rather than generic ripples, and it is also why the basis only compresses images that resemble the training set.

Can we run the decoder backwards to invent new clothes?

The reconstruction step $\hat{x} = \bar{x} + \sum_j z_j w_j$ turns 50 numbers into an image, and it does not care where those numbers came from. So here is a tempting idea: instead of taking $z$ from a real image, make $z$ up, and see what comes out.

For this to work, the made-up $z$ has to look like the $z$ of a real image. So first we look at how the real scores are distributed.

sub = np.random.default_rng(189).choice(n, 4000, replace=False)
scores_2d = pd.DataFrame({'z1': Z_img[sub, 0], 'z2': Z_img[sub, 1], 'class': labels[sub]})

px.scatter(scores_2d, x='z1', y='z2', color='class', opacity=0.6,
           title='The first two scores of 4,000 images', width=850, height=600)

This is not one cloud. Footwear sits on the left, tops sit on the upper right, trousers hang below, bags sit on top. The score distribution is multimodal, and there are wide empty regions between the groups.

Let us ignore that for a moment and do the simplest thing: fit a single Gaussian to the 50 scores, draw from it, and decode.

k = 50
Zk = Z_img[:, :k]

def sample_images(Z_ref, n_samples, rng):
    """Fit one Gaussian to the scores in Z_ref, draw from it, and decode into images."""
    draws = rng.multivariate_normal(Z_ref.mean(axis=0), np.cov(Z_ref.T), n_samples)
    return draws @ pca_img.components_[:Z_ref.shape[1]] + mean_image

fake = sample_images(Zk, 16, np.random.default_rng(0))
show_images(np.clip(fake, 0, 255).reshape(-1, 28, 28), ncols=8, max_images=16)

These are garment-shaped smudges. Several are two items at once: a sleeve fading into a trouser leg, a shoe ghosted over a shirt. They have the statistics of the dataset without being plausible members of it.

The scatter plot above explains why. A single Gaussian is one blob, so most of its mass lands in the empty space between the clusters, and a point halfway between "sneaker" and "pullover" decodes to a superposition of the two. We can check that the samples really are landing in unoccupied territory by measuring how far each one is from the nearest real image.

rng_nn = np.random.default_rng(1)
reference = Zk[rng_nn.choice(n, 8000, replace=False)]
real_pts  = Zk[rng_nn.choice(n, 200, replace=False)]
fake_pts  = rng_nn.multivariate_normal(Zk.mean(axis=0), np.cov(Zk.T), 200)

def nn_distance(query, reference):
    """Distance from each query point to its closest neighbour in reference."""
    sq = ((query[:, None, :] - reference[None, :, :]) ** 2).sum(axis=2)
    return np.sqrt(sq.min(axis=1))

print(f"real image  -> nearest real image: {np.median(nn_distance(real_pts, reference)):.0f}")
print(f"fake sample -> nearest real image: {np.median(nn_distance(fake_pts, reference)):.0f}")

real image  -> nearest real image: 584
fake sample -> nearest real image: 1292

A real image has a real neighbour roughly twice as close. The samples are not near the data; they are in the gaps.

The fix follows directly from the diagnosis. The problem was fitting one blob to ten clusters, so we fit one Gaussian per class instead, still in the same 50-dimensional score space, and still decoding with the same components.

rng_cls = np.random.default_rng(0)
chosen = [0, 7, 8, 1]   # T-shirt/top, Sneaker, Bag, Trouser

per_class = np.vstack([sample_images(Zk[targets == c], 4, rng_cls) for c in chosen])

show_images(np.clip(per_class, 0, 255).reshape(-1, 28, 28), ncols=4, max_images=16,
            labels=[class_dict[c] for c in chosen for _ in range(4)])

These work. Each row is four garments that do not exist in the dataset, and they are recognizably sneakers, t-shirts, bags, and trousers, with varied heights, widths, and shading. They are blurry, because 50 components cannot represent a sharp edge and because a Gaussian is still only an approximation of a class, but they are plausible items rather than superpositions.

So the honest answer is: the decoder is fine, and the hard part is knowing which $z$ to feed it. Sampling works exactly as well as our model of the score distribution does.

One caution about what the subspace can do

It is tempting to read the score space as a space of concepts, where moving from one image to another should morph a shirt into a shoe. It cannot, and the reason is that the map from $z$ to pixels is linear. Interpolating between two images in score space is algebraically identical to cross-fading the two images in pixel space.

a = np.where(targets == 0)[0][0]   # a t-shirt
b = np.where(targets == 7)[0][0]   # a sneaker

t = np.linspace(0, 1, 8)[:, None]
path = (1 - t) * Zk[a] + t * Zk[b]
blend = path @ pca_img.components_[:k] + mean_image

show_images(np.clip(blend, 0, 255).reshape(-1, 28, 28), ncols=8, max_images=8)

# The same path, computed instead by blending pixels and then projecting. Identical.
pixel_blend = (1 - t) * X_img[a] + t * X_img[b]
projected = (pixel_blend - mean_image) @ pca_img.components_[:k].T @ pca_img.components_[:k] + mean_image

np.abs(blend - projected).max()

np.float64(4.547473508864641e-13)

The shirt does not become a shoe; it dissolves while a shoe appears underneath. Linearity is what makes PCA cheap to fit, easy to interpret, and provably optimal in the sense we are about to define, and it is also precisely what stops it from being a generative model of images. Getting a genuine morph requires a decoder that is not restricted to a linear subspace.

Summary of this section

• Each image is a point in $\mathbb{R}^{784}$; PCA finds a $k$-dimensional subspace it nearly lies in.
• The components are pictures, and reconstruction is the mean image plus a weighted sum of them.
• $k = 50$ compresses the dataset 15x with the content of the images intact.
• Decoding invented scores does generate new clothing, but only once the score distribution is modeled per class. A single Gaussian samples the empty space between clusters.
• The subspace is linear, so interpolation is a cross-fade, not a morph.