Lecture 25 – Data 100, Summer 2025

Data 100, Summer 2025

Acknowledgments Page

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import yaml
from datetime import datetime
from ds100_utils import *
import plotly.express as px

# Reduce number of sigfigs shown by numpy
np.set_printoptions(precision=2, suppress=True)

# Reduce number of sigfigs shown by pandas
pd.set_option('display.float_format', lambda x: '%.2f' % x)

PCA with SVD

Looking at this rectangle data, we can see that it is rank 3, since perimeter is a linear combination of width and height.

Area is not a linear combination of width and height, but we might surmise that area does not provide a lot of additional information beyond width and height. Let's see if PCA picks up on this.

rectangle = pd.read_csv("data/rectangle_data.csv")
rectangle

	width	height	area	perimeter
0	8	6	48	28
1	2	4	8	12
2	1	3	3	8
3	9	3	27	24
4	9	8	72	34
...	...	...	...	...
95	8	5	40	26
96	8	7	56	30
97	1	4	4	10
98	1	6	6	14
99	2	6	12	16

100 rows × 4 columns

Step 1: Center the Data Matrix $X$

Keep in mind that sklearn centers data by default when fitting PCA. Here, we are doing the linear algebra by hand.

X_centered = rectangle - np.mean(rectangle, axis = 0)
X_centered.head(10)

	width	height	area	perimeter
0	2.97	1.35	24.78	8.64
1	-3.03	-0.65	-15.22	-7.36
2	-4.03	-1.65	-20.22	-11.36
3	3.97	-1.65	3.78	4.64
4	3.97	3.35	48.78	14.64
5	-2.03	-3.65	-20.22	-11.36
6	-1.03	-2.65	-15.22	-7.36
7	0.97	0.35	6.78	2.64
8	1.97	-3.65	-16.22	-3.36
9	2.97	-2.65	-7.22	0.64

In situations where the units are on different scales, it is useful to normalize (i.e., standardize) the data before performing SVD. This can be done by dividing each column by its standard deviation.

• This puts every column on a standard deviation scale. A value of 1 implies the entry is 1 standard deviation higher than its mean.

X = X_centered / np.std(X_centered, axis = 0)
X.head(10)

	width	height	area	perimeter
0	1.07	0.58	1.35	1.21
1	-1.09	-0.28	-0.83	-1.03
2	-1.45	-0.71	-1.10	-1.59
3	1.43	-0.71	0.21	0.65
4	1.43	1.45	2.65	2.05
5	-0.73	-1.58	-1.10	-1.59
6	-0.37	-1.15	-0.83	-1.03
7	0.35	0.15	0.37	0.37
8	0.71	-1.58	-0.88	-0.47
9	1.07	-1.15	-0.39	0.09

Step 2: Get the SVD of standardized $X$

The np.linalg.svd function computes the SVD of an inputted X matrix.

U, S, Vt = np.linalg.svd(X, full_matrices = False)

full_matrices = False truncates the number of columns of U to the rank of X to avoid unnecessary computation. PCA does not use more columns of U than the rank of X. This is sometimes called the "economy" SVD. The slides use the dimensions of the economy SVD. Don't worry about these details for Data 100! Just include the argument.

SVD dimensions:

print("Shape of U", U.shape)
print("Shape of S", S.shape)
print("Shape of Vt", Vt.shape)

Shape of U (100, 4)
Shape of S (4,)
Shape of Vt (4, 4)

print('First 10 rows of U. The 4 cols are the latent features but expressed as length 1 vectors.')
print(U[:10, :])
print()

print('S. The 4 singular values that "scale up" the 4 cols of U into the 4 latent features (Z).')
print(S)
print()

print('Vt. The 4 rows are the principal components. "Recipes" for combining the 4 real features into each of the 4 latent features. Rows and columns are unit vectors.')
print(Vt)
print()

First 10 rows of U. The 4 cols are the latent features but expressed as length 1 vectors.
[[-0.13  0.02  0.05  0.61]
 [ 0.1  -0.05  0.04  0.04]
 [ 0.14 -0.04  0.12 -0.06]
 [-0.05  0.15 -0.1  -0.08]
 [-0.23 -0.02  0.18 -0.11]
 [ 0.15  0.07  0.14 -0.05]
 [ 0.1   0.06  0.06 -0.04]
 [-0.04  0.01  0.   -0.07]
 [ 0.06  0.16 -0.1  -0.1 ]
 [ 0.01  0.16 -0.11 -0.07]]

S. The 4 singular values that "scale up" the 4 cols of U into the 4 latent features (Z).
[16.99 10.13  2.94  0.  ]

