Eigenvectors, and What PCA Is Doing

Context

An eigenvector of a square matrix is a direction the matrix only stretches: multiplied by the matrix, it comes back as itself times one number, its eigenvalue. PCA (principal component analysis) turns a table's centered columns onto the directions they spread along most, which for standardized columns are the eigenvectors of their correlation matrix, so the first new column runs along the direction the listings spread most and each one after it along the most of what is left.

A pricing team uses it to draw 3,000 listings on one chart whose two axes keep three quarters of six columns' spread, to see which columns rise and fall as one before a model is handed twins, and to give a model a few columns that do not move together in place of many that do. The part of a listing the first few directions cannot rebuild points at listings whose columns break the pattern the rest share.

The failures print plausible numbers. On unscaled columns PCA returns the price alone as its first direction, with a share of 1.0 and no warning. np.linalg.eigh lists the smallest eigenvalue first, so code that reads its first column takes the weakest direction, and NumPy and scikit-learn can print one direction with opposite signs, so a report that reads meaning into a sign contradicts itself. Lesson 3.8.1 grouped the listings by distance, Lesson 2.3.4 projected one vector onto another, and this lesson projects every listing onto the directions its table spreads along most.

Goal

• Verify that a vector is an eigenvector of a matrix by multiplying, and read its eigenvalue off the product.
• Build the correlation matrix of standardized columns as Z.T @ Z / len(Z) and check it against pandas' corr of the filled columns.
• Compute eigenvalues and eigenvectors with np.linalg.eigh, order them largest first, and check every pair at once.
• Match PCA's components_, shares, explained variances and transform to eigh and to 2.3.4's projection, signs and divisors included.
• Rebuild the standardized table from fewer components, show that its error is the share left out, and explain precisely why unscaled columns leave the price alone.

Prerequisites

Neither is required, and the lesson works out every number in front of you. They are here if you want eigenvectors drawn, or PCA worked by hand, before you start.

• YouTube - "Eigenvectors and eigenvalues | Chapter 14, Essence of linear algebra" by 3Blue1Brown. A matrix drawn as a transformation of the plane, and the vectors it leaves on their own line, stretched or squished.
• YouTube - "StatQuest: Principal Component Analysis (PCA), Step-by-Step" by StatQuest with Josh Starmer. PCA on two genes measured in six mice from 3:23, the second component at a right angle to the first at 12:56, and each component's percent of the variation at 15:03.

Table of Contents

• Directions a Matrix Only Stretches
• The Correlation Matrix of Standardized Columns
• np.linalg.eigh
• PCA as a Rotation
• Keeping Fewer Components

Theory

Nothing to install: np.linalg.eigh is part of NumPy, and PCA, from sklearn.decomposition, arrived with Lesson 3.1.1's scikit-learn. Every Theory block below goes into one script, pca.py, saved in ~/dev/bootcamp beside listings.csv and run with python3 pca.py, each block added under the one before it.

Directions a Matrix Only Stretches

What. Multiplying a vector by a square matrix, Lesson 2.3.2's @, generally turns it and changes its length. An eigenvector of the matrix is a direction the matrix only stretches: A @ v equals lam * v for one number lam, the eigenvalue, which is the stretch, and a flip as well when it is negative. Why. Along its eigenvectors a matrix only scales, so a vector split into pieces along them goes through the matrix as each piece times its own eigenvalue. For the next subtitle's correlation matrix, those directions are the ones the listings spread along. When. Whenever a square matrix is taken apart into the directions it stretches, the correlation matrix PCA works on above all, and whenever a claimed eigenvector is checked, with one multiplication.

import numpy as np

A = np.array([[2, 1], [1, 2]])
print(A @ np.array([1, 0]))   # [2 1]
print(A @ np.array([1, 1]))   # [3 3]
print(A @ np.array([1, -1]))  # [ 1 -1]

[1, 0] came back as [2, 1], longer and turned, so it is no eigenvector. [1, 1] came back three times itself, an eigenvector with eigenvalue 3, and [1, -1] came back as itself, an eigenvector with eigenvalue 1. The test is that every entry of A @ v is the same multiple of the matching entry of v. Any multiple of an eigenvector, zero aside, is one too:

print(A @ np.array([2, 2]))    # [6 6]
print(A @ np.array([-1, -1]))  # [-3 -3]
print(A @ np.array([-1, 1]))   # [-1  1]

An eigenvector is really a line through the center, and that line holds two vectors of length 1, one pointing each way. A function that returns one has to pick a sign, and nothing in the matrix decides which. Split a vector into pieces along the eigenvectors, and the matrix stretches each piece by its own eigenvalue:

up = np.array([1, 1])                 # eigenvalue 3
down = np.array([1, -1])              # eigenvalue 1
print(0.5 * up + 0.5 * down)          # [1. 0.]
print(0.5 * 3 * up + 0.5 * 1 * down)  # [2. 1.]
print(up @ down)                      # 0
print(np.array_equal(A, A.T))         # True

[1, 0] is half of up plus half of down, and A tripled the first half and kept the second, which adds up to the first block's [2, 1]. up @ down is 0, Lesson 2.3.1's orthogonality, and that is no accident: A equals its transpose, which makes it a symmetric matrix (a square matrix whose entry [i, j] equals its entry [j, i]), and the eigenvectors of a symmetric matrix can always be chosen at right angles to each other.

DON'T:

