What PCA Keeps

The reading shows a three-dimensional cloud compressed into two coordinates. We’ll use that same cloud, fit PCA ourselves, and check what changes when we retain fewer directions.

First, center the columns. The directions must be learned from the centered data, and new observations must use those same means and directions.

from pathlib import Path
import json
import numpy as np
import matplotlib.pyplot as plt

plt.style.use(str(Path('../figure_style.mplstyle')))
data = json.loads(Path('pca-cloud.json').read_text())
observations = np.asarray(data['points'])
mean = observations.mean(axis=0)
X = observations - mean
U, singular_values, Vt = np.linalg.svd(X, full_matrices=False)
print(f'{X.shape[0]} observations, {X.shape[1]} features')
print('Singular values:', np.round(singular_values, 3))

144 observations, 3 features
Singular values: [22.294 11.032  1.827]

Encode and reconstruct

Change k to 1, 2, or 3 and rerun. The rows of Z contain the retained coordinates; X_hat puts those coordinates back into the original feature space. Three coordinates allow exact reconstruction here. With fewer, the residual should be perpendicular to every retained direction.

k = 2
V_k = Vt[:k].T
Z = X @ V_k
X_hat = Z @ V_k.T
residual = X - X_hat

assert np.allclose(Z, U[:, :k] * singular_values[:k])
assert np.allclose(residual @ V_k, 0)
error = np.linalg.norm(residual, 'fro') ** 2
tail_energy = np.sum(singular_values[k:] ** 2)
assert np.isclose(error, tail_energy)
print(f'Embedding shape: {Z.shape}; reconstruction shape: {X_hat.shape}')
print(f'Reconstruction error: {error:.4f}; discarded energy: {tail_energy:.4f}')

Embedding shape: (144, 2); reconstruction shape: (144, 3)
Reconstruction error: 3.3381; discarded energy: 3.3381

selected = data['selected']
fig = plt.figure(figsize=(9, 3.8))
left = fig.add_subplot(121, projection='3d')
left.scatter(*X.T, color='#44AA99', s=12, alpha=.7, label='Original')
left.scatter(*X_hat.T, color='#332288', s=12, marker='x', label='Reconstructed')
left.plot(*np.vstack([X[selected], X_hat[selected]]).T, color='#CC6677', lw=2)
left.set(xlabel='Feature 1', ylabel='Feature 2', zlabel='Feature 3', title=f'Reconstruction with {k} coordinates')
left.legend(fontsize=8)
right = fig.add_subplot(122)
counts = np.arange(1, X.shape[1] + 1)
errors = [np.linalg.norm(X - X @ Vt[:j].T @ Vt[:j], 'fro')**2 for j in counts]
tails = [np.sum(singular_values[j:]**2) for j in counts]
right.plot(counts, tails, color='#332288', label='Discarded singular-value energy')
right.scatter(counts, errors, color='#44AA99', zorder=3, label='Measured error')
right.set(xlabel='Retained coordinates', ylabel='Squared reconstruction error', xticks=counts, title='EYM predicts the error')
right.legend(fontsize=8)
fig.tight_layout()
plt.show()

Compare distances

PCA never increases a distance, but it can shrink different distances by different amounts. Each dot below is one pair of observations. Points on the diagonal retain their distance; points below it lose distance.

We’ll check the full identity, not just the inequality: original squared distance equals embedding squared distance plus discarded squared distance.

original_sq = np.sum((X[:, None, :] - X[None, :, :])**2, axis=2)
embedded_sq = np.sum((Z[:, None, :] - Z[None, :, :])**2, axis=2)
discarded_sq = np.sum((residual[:, None, :] - residual[None, :, :])**2, axis=2)
assert np.allclose(original_sq, embedded_sq + discarded_sq)
assert np.all(embedded_sq <= original_sq + 1e-10)

pairs = np.triu_indices(len(X), k=1)
fig, ax = plt.subplots(figsize=(5, 4))
ax.scatter(np.sqrt(original_sq[pairs]), np.sqrt(embedded_sq[pairs]), s=9, alpha=.4, color='#44AA99')
limit = np.sqrt(original_sq.max())
ax.plot([0, limit], [0, limit], color='#332288', lw=1, label='Distance preserved')
ax.set(xlabel='Original distance', ylabel='Embedding distance', title=f'Pairwise distances with {k} coordinates', aspect='equal')
ax.legend()
fig.tight_layout()
plt.show()

def nearest_neighbors(squared_distances):
    distances = squared_distances.copy()
    np.fill_diagonal(distances, np.inf)
    return distances.argmin(axis=1)

changed = nearest_neighbors(original_sq) != nearest_neighbors(embedded_sq)
print(f'Nearest neighbor changed for {changed.sum()} of {len(X)} observations.')

Nearest neighbor changed for 30 of 144 observations.

Try keeping one coordinate, then all three. Which statements hold for every pair, and which only describe the total error? Could the reconstruction error be small while the identity of a nearest neighbor changes?

For a new observation, subtract mean, multiply by V_k, and reconstruct with V_k.T. Add mean back only when returning to the original feature units.