Vt. The 4 rows are the principal components. "Recipes" for combining the 4 real features into each of the 4 latent features. Rows and columns are unit vectors.
[[-0.44 -0.38 -0.57 -0.58]
 [ 0.65 -0.76 -0.02  0.02]
 [-0.28 -0.28  0.82 -0.41]
 [ 0.55  0.46 -0.   -0.7 ]]

$S$ is a little different in NumPy. Since the only useful values in the diagonal matrix $S$ are the singular values on the diagonal axis, only those values are returned and they are stored in an array.

If we want the diagonal elements:

# np.diag makes a diagonal matrix from the vector S
Sm = np.diag(S)
Sm

array([[16.99,  0.  ,  0.  ,  0.  ],
       [ 0.  , 10.13,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  2.94,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ]])

Computing the contribution to the total variance:

pd.DataFrame(S**2 / np.sum(S**2))

	0
0	0.72
1	0.26
2	0.02
3	0.00

Now we see that 72% and 26% of the variance is in the first two PC dimensions, respectively, which makes sense since rectangles are largely described by height and length.

• Area is not a linear combination of height and length, so its contribution is non-zero but very small.

• Perimeter is a linear combination of height and length, so its corresponding singular value is 0.

The information below is only relevant if you print out all digits with numpy. We set an option at the top of the notebook to only shown two decimal places.

Hmm, looks like are four diagonal entries are not zero. What happened?

It turns out there were some numerical rounding errors, but the last value is so small ($10^{-15}$) that it's practically $0$.

np.isclose(S[3], 0)

np.True_

S.round(5)

array([16.99, 10.13,  2.94,  0.  ])

pd.DataFrame(np.round(np.diag(S),3))

	0	1	2	3
0	16.99	0.00	0.00	0.00
1	0.00	10.13	0.00	0.00
2	0.00	0.00	2.94	0.00
3	0.00	0.00	0.00	0.00

Step 3 Computing Approximations to the Data

Let's try to approximate the data X in two dimensions.

Using $Z = X * V$

Recall that the columns of Z are the latent features.

The first column of Z is the latent feature with the largest variance, and the second column of Z is the latent feature with the second largest variance that is orthogonal to the first column.

In this example, Z has the same dimensions as the first two columns of X.

# We can construct Z using the V matrix (transpose Vt!)
# The columns of V are the PCs, so the rows of Vt are the PCs.

print('X (truncated):')
print(X.head())
print()

print('Vt:')
print(Vt)
print()

X (truncated):
   width  height  area  perimeter
0   1.07    0.58  1.35       1.21
1  -1.09   -0.28 -0.83      -1.03
2  -1.45   -0.71 -1.10      -1.59
3   1.43   -0.71  0.21       0.65
4   1.43    1.45  2.65       2.05

Vt:
[[-0.44 -0.38 -0.57 -0.58]
 [ 0.65 -0.76 -0.02  0.02]
 [-0.28 -0.28  0.82 -0.41]
 [ 0.55  0.46 -0.   -0.7 ]]

# Construct Z using only the first two PCs
Z = X.to_numpy() @ Vt.T[:,:2]
pd.DataFrame(Z).head(10)

	0	1
0	-2.17	0.25
1	1.66	-0.50
2	2.46	-0.42
3	-0.86	1.48
4	-3.88	-0.18
5	2.47	0.71
6	1.67	0.62
7	-0.64	0.11
8	1.06	1.67
9	0.13	1.58

Using $Z = U * S$

Recall that $Z = XV = (USV^T)V = US$, since V is orthonormal.

print('First two columns of U (truncated):')
print(U[:10, :2])
print()

print('S:')
np.diag(S[:2])

First two columns of U (truncated):
[[-0.13  0.02]
 [ 0.1  -0.05]
 [ 0.14 -0.04]
 [-0.05  0.15]
 [-0.23 -0.02]
 [ 0.15  0.07]
 [ 0.1   0.06]
 [-0.04  0.01]
 [ 0.06  0.16]
 [ 0.01  0.16]]

S:

array([[16.99,  0.  ],
       [ 0.  , 10.13]])

Construct Z using the first two columns of U and the first two singular values:

Z = U[:, :2] @ np.diag(S[:2])
print(Z.shape)
pd.DataFrame(Z).head(10)

(100, 2)

	0	1
0	-2.17	0.25
1	1.66	-0.50
2	2.46	-0.42
3	-0.86	1.48
4	-3.88	-0.18
5	2.47	0.71
6	1.67	0.62
7	-0.64	0.11
8	1.06	1.67
9	0.13	1.58

The columns of U are just the normalized (i.e., length 1) columns of Z:

# Normalize first column of Z using L2 norm
length_of_col_1 = np.sqrt(np.sum(Z[:, 0]**2))
normed_z = Z[:, 0] / length_of_col_1
print(normed_z[:10])

print(U[:10, 0])