• DON'T call a vector an eigenvector because the matrix lengthened it. [1, 0] came back as [2, 1], longer and turned, and only a vector that stays on its own line is one.
• DON'T expect one eigenvector per eigenvalue. [1, 1], [2, 2] and [-1, -1] are all eigenvectors of A for 3, so even at length 1 there are two, pointing opposite ways.

The Correlation Matrix of Standardized Columns

What. On standardized columns, whose means are 0 and variances 1, the covariance matrix (the square table holding, for every pair of columns, the average product of their deviations from their means, each column's variance on the diagonal) is Z.T @ Z / len(Z), and it is Lesson 2.7.1's correlation matrix. Why. It is the matrix PCA takes apart. Entry [i, j] says how closely columns i and j rise and fall together across the listings, so the directions it stretches most are the directions the listings spread along most. When. Before PCA, which works on it whether you build it or not, and whenever every pair's correlation is wanted in one array.

The table is the one Lesson 3.8.1 clustered: with no target to check held-out rows against, none are held out, so the imputer and the scaler are fitted on all 3,000 listings.

import pandas as pd
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
cols = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km", "price"]
prep = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
Z = prep.fit_transform(listings[cols])
print(Z.shape)                       # (3000, 6)
print(round(prep[-1].scale_[0], 2))  # 43.18

prep[-1] is the fitted scaler, and its first scale_ is one standard deviation of area, 43.18 square meters, the unit Z's first column counts in. The covariance matrix, and pandas' correlation matrix beside it:

C = Z.T @ Z / len(Z)
print(C.shape)                  # (6, 6)
print(np.round(np.diag(C), 6))  # [1. 1. 1. 1. 1. 1.]
filled = pd.DataFrame(prep[0].transform(listings[cols]), columns=cols)
corr = filled.corr()
print(np.allclose(C, corr))     # True
print(corr.round(3))

             area_sqm  rooms  bathrooms  year_built  distance_km  price
area_sqm        1.000  0.900      0.761       0.073        0.277  0.725
rooms           0.900  1.000      0.835       0.057        0.265  0.681
bathrooms       0.761  0.835      1.000       0.048        0.230  0.569
year_built      0.073  0.057      0.048       1.000        0.289  0.162
distance_km     0.277  0.265      0.230       0.289        1.000  0.019
price           0.725  0.681      0.569       0.162        0.019  1.000

Entry [i, j] of Z.T @ Z is column i dotted with column j, Lesson 2.3.2's feature-by-feature table, so C holds the average product of every pair of standardized columns. prep[0] is the fitted imputer, whose transform returns the filled columns unscaled, and np.diag, handed a square matrix instead of Lesson 2.3.2's 1-D array, returns its diagonal.

The divisor is n because of how Z was made. Lesson 2.7.1's Pearson r was the dot product of two columns standardized with pandas' ddof=1, divided by n minus 1, and the scaler standardizes with NumPy's ddof=0, Lesson 3.2.1's finding, so here the divisor is n: the one inside the standardization and the one outside must agree, and then they cancel.

# False: the scaler divided by n, not n - 1
print(np.allclose(Z.T @ Z / (len(Z) - 1), corr))
print(np.allclose(C, C.T))  # True

Dividing by n minus 1 overshoots the correlation matrix by 3,000 over 2,999. C is symmetric, as every covariance matrix is, because column i dotted with column j is column j dotted with column i. Read row by row, area and rooms correlate at 0.9 and rooms and bathrooms at 0.835, Lesson 3.3.2's twins, and the price rises with all three, at 0.569 to 0.725. Build year and distance correlate with each other at 0.289 and with the other four at 0.019 to 0.277: four columns rise and fall together, and two mostly go their own way.

DON'T:

• DON'T check C against listings[cols].corr(). It keeps each pair's own complete rows, Lesson 2.7.1's rule, so area and rooms come out at 0.924 without the 125 filled areas against the filled table's 0.900, and np.allclose prints False.
• DON'T divide the scaler's output by len(Z) - 1. The scaler divided by n, and the matrix then misses the correlation by 3,000 over 2,999.

np.linalg.eigh

What. np.linalg.eigh(C) returns the eigenvalues and eigenvectors of a symmetric matrix as a pair: the eigenvalues ascending in a 1-D array, and the eigenvectors as the columns of a 2-D array, each of length 1, column j belonging to eigenvalue j. Why. It is built for symmetric matrices: the eigenvalues come back real and sorted and the eigenvectors at right angles, everything PCA needs. When. Whenever the matrix is symmetric, a covariance or correlation matrix above all, and never otherwise.

vals, vecs = np.linalg.eigh(C)
print(np.round(vals, 3))     # [0.086 0.191 0.365 0.806 1.214 3.338]
print(round(vals.sum(), 6))  # 6.0
print(vecs.shape)            # (6, 6)

A matrix's eigenvalues always add up to its diagonal's sum, and C's diagonal is six variances of 1, so 6 is the total spread of the six standardized columns, which the eigenvalues divide among six perpendicular directions. The largest takes 3.338 of it.

top = vecs[:, -1]                          # the column of 3.338
print(np.round(C @ top / top, 4))
print(np.allclose(C @ vecs, vecs * vals))  # True
print(np.round(np.linalg.norm(vecs, axis=0), 6))
# True: every column at a right angle to the rest
print(np.allclose(vecs.T @ vecs, np.eye(6)))

[3.3378 3.3378 3.3378 3.3378 3.3378 3.3378]
[1. 1. 1. 1. 1. 1.]

Every entry of C @ top is 3.3378 times the matching entry of top, the first subtitle's test on six entries. vecs * vals multiplies column j by vals[j], Lesson 2.2.4's broadcasting, so one np.allclose checks all six pairs, and the norms along axis=0, one per column, show every eigenvector at length 1.

top_two = pd.DataFrame(vecs[:, [-1, -2]], index=cols,
                       columns=["largest", "second"])
print(top_two.round(3))

             largest  second
area_sqm      -0.516  -0.068
rooms         -0.520  -0.086
bathrooms     -0.480  -0.090
year_built    -0.091   0.719
distance_km   -0.182   0.664
price         -0.436  -0.145

The largest direction weighs area, rooms, bathrooms and price at -0.436 to -0.52, all one sign, and build year and distance far less: the four columns that rose and fell together in C do so along it. The second, 1.214, weighs build year at 0.719 and distance at 0.664. Each weight multiplies a column counted in its own standard deviations, never in square meters or years.

eigh promises ascending eigenvalues and eigenvectors of length 1, and nothing about which way each one points, so on another machine a column can print with every sign flipped. Every sign eigh prints in this lesson can differ on yours, and no eigenvalue and no line can. eigh also takes the symmetry on trust:

lopsided = np.array([[2.0, 5.0], [1.0, 2.0]])
lop_vals, lop_vecs = np.linalg.eigh(lopsided)
print(np.round(lop_vals, 6))          # [1. 3.]
v = lop_vecs[:, 1]                    # its vector for 3
print(np.round(lopsided @ v / v, 6))  # [7. 3.]

eigh read only the lower triangle, the entries on and below the diagonal, so it never saw the 5 and answered for [[2, 1], [1, 2]], the first subtitle's A. Its vector for 3 is stretched 7 times in one entry and 3 times in the other, no eigenvector of lopsided at all, and nothing warned.

DON'T:

• DON'T take vecs[:, 0] as the main direction. eigh sorts ascending, so column 0 belongs to 0.086, the smallest eigenvalue, and the largest is vecs[:, -1].
• DON'T read an eigenvector from a row. vecs[-1] holds the price's weight in each of the six eigenvectors, and C @ vecs[-1] / vecs[-1] came out here as six different numbers, from -6.203 to 1.218.
• DON'T hand eigh a matrix that is not symmetric. It read lopsided's lower triangle and returned 1 and 3 without a warning.
• DON'T test an eigenvector with a 0 in it by dividing. For D = np.diag([2.0, 3.0]) and its eigenvector v = np.array([1.0, 0.0]), D @ v / v prints [ 2. nan] under RuntimeWarning: invalid value encountered in divide, while np.allclose(D @ v, 2.0 * v) prints True.

PCA as a Rotation

What. PCA (principal component analysis) turns centered columns onto the directions they spread along most, which for standardized columns such as Z are the eigenvectors of their correlation matrix, largest eigenvalue first. PCA().fit(Z) keeps those eigenvectors as the rows of components_, the principal components (the new axes, each a direction of length 1 across the six columns), and transform returns each listing's scores, its coordinates along them. Perpendicular axes of length 1 move no listing nearer to the center or farther from it: the turn is a rotation. Why. In the new axes the scores no longer move together, the first carries as much of the spread as any direction can and each next one as much of what is left, so the first few describe a listing with the least lost. When. Whenever columns that move together are to be described by fewer directions: a chart of thousands of listings on two axes, a look at which columns rise as one, or a few inputs for a model.

from sklearn.decomposition import PCA

pca = PCA().fit(Z)
first_two = pd.DataFrame(pca.components_[:2], columns=cols)
print(first_two.round(3))

   area_sqm  rooms  bathrooms  year_built  distance_km  price
0     0.516  0.520       0.48       0.091        0.182  0.436
1    -0.068 -0.086      -0.09       0.719        0.664 -0.145

Row 0 is eigh's largest column with every sign flipped here, and row 1 its second as it was. All six at once:

order = vecs[:, ::-1].T                     # largest first, as rows
print(np.allclose(pca.components_, order))  # False
# True: the same lines once the signs are dropped
print(np.allclose(np.abs(pca.components_), np.abs(order)))

The first check compares the signs eigh chose and prints False here, and the second compares only the lines and holds on any machine. What PCA does on this table is what this lesson did: at svd_solver="auto", a table with fewer than 1,000 columns and more than ten times as many rows gets the "covariance_eigh" solver, which builds the covariance matrix, calls eigh on it and reverses the order. It then flips each component so that its weight largest in size is positive, 0.52 on rooms in row 0, scikit-learn's way of printing the same answer on every run and a fact about the library, not the listings.

ratio = pca.explained_variance_ratio_
print(np.round(ratio, 3))     # [0.556 0.202 0.134 0.061 0.032 0.014]
# True: each eigenvalue over their sum
print(np.allclose(ratio, vals[::-1] / vals.sum()))
print(np.round(pca.explained_variance_ / vals[::-1], 8))
print(len(Z) / (len(Z) - 1))  # 1.0003334444814937

[1.00033344 1.00033344 1.00033344 1.00033344 1.00033344 1.00033344]

The first component carries 0.556 of the six columns' spread and the second 0.202. explained_variance_ holds each component's explained variance (the variance of the listings' scores along it), and it divides by n minus 1 where C divided by n, so each value is its eigenvalue times 3,000 over 2,999. A share divides variances that all carry that factor, so it cancels.

scores = pca.transform(Z)
print(scores.shape)            # (3000, 6)
print(np.round(pca.mean_, 6))  # [ 0. -0. -0.  0.  0.  0.]
# True: one product scores every listing
print(np.allclose(scores, Z @ pca.components_.T))

transform subtracts mean_, the column means the fit learned, 0 on Z up to the float noise printed as -0., and projects what is left onto each component. A component has length 1, so each score is Lesson 2.3.4's scalar projection, a . b / |b| with |b| equal to 1, how far the listing reaches along the component.

comps = pca.components_
# True: perpendicular, each of length 1
print(np.allclose(comps @ comps.T, np.eye(6)))
before = np.linalg.norm(Z, axis=1)      # each listing's distance
after = np.linalg.norm(scores, axis=1)  # from the center
print(np.allclose(before, after))       # True
spread = scores.T @ scores / len(Z)
# True: the eigenvalues, and zeros off the diagonal
print(np.allclose(spread, np.diag(vals[::-1])))

Every listing keeps its distance from the center, Lesson 2.3.3's norm per row, which is what makes PCA a rotation, and flipping a component's sign mirrors that one axis and keeps every distance too. The scores' covariance matrix holds the eigenvalues on its diagonal and zeros elsewhere: on the listings PCA was fitted on, no two columns of scores move together. That no direction spreads them further than the first component can be put to 1,000 random directions:

rng = np.random.default_rng(42)
dirs = rng.normal(size=(1000, 6))            # 1,000 random directions
dirs = dirs / np.linalg.norm(dirs, axis=1)[:, np.newaxis]
along = Z @ dirs.T                           # each listing on each one
print(round((along ** 2).mean(axis=0).max(), 3))
print(round((scores[:, 0] ** 2).mean(), 3))  # 3.338

3.111

Each row of dirs is divided by its length, np.newaxis (Lesson 2.2.4) lining the 1,000 lengths up with the rows, and because Z is centered, a direction's mean squared score is the spread along it. The best random direction reaches 3.111, the first component 3.338. Six dimensions cannot be drawn and two can, so on area and price alone:

two = Z[:, [0, 5]]  # standardized area, price
flat = PCA().fit(two)
# the first subtitle's [1, 1] and [1, -1], at length 1
diag = np.array([1, 1]) / np.sqrt(2)
across = np.array([1, -1]) / np.sqrt(2)
lines = np.array([diag, across])
print(np.round(lines, 3))
# True: PCA found the two diagonals
print(np.allclose(np.abs(flat.components_), np.abs(lines)))
r = C[0, 5]         # area with price
print(round(r, 3))  # 0.725
print(round((1 + r) / 2, 3), round((1 - r) / 2, 3))
print(np.round(flat.explained_variance_ratio_, 3))
print(np.round(np.sqrt(flat.explained_variance_), 3))

[[ 0.707  0.707]
 [ 0.707 -0.707]]
0.863 0.137
[0.863 0.137]
[1.314 0.524]

Two standardized columns have the correlation matrix [[1, r], [r, 1]], which sends [1, 1] to 1 + r times itself and [1, -1] to 1 - r times itself, whatever r is, the first subtitle's A with 1 + r and 1 - r for its 3 and 1. So the components are the two diagonals at length 1, 0.707 being 1 over the square root of 2, and on two columns that rise together the rising diagonal carries (1 + r) / 2 of the spread. np.sqrt turns the explained variances, n minus 1 and all, into standard deviations: 1.314 along the first component and 0.524 across it.

# [134.0, 722000.0]
print(listings[["area_sqm", "price"]].iloc[16].tolist())
print(np.round(two[16], 3))        # [0.82  1.748]
first = two[16] @ diag
print(round(first, 3))             # 1.815
print(round(two[16] @ across, 3))  # -0.656
print(np.round(first * diag, 3))   # [1.284 1.284]

The listing at position 16, a townhouse of 134.0 square meters at 722,000, sits above the rising diagonal. Its first score, 1.815, is its projection onto the diagonal, and first * diag is where that projection lands, the foot of a perpendicular from the listing, Lesson 2.3.4's right angle. Its second score, -0.656, is its distance from the diagonal, negative because it sits on the side away from where [0.707, -0.707] points.

DON'T:

• DON'T read meaning into a component's sign. eigh returned the first direction here with every weight negative and PCA with every weight positive, one line through the center.
• DON'T line components_ up against vecs as eigh returned them. components_ holds rows, largest first, and vecs[:, ::-1].T is the form that matches, up to sign.
• DON'T read a component as one of the original columns, or a weight as a rate in listing units. The first component weighs four columns at 0.436 to 0.52 at once, and its 0.516 on area raises the first score by 0.516 when area rises by one standard deviation, 43.18 square meters, with the other five held fixed.

Keeping Fewer Components

What. PCA(n_components=2) keeps the first two components: transform returns two scores per listing, and inverse_transform rebuilds the six standardized columns from them, scores @ components_ + mean_. On the scaler's columns, the rebuild's mean squared error is the share of the spread the left-out components carry. Why. Two numbers per listing in place of six is a chart, a smaller table or fewer inputs for a model, and the share left out says what that summary costs before anything is built on it. When. Whenever a table is cut to its first few components: to choose how many to keep from what they leave out, and to see which columns and which listings the summary serves worst.

pca2 = PCA(n_components=2).fit(Z)
kept = pca2.explained_variance_ratio_.sum()
print(round(kept, 4))                 # 0.7586
short = pca2.transform(Z)
print(short.shape)                    # (3000, 2)
rebuilt = pca2.inverse_transform(short)
print(rebuilt.shape)                  # (3000, 6)
by_hand = short @ pca2.components_ + pca2.mean_
print(np.allclose(rebuilt, by_hand))  # True

mse = ((Z - rebuilt) ** 2).mean()
print(round(mse, 6))                 # 0.241361
print(round(1 - kept, 6))            # 0.241361
print(round(vals[:4].sum() / 6, 6))  # 0.241361

A listing's row minus its rebuild is its residual off the plane of the two components, Lesson 2.3.4's residual, made of its scores on the four components left out. The rotation kept every length, so the squared residual is those four squared scores added, and over the listings each squared score averages to its eigenvalue: the mean over all 18,000 entries is the four smallest eigenvalues over 6, the share left out. The match is exact because the scaler gave every column a mean square of 1, so the six columns hold 6 between them.

per_column = 1 - ((Z - rebuilt) ** 2).mean(axis=0)
print(pd.Series(per_column, index=cols).round(3))

area_sqm       0.896
rooms          0.912
bathrooms      0.780
year_built     0.656
distance_km    0.647
price          0.661
dtype: float64

Column by column, the two components rebuild from 0.647 of distance to 0.912 of rooms, and 0.759 is the average of the six. The price's 0.661 is the share of that one column rebuilt, neither 0.759 nor how well the components would predict a price.

scaler = prep[-1]
back = scaler.inverse_transform(rebuilt[[16]])[0]
listed = listings[cols].iloc[16]
print(pd.DataFrame({"listed": listed, "rebuilt": back}).round(1))

               listed   rebuilt
area_sqm        134.0     120.2
rooms             4.0       4.4
bathrooms         1.0       2.1
year_built     2020.0    2008.4
distance_km       6.8       8.1
price        722000.0  446461.2

The scaler's inverse_transform, Lesson 3.2.1's, turns standardized rows back into listing units. Two components rebuild listing 16 at 120.2 square meters, 2.1 bathrooms and 446,461: its 722,000 and its single bathroom sit in the four components left out.

raw = prep[0].transform(listings[cols])  # filled, not scaled
raw_pca = PCA().fit(raw)
print(np.round(raw_pca.explained_variance_ratio_, 3))
print(np.round(raw_pca.components_[0], 3))
print(pd.Series(scaler.scale_, index=cols).round(2))

[1. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 1.]
area_sqm           43.18
rooms               1.71
bathrooms           0.81
year_built         22.12
distance_km         3.04
price          195267.73
dtype: float64

Unscaled, the first component is the price alone, with a share of 1.0 to three decimals. PCA centers each column and never scales one, so when one column's spread in its own units dwarfs the rest, the direction of most spread is that column: 195,267.73 for the price against 43.18 square meters of area, a gap a variance squares. It is the trap Lesson 3.8.1's k-means met, whose clusters on unscaled columns split the listings by price alone.

DON'T:

• DON'T run PCA on columns measured in different units without scaling them. On the filled listings its first component came out as the price alone, with a share of 1.0 to three decimals.
• DON'T read 0.759 as the share of every column two components keep, or as what they predict. They kept 0.912 of rooms and 0.661 of the price, and what components predict about a price is a question for held-out rows, Lesson 3.8.3's.
• DON'T read inverse_transform's output in listing units. It rebuilt standard deviations, and listing 16's 722,000 came back as 446,461 only through the scaler's inverse_transform.

Lab Exercises

Lab Exercise 1: Eigenvectors by Hand, Then Three Columns That Barely Move Together

In eigen_check.py in ~/dev/bootcamp, multiply M = np.array([[5, 2], [2, 2]]) by each of [2, 1], [1, -2], [1, 1] and [-4, -2] in one loop, printing each vector beside its product. Say in a comment which are eigenvectors of M and with which eigenvalue, and show that the two you found sit at a right angle.

Then, in three_columns.py, standardize year_built, distance_km and price with the lesson's median fill and scaler on all 3,000 listings. Build their correlation matrix, check it against pandas' corr of the filled columns and print it rounded to three decimals. Run np.linalg.eigh on it, print the eigenvalues largest first and their sum, check all three pairs and the right angles, and print the largest eigenvalue's eigenvector with the column names. Print that eigenvalue's share of the total, and compare it in a comment with the six columns' 3.338 of 6.

Lab Exercise 2: Four Columns That Move Together

In four_columns.py, standardize area_sqm, rooms, bathrooms and price the same way. Fit PCA() on them and run eigh on their correlation matrix, show that components_ matches eigh's eigenvectors once the signs are dropped, and print the first component from each. Print the variance shares, and check that the components are perpendicular and of length 1 and that every listing keeps its distance from the center.

Then keep one component. Print the share it keeps, the mean squared error of the table rebuilt from it and the share of each column it rebuilds. Score the listing at position 2 on that component by hand with a dot product and with transform, and print the listing beside its rebuild in listing units. Finally, fit one component on the filled but unscaled columns and print it and its share.

import numpy as np

M = np.array([[5, 2], [2, 2]])
for v in ([2, 1], [1, -2], [1, 1], [-4, -2]):
    v = np.array(v)
    print(v, M @ v)
first, second = np.array([2, 1]), np.array([1, -2])
print(first @ second)  # 0

# [2, 1] and [-4, -2] come back 6 times themselves, eigenvalue 6, and
# [1, -2] comes back as itself, eigenvalue 1. [1, 1] comes back as
# [7, 4], 7 times one entry and 4 times the other: no eigenvector.

[2 1] [12  6]
[ 1 -2] [ 1 -2]
[1 1] [7 4]
[-4 -2] [-24 -12]

import numpy as np
import pandas as pd
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
three = ["year_built", "distance_km", "price"]
prep3 = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
Z3 = prep3.fit_transform(listings[three])
C3 = Z3.T @ Z3 / len(Z3)
filled3 = pd.DataFrame(prep3[0].transform(listings[three]), columns=three)
print(np.allclose(C3, filled3.corr()))  # True
print(np.round(C3, 3))

[[1.    0.289 0.162]
 [0.289 1.    0.019]
 [0.162 0.019 1.   ]]

vals3, vecs3 = np.linalg.eigh(C3)
pairs = np.allclose(C3 @ vecs3, vecs3 * vals3)
at_right_angles = np.allclose(vecs3.T @ vecs3, np.eye(3))
print(np.round(vals3[::-1], 3))           # [1.339 0.984 0.677]
print(round(vals3.sum(), 6))              # 3.0
print(pairs, at_right_angles)             # True True
print(pd.Series(vecs3[:, -1], index=three).round(3))
print(round(vals3[-1] / vals3.sum(), 3))  # 0.446

# 1.339 of 3 is 0.446 of the spread, against 3.338 of 6, 0.556, for the
# six columns: three columns that barely move together leave no
# direction much longer than the others.

year_built    -0.698
distance_km   -0.615
price         -0.367
dtype: float64

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
four = ["area_sqm", "rooms", "bathrooms", "price"]
prep4 = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
Z4 = prep4.fit_transform(listings[four])
vals4, vecs4 = np.linalg.eigh(Z4.T @ Z4 / len(Z4))
pca4 = PCA().fit(Z4)
comps4 = pca4.components_
order4 = vecs4[:, ::-1].T           # largest first, as rows
print(np.allclose(comps4, order4))  # False
# True
print(np.allclose(np.abs(comps4), np.abs(order4)))
print(np.round(comps4[0], 3))       # [0.524 0.529 0.489 0.454]
print(np.round(order4[0], 3))       # [-0.524 -0.529 -0.489 -0.454]
# [0.812 0.114 0.053 0.022]
print(np.round(pca4.explained_variance_ratio_, 3))
# True
print(np.allclose(comps4 @ comps4.T, np.eye(4)))
scores4 = pca4.transform(Z4)
before = np.linalg.norm(Z4, axis=1)
after = np.linalg.norm(scores4, axis=1)
print(np.allclose(before, after))   # True

one = PCA(n_components=1).fit(Z4)
kept1 = one.explained_variance_ratio_.sum()
rebuilt4 = one.inverse_transform(one.transform(Z4))
mse4 = ((Z4 - rebuilt4) ** 2).mean()
print(round(kept1, 4))  # 0.8116
print(round(mse4, 4))   # 0.1884
per_col4 = 1 - ((Z4 - rebuilt4) ** 2).mean(axis=0)
print(pd.Series(per_col4, index=four).round(3))

area_sqm     0.892
rooms        0.909
bathrooms    0.778
price        0.668
dtype: float64

by_hand = Z4[2] @ one.components_[0]
print(round(by_hand, 3))                    # 2.235
print(np.round(one.transform(Z4[[2]]), 3))  # [[2.235]]
back = prep4[-1].inverse_transform(rebuilt4[[2]])[0]
listed = listings[four].iloc[2]
print(pd.DataFrame({"listed": listed, "rebuilt": back}).round(1))
raw4 = prep4[0].transform(listings[four])
raw_one = PCA(n_components=1).fit(raw4)
print(np.round(raw_one.components_[0], 3))  # [0. 0. 0. 1.]
# [1.]
print(np.round(raw_one.explained_variance_ratio_, 3))

             listed   rebuilt
area_sqm      143.6     149.2
rooms           5.0       5.6
bathrooms       3.0       2.6
price      582000.0  578660.6

Theory Assessment

Q1: Four Products of One Matrix

For B = np.array([[4, 1], [1, 4]]), say without running anything what B @ v returns for v equal to [1, 1], [1, -1], [1, 0] and [-3, -3], which of the four are eigenvectors of B and with which eigenvalues, and how the product for [1, 0] follows from two of the others. Then explain why eigh may return either [0.707, -0.707] or [-0.707, 0.707] for one eigenvalue, and why B's two eigenvectors must sit at a right angle.

Hint: Write [1, 0] as half of [1, 1] plus half of [1, -1]. For the sign, count the vectors of length 1 on one line through the center.

Q2: What C Holds and What eigh Returns

Explain why C's diagonal is all 1s, why C equals its transpose, and why listings[cols].corr() gives area and rooms 0.924 where C gives 0.900. Then say what vals[0], vecs[:, 0], vecs[:, -1] and vecs[-1] each hold, why the six eigenvalues add up to 6, and why eigh returned 1 and 3 for lopsided without a word.

Hint: Every entry of C is two columns dotted together and divided by n. For eigh, ask which way it sorts, which axis holds one eigenvector, and which entries of a matrix it reads.

Q3: Six Components, Then Two

PCA().fit(Z) returned the first component with every weight positive where eigh's eigenvector had every weight negative, its explained_variance_ came out 1.00033344 times each eigenvalue, and two components rebuilt the table with a mean squared error of 0.241361. Explain each. Then explain what transform computes for one listing, why every listing keeps its distance from the center, and why two components keep 0.661 of the price when they keep 0.759 of the table.

Hint: One fact about eigenvectors, one divisor and one right angle cover the first three. For the price, write a standardized column as the scores times the components' weights on it, and ask what each component adds to that column's variance of 1.

Q4: The Same Table, Standardized by pandas

A teammate standardizes the six filled columns with pandas instead of the scaler, then repeats the lesson. Read the script without running it.

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer

listings = pd.read_csv("listings.csv")
cols = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km", "price"]
imputer = SimpleImputer(strategy="median").set_output(transform="pandas")
filled = imputer.fit_transform(listings[cols])
Zp = ((filled - filled.mean()) / filled.std()).to_numpy()
Cp = Zp.T @ Zp / len(Zp)
vals_p, vecs_p = np.linalg.eigh(Cp)
print(np.allclose(Cp, filled.corr()))
print(round(vals_p.sum(), 3), round(vals_p[0], 3))
pca_p = PCA().fit(Zp)
print(round(pca_p.explained_variance_.sum(), 3))
print(np.round(pca_p.explained_variance_ratio_, 3))
pair_p = PCA(n_components=2).fit(Zp)
rebuilt_p = pair_p.inverse_transform(pair_p.transform(Zp))
left_out = 1 - pair_p.explained_variance_ratio_.sum()
print(round(((Zp - rebuilt_p) ** 2).mean(), 4), round(left_out, 4))

The teammate claims five things: Cp is the correlation matrix, its eigenvalues add up to 6 as the lesson's did, vals_p[0] belongs to the main direction because it prints first, the shares differ from the lesson's because the standardizing differed, and the rebuild's error equals the share left out. Say what each print shows and which claims hold.

Hint: pandas' std divides by n minus 1, so work out each standardized column's mean square under n. Then follow that one factor through every number, asking which divide by n, which by n minus 1, and which are a ratio where it cancels.

import numpy as np

B = np.array([[4, 1], [1, 4]])
for v in ([1, 1], [1, -1], [1, 0], [-3, -3]):
    print(B @ np.array(v))
print(np.linalg.eigh(B))

[5 5]
[ 3 -3]
[4 1]
[-15 -15]
EighResult(eigenvalues=array([3., 5.]), eigenvectors=array([[-0.70710678,  0.70710678],
       [ 0.70710678,  0.70710678]]))

w = pca.components_[:2]  # the two kept rows
by_weights = (vals[::-1][:2, np.newaxis] * w ** 2).sum(axis=0)
print(pd.Series(by_weights, index=cols).round(3))
print(np.round(pca.components_[:, 5], 3))

area_sqm       0.896
rooms          0.912
bathrooms      0.780
year_built     0.656
distance_km    0.647
price          0.661
dtype: float64
[ 0.436 -0.145  0.445 -0.621 -0.447  0.068]

False
5.998 0.086
6.0
[0.556 0.202 0.134 0.061 0.032 0.014]
0.2413 0.2414

print(np.round(np.diag(Cp)[:2], 5))  # [0.99967 0.99967]
print(round(2999 / 3000, 5))         # 0.99967

Assignment

Assignment: The Listings Four Components Cannot Rebuild

Before a pricing model sees listings.csv, the data team wants a check that finds the listings whose columns break the pattern the rest share, so a person looks at them first. Build it in rebuild_check.py on the lesson's six columns, fill and scaler, from a PCA that keeps the fewest components whose shares reach 0.9. PCA(n_components=0.9) does that: a float between 0 and 1 in place of a count keeps components until their shares reach it, n_components_ reports how many the fit kept, and n_components keeps the 0.9 you gave.

Give every listing a rebuild error, the mean of its six squared misses, and show that the errors average to the share left out. Report the five worst with their listed area, rooms, bathrooms and price, their error, the column with the largest squared miss, and passes_z, True when Lesson 2.5.4's |z| > 3 check flags none of their six standardized values. Print the areas the components rebuild for those five in square meters, and the mean error of the 125 listings with no area beside the rest's.

Then test what the report suggests. Refill each missing area with the median area of the listings with the same number of rooms, Lesson 2.6.2's groupby and transform, on all 3,000 listings like every step in this lesson, and run the check again. Print the three refilled areas, their errors before and after, the 125's mean error, and the new five worst. End with a comment on what the first five had in common, what a column-by-column check missed, and what the refill showed.

Hints:

• np.argsort returns positions in error, so read the same positions from the table with .iloc, Lesson 2.4.2's distinction, and index Z and the rebuild with them directly.
• Write the check as a function of the table, so the refilled table goes through exactly the same fill, scaler and PCA(n_components=0.9), and return everything you print from.

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
cols = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km", "price"]


def rebuild_check(table):
    prep = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
    Z = prep.fit_transform(table[cols])
    check = PCA(n_components=0.9).fit(Z)
    rebuilt = check.inverse_transform(check.transform(Z))
    # one error per listing, the mean of its six squared misses
    error = ((Z - rebuilt) ** 2).mean(axis=1)
    return error, check, prep, Z, rebuilt


error, check, prep, Z, rebuilt = rebuild_check(listings)
kept = check.explained_variance_ratio_.sum()
print(check.n_components, check.n_components_)
print(round(kept, 4))
print(round(error.mean(), 4), round(1 - kept, 4))

0.9 4
0.9538
0.0462 0.0462

worst = np.argsort(error)[::-1][:5]  # positions, worst first
misses = (Z[worst] - rebuilt[worst]) ** 2
shown = ["area_sqm", "rooms", "bathrooms", "price"]
report = listings[shown].iloc[worst].copy()
report["error"] = error[worst].round(3)
report["largest_miss"] = [cols[i] for i in misses.argmax(axis=1)]
report["passes_z"] = (np.abs(Z[worst]) <= 3).all(axis=1)
print(report)
back = prep[-1].inverse_transform(rebuilt[worst])
print(np.round(back[:, 0], 1))
no_area = listings["area_sqm"].isna().to_numpy()
print(no_area.sum(), round(error[no_area].mean(), 3))
print(round(error[~no_area].mean(), 3))

      area_sqm  rooms  bathrooms    price  error largest_miss  passes_z
534      117.1      4          2  1930000  1.954     area_sqm     False
2675       NaN      9          4   941000  1.206     area_sqm     False
2775       NaN      8          4   868000  1.135     area_sqm      True
1497       NaN      7          4   797000  0.836     area_sqm      True
1745      63.2      3          1  1255000  0.824     area_sqm     False
[211.5 186.8 183.1 167.6 135. ]
125 0.154
0.042

room_median = listings.groupby("rooms")["area_sqm"].transform("median")
refilled = listings.assign(area_sqm=listings["area_sqm"].fillna(room_median))
filled_rows = [2675, 2775, 1497]  # positions of the three
print(refilled["area_sqm"].iloc[filled_rows].tolist())
error2, check2, prep2, Z2, rebuilt2 = rebuild_check(refilled)
print(check2.n_components_, round(check2.explained_variance_ratio_.sum(), 4))
print(np.round(error[filled_rows], 3), np.round(error2[filled_rows], 3))
print(round(error2[no_area].mean(), 3), round(error2[~no_area].mean(), 3))
worst2 = np.argsort(error2)[::-1][:5]
print(listings[shown].iloc[worst2])
print(np.round(error2[worst2], 3))

# All five missed most in area. Three had no area, and the median fill
# wrote 89.2 square meters beside 7 to 9 rooms, where the components
# rebuild 167.6 to 186.8, and two are priced far above their size.
# 2775 and 1497 have no |z| above 3 in any column: each value is
# ordinary alone, and only the combination breaks the pattern. Refilled
# from listings with the same rooms, the three fell to 0.005, 0.036 and
# 0.088: the median fill was what the check had found.

[230.0, 211.25, 184.0]
4 0.9595
[1.206 1.135 0.836] [0.005 0.036 0.088]
0.034 0.041
      area_sqm  rooms  bathrooms    price
534      117.1      4          2  1930000
1745      63.2      3          1  1255000
1389      67.0      3          1  1233000
1276     111.4      4          2  1377000
1310     117.9      5          3   908000
[1.811 0.739 0.654 0.584 0.542]

Additional Resources

Optional, and nothing later depends on either. Worth opening for the other forms of PCA, or for the same components found from the table itself.

• Official docs - scikit-learn 1.9 User Guide, "Decomposing signals in components (matrix factorization problems)", the section "Principal component analysis (PCA)". whiten=True, which rescales each component's scores to a variance of 1, IncrementalPCA for a table too large to hold in memory, and the randomized solver that finds only the first few components.
• Wikipedia - "Principal component analysis", the section "Singular value decomposition". The components and scores read from the table's singular value decomposition with no covariance matrix formed, which the article calls the standard way to compute PCA unless only a few components are needed.

Summary

An eigenvector turned out to be the plainest thing a matrix can do to a direction: leave it on its line and scale it. A symmetric matrix splits every vector into pieces along perpendicular directions and stretches each by its own number, so a table of entries becomes a short list of directions and stretches.

PCA is that idea aimed at a table. The correlation matrix holds how every pair of columns moves together, on standardized columns its eigenvectors are the directions the listings spread along, and its eigenvalues divide the total spread among them. Turning the axes onto those directions moves no listing, only the numbers that describe it, and on these listings more than half the spread lay in one direction, where area, rooms, bathrooms and price rise as one.

Most of what the lesson ran into was convention: a sign chosen so that a library prints the same answer twice, an order sorted one way here and the other way there, a divisor of n or of n minus 1. None of them moves a direction and each moves a printed number, so every comparison was made up to sign, in a stated order, with the divisor named. The scale of the columns is no convention, because it moves the directions themselves, and the share fewer components keep belongs to the whole table, never to each column, and never to what a model could predict.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

C = Z.T @ Z / len(Z)                 # correlation, scaler's Z
vals, vecs = np.linalg.eigh(C)       # ascending, columns
np.allclose(C @ vecs, vecs * vals)   # every pair at once
pca = PCA().fit(Z)                   # rows, largest first
np.allclose(np.abs(pca.components_), np.abs(vecs[:, ::-1].T))
scores = pca.transform(Z)            # Z @ components_.T
pca2 = PCA(n_components=2).fit(Z)
rebuilt = pca2.inverse_transform(pca2.transform(Z))
((Z - rebuilt) ** 2).mean()          # the share left out
prep[-1].inverse_transform(rebuilt)  # back to listing units