[-0.13  0.1   0.14 -0.05 -0.23  0.15  0.1  -0.04  0.06  0.01]
[-0.13  0.1   0.14 -0.05 -0.23  0.15  0.1  -0.04  0.06  0.01]

This implies that the singular values are just the length of the column vectors of Z:

# length of first column of Z (L2 norm)
print(np.sqrt(np.sum(Z[:, 0]**2)))

# Identical to code above
print(np.linalg.norm(Z[:, 0]))

# Print first singular value
print(S[0])

16.989286229589222
16.989286229589222
16.98928622958922

We get the same results if we fit PCA with scikit-learn:

from sklearn.decomposition import PCA

# This code computes first two columns of Z (i.e., the first two latent features)
# And, yes, this whole lecture can be summarized by these two lines of code! 

# Initialize a PCA model object with 2 components
pca = PCA(2)

# Fit the PCA model to the data
pd.DataFrame(pca.fit_transform(X)).head(10)

	0	1
0	2.17	-0.25
1	-1.66	0.50
2	-2.46	0.42
3	0.86	-1.48
4	3.88	0.18
5	-2.47	-0.71
6	-1.67	-0.62
7	0.64	-0.11
8	-1.06	-1.67
9	-0.13	-1.58

The Z we computed is identical to the one from sklearn:

pd.DataFrame(Z).head(10)

	0	1
0	-2.17	0.25
1	1.66	-0.50
2	2.46	-0.42
3	-0.86	1.48
4	-3.88	-0.18
5	2.47	0.71
6	1.67	0.62
7	-0.64	0.11
8	1.06	1.67
9	0.13	1.58

Notice that the covariance matrix of Z is diagonalized, since the latent features are uncorrelated, unlike the original features.

• In other words, the off-diagonal elements are 0 since the covariance between features is 0.

• The diagonal elements are the variance of each latent feature

print('Covariance matrix of Z is diagonalized, since latent features are uncorrelated:')
print(pd.DataFrame(np.cov(Z.T)))
print()

print('Covariance matrix of X is NOT diagonalized, since original features are correlated:')
print(pd.DataFrame(np.cov(X.T)))

Covariance matrix of Z is diagonalized, since latent features are uncorrelated:
     0    1
0 2.92 0.00
1 0.00 1.04

Covariance matrix of X is NOT diagonalized, since original features are correlated:
      0     1    2    3
0  1.01 -0.03 0.69 0.77
1 -0.03  1.01 0.62 0.63
2  0.69  0.62 1.01 0.94
3  0.77  0.63 0.94 1.01

Lower Rank Approximation of X

Let's now try to recover X from our approximation.

In other words, we do the reverse transformation: Transform our latent features back to the original scale of the data, and then see how close we are to the original data.

rectangle.head()

	width	height	area	perimeter
0	8	6	48	28
1	2	4	8	12
2	1	3	3	8
3	9	3	27	24
4	9	8	72	34

# Use two principal components
k = 2

U, S, Vt = np.linalg.svd(X, full_matrices = False)

## Construct the latent features
Z = U[:,:k] @ np.diag(S[:k])

## Approximate the original rectangle using the latent features Z and the principle components.
# Remember that X = USVt = ZVt. If we only use the first two columns of U, 
# first two singular values, and first two principal components, this
# equation becomes an approximation of the original data X. 
rectangle_hat = pd.DataFrame(Z @ Vt[:k, :], columns = rectangle.columns)

## Scale and shift the factors back to the original coordinate system.
# Recall that we standardized the original data by subtracting the mean
# and dividing by the SD. We do this in reverse to get back to the natural scale.
rectangle_hat = rectangle_hat * np.std(rectangle, axis = 0) + np.mean(rectangle, axis = 0)

print("Shape of approximated data:", rectangle_hat.shape)
rectangle_hat.head(10)

Shape of approximated data: (100, 4)

	width	height	area	perimeter
0	8.11	6.09	45.86	28.41
1	2.10	4.09	6.01	12.38
2	1.29	3.24	-2.47	9.05
3	8.76	2.80	31.54	23.13
4	9.41	8.34	64.12	35.51
5	3.32	1.26	-3.08	9.17
6	4.14	2.11	5.40	12.50
7	6.01	5.01	29.86	22.03
8	6.77	0.81	11.31	15.17
9	7.73	1.78	21.10	19.02

## Plot the data
fig = px.scatter_3d(rectangle, x="width", y="height", z="area",
                    width=800, height=600)
fig.add_scatter3d(x=rectangle_hat["width"], 
                  y=rectangle_hat["height"], 
                  z=rectangle_hat["area"], 
                  mode="markers", name = "approximation")

fig.update_layout(scene=dict(
  xaxis=dict(title=dict(font=dict(size=22))),
  yaxis=dict(title=dict(font=dict(size=22))),
  zaxis=dict(title=dict(font=dict(size=22)))
))