# Principal Component Analysis and 3D Ellipsoids
Related notes:
- [[MATH310S26-Day1-Notes]]
- [[MATH310S26-Day2-Notes]]
- [[MATH310S26-Day3-WorkdayMaterials (Two-by-two matrices, normal equations, and ordinary least squares)]]
- [[MATH310S26-Day4-Notes]]
- [[MATH310S26-Day5-Notes]]
## Day 3 recap (see [[MATH310S26-Day3-WorkdayMaterials (Two-by-two matrices, normal equations, and ordinary least squares)]])
- Used [normal equations](https://en.wikipedia.org/wiki/Normal_equation) to find best-fit line for 2D scatter data via [linear regression](https://en.wikipedia.org/wiki/Linear_regression) (see [[MATH310S26-Day3-WorkdayMaterials (Two-by-two matrices, normal equations, and ordinary least squares)]])
- Worked with 7,542 points in 3D space from S26_Jan_16_IW_Data.csv
- Projected 3D data onto 2D (columns 1 and 2) for linear regression
## Day 4 recap (see [[MATH310S26-Day4-Notes]])
- Reviewed least squares and [normal equations](https://en.wikipedia.org/wiki/Normal_equation) for $y=\beta_0+\beta_1x$ via orthogonal projection.
- Extended feature sets (polynomials in $x$) while remaining linear in parameters; cautioned about [overfitting](https://en.wikipedia.org/wiki/Overfitting) and [Runge's phenomenon](https://en.wikipedia.org/wiki/Runge%27s_phenomenon).
- Defined $R^2 = 1 - \frac{RSS}{TSS}$ and interpreted $R^2=1$ (perfect fit) and $R^2=0$ (no better than predicting $\bar y$).
- Introduced mean centering: $\tilde x = x-\bar x$, $\tilde y = y-\bar y$.
- Showed effects on normal equations: $\tilde X^T\tilde X$ becomes diagonal, intercept vanishes ($\beta_0=0$), slope and $R^2$ unchanged.
- Previewed eigenanalysis/PCA for “flat like a pancake” 3D data.
## Day 5 recap (see [[MATH310S26-Day5-Notes]])
- Revisited mean centering with center-of-mass intuition; best-fit line shifts but slope is unchanged.
- Built the 2D [covariance matrix](https://en.wikipedia.org/wiki/Covariance_matrix) from dot products (variances on diagonal, covariances off-diagonal).
- Noted properties: symmetric, real/nonnegative eigenvalues, orthogonal eigenvectors.
- Linked eigenanalysis to data geometry: eigenvectors give principal directions; axis lengths scale with $\sqrt{\lambda}$.
- Toy dataset PCA: $\lambda_1\approx7.8956$, $\lambda_2\approx0.0211$ → variance explained ≈99.73% vs 0.27%.
- 3D dataset PCA: PC1 ≈58%, PC2 ≈41%, PC3 <1% → strong justification for 3D→2D [dimensionality reduction](https://en.wikipedia.org/wiki/Dimensionality_reduction).
- Ellipsoid visualization and σ-level coverage; symmetry-breaking demo (circle→ellipse) clarifies principal directions.
- Launched “Mystery Data” weekend challenge using real-world measurements.
## Day 6 goals
- Find the directions of maximum variance in 3D data ([principal components](https://en.wikipedia.org/wiki/Principal_component_analysis))
- Quantify how much [variance](https://en.wikipedia.org/wiki/Variance) each direction explains
- Understand [ellipsoid](https://en.wikipedia.org/wiki/Ellipsoid) coverage in 3D space
- See why [dimensionality reduction](https://en.wikipedia.org/wiki/Dimensionality_reduction) from 3D to 2D is justified
# Finding principal components through eigendecomposition
[Principal Component Analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis) finds orthogonal directions that maximize variance. The first PC captures the most variance, the second PC captures the most remaining variance (orthogonal to the first), and so on. This is obtained via [eigendecomposition](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix) of the covariance matrix.
For our 3D data:
1. Center the data by subtracting the mean
2. Compute the [covariance matrix](https://en.wikipedia.org/wiki/Covariance_matrix)
3. Find [eigenvalues and eigenvectors](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) of the covariance matrix
4. Eigenvalues = variance along each principal component
5. Eigenvectors = directions of principal components
## Python (NumPy)
Step 1. Load the 3D data
```python
import numpy as np
data = np.loadtxt("S26_Jan_16_IW_Data.csv", delimiter=",")
print(data.shape) # expect (7543, 3)
```
Purpose: Read our 3D dataset with no headers.
Step 2. Mean-center the data
```python
mean = np.mean(data, axis=0)
X_centered = data - mean
print(f"Mean: [{mean[0]:.1f}, {mean[1]:.1f}, {mean[2]:.1f}]")
```
Purpose: PCA requires centered data. The mean should be near [69.5, 72.0, 100.0].
Step 3. Compute the covariance matrix
```python
n = data.shape[0]
cov = (X_centered.T @ X_centered) / (n - 1)
print(f"Covariance matrix shape: {cov.shape}") # expect (3, 3)
```
Purpose: The 3×3 covariance matrix captures relationships between dimensions.
Step 4. Find eigenvalues and eigenvectors
```python
eigenvalues, eigenvectors = np.linalg.eigh(cov)
# Sort by eigenvalue (descending)
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print(f"Eigenvalues: {eigenvalues.round(1)}")
```
Purpose: Eigenvalues tell us variance along each PC. Expect approximately [101, 72, 2].
Step 5. Calculate variance explained
```python
var_explained = eigenvalues / eigenvalues.sum() * 100
cumulative = np.cumsum(var_explained)
for i in range(3):
print(f"PC{i+1}: {var_explained[i]:.1f}% (total: {cumulative[i]:.1f}%)")
```
Purpose: See that PC1 and PC2 capture ~99% of variance, justifying 3D→2D reduction. See [explained variation](https://en.wikipedia.org/wiki/Explained_variation).
Step 6. Simple scatter plot of centered data
```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_centered[:, 0], X_centered[:, 1], X_centered[:, 2],
s=1, alpha=0.3, color='blue')
ax.set_xlabel('X1-mean'); ax.set_ylabel('X2-mean'); ax.set_zlabel('X3-mean')
ax.set_title('Centered Data Cloud')
plt.show()
```
Purpose: Visualize the centered cloud of points before adding structure.
Step 7. Add principal axes (eigenvector lines)
```python
# Replot the scatter
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_centered[:, 0], X_centered[:, 1], X_centered[:, 2],
s=1, alpha=0.2, color='blue')
# Add eigenvector lines scaled by ±sqrt(eigenvalue)
colors = ['red', 'green', 'orange']
for i in range(3):
v = eigenvectors[:, i] * np.sqrt(eigenvalues[i])
ax.plot([-v[0], v[0]], [-v[1], v[1]], [-v[2], v[2]],
color=colors[i], linewidth=3, label=f'PC{i+1}')
ax.legend()
ax.set_xlabel('X1-mean'); ax.set_ylabel('X2-mean'); ax.set_zlabel('X3-mean')
plt.show()
```
Purpose: Visualize principal directions scaled by their [standard deviations](https://en.wikipedia.org/wiki/Standard_deviation).
Step 8. Add 2σ ellipsoid (transparent)
```python
# Create ellipsoid surface using parametric equations
u = np.linspace(0, 2 * np.pi, 30)
v = np.linspace(0, np.pi, 20)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Scale by 2σ (k=2) along principal axes
k = 2 # For 2-sigma ellipsoid
ellipsoid = np.zeros((3, x.shape[0], x.shape[1]))
for i in range(x.shape[0]):
for j in range(x.shape[1]):
point = np.array([x[i,j], y[i,j], z[i,j]]) * k
# Transform by eigenvalues and eigenvectors
scaled = eigenvectors @ (point * np.sqrt(eigenvalues))
ellipsoid[:, i, j] = scaled
# Plot scatter, axes, and ellipsoid
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_centered[:, 0], X_centered[:, 1], X_centered[:, 2],
s=1, alpha=0.2, color='blue')
# Principal axes
for i in range(3):
v = eigenvectors[:, i] * np.sqrt(eigenvalues[i]) * k
ax.plot([-v[0], v[0]], [-v[1], v[1]], [-v[2], v[2]],
color=colors[i], linewidth=2)
# Ellipsoid surface
ax.plot_surface(ellipsoid[0], ellipsoid[1], ellipsoid[2],
alpha=0.2, color='yellow')
ax.set_xlabel('X1-mean'); ax.set_ylabel('X2-mean'); ax.set_zlabel('X3-mean')
ax.set_title('2σ Ellipsoid')
plt.show()
```
Purpose: Visualize the 2σ boundary that should contain ~74% of points under a [multivariate normal distribution](https://en.wikipedia.org/wiki/Multivariate_normal_distribution).
Step 9. Count points in 3D ellipsoids
```python
# Mahalanobis distance for ellipsoid membership
cov_inv = np.linalg.inv(cov)
distances_sq = np.sum((X_centered @ cov_inv) * X_centered, axis=1)
for k in [1, 2, 3]:
inside = np.sum(distances_sq <= k**2)
percent = inside / n * 100
print(f"{k}σ ellipsoid contains {percent:.1f}% of points")
```
Purpose: Verify coverage using the [Mahalanobis distance](https://en.wikipedia.org/wiki/Mahalanobis_distance): 1σ ~21%, 2σ ~74%, 3σ ~97% in 3D.
## R
Step 1. Load the 3D data
```r
data <- read.csv("S26_Jan_16_IW_Data.csv", header = FALSE)
names(data) <- c("X1", "X2", "X3")
dim(data) # expect 7543 3
```
Purpose: Load the headerless CSV and name the columns.
Step 2. Mean-center the data
```r
mean_vec <- colMeans(data)
X_centered <- scale(data, center = TRUE, scale = FALSE)
round(mean_vec, 1) # expect near [69.5, 72.0, 100.0]
```
Purpose: Center each column by subtracting its mean.
Step 3. Compute the covariance matrix
```r
cov_mat <- cov(X_centered)
dim(cov_mat) # expect 3 3
```
Purpose: R’s cov() uses the unbiased estimator (n-1 divisor).
Step 4. Find eigenvalues and eigenvectors
```r
eigen_decomp <- eigen(cov_mat)
eigenvalues <- eigen_decomp$values # already sorted descending
eigenvectors <- eigen_decomp$vectors
round(eigenvalues, 1) # expect approximately [101, 72, 2]
```
Purpose: R’s eigen() returns values sorted in descending order.
Step 5. Calculate variance explained
```r
var_explained <- eigenvalues / sum(eigenvalues) * 100
cumulative <- cumsum(var_explained)
for(i in 1:3) {
cat(sprintf("PC%d: %.1f%% (total: %.1f%%)\n",
i, var_explained[i], cumulative[i]))
}
```
Purpose: Confirm that first two PCs explain ~99% of variance.
Step 6. Simple 3D scatter plot
```r
library(scatterplot3d)
s3d <- scatterplot3d(X_centered[,1], X_centered[,2], X_centered[,3],
pch = 16, cex.symbols = 0.1, color = "blue",
main = "Centered Data Cloud",
xlab = "X1-mean", ylab = "X2-mean", zlab = "X3-mean",
angle = 40)
```
Purpose: Visualize the centered cloud of points before adding structure.
Step 7. Add principal axes (eigenvector lines)
```r
# Create new plot with data
s3d <- scatterplot3d(X_centered[,1], X_centered[,2], X_centered[,3],
pch = 16, cex.symbols = 0.05, color = rgb(0,0,1,0.2),
xlab = "X1-mean", ylab = "X2-mean", zlab = "X3-mean",
angle = 40)
# Add eigenvector lines scaled by ±sqrt(eigenvalue)
colors_axes <- c("red", "green", "orange")
for(i in 1:3) {
v <- eigenvectors[, i] * sqrt(eigenvalues[i])
# Convert 3D to 2D coordinates for plotting
from <- s3d$xyz.convert(-v[1], -v[2], -v[3])
to <- s3d$xyz.convert(v[1], v[2], v[3])
segments(from$x, from$y, to$x, to$y,
col = colors_axes[i], lwd = 3)
}
legend("topright", legend = paste0("PC", 1:3),
col = colors_axes, lwd = 3)
```
Purpose: Visualize principal directions scaled by their standard deviations.
Step 8. Add 2σ ellipsoid (using ellipse3d)
```r
library(rgl) # For 3D ellipsoid
# Create 2σ ellipsoid
k <- 2
ellipsoid_cov <- cov_mat * k^2
# For static plot approximation, show ellipse projections
s3d <- scatterplot3d(X_centered[,1], X_centered[,2], X_centered[,3],
pch = 16, cex.symbols = 0.05, color = rgb(0,0,1,0.2),
xlab = "X1-mean", ylab = "X2-mean", zlab = "X3-mean",
angle = 40, main = "2σ Ellipsoid")
# Add principal axes at 2σ scale
for(i in 1:3) {
v <- eigenvectors[, i] * sqrt(eigenvalues[i]) * k
from <- s3d$xyz.convert(-v[1], -v[2], -v[3])
to <- s3d$xyz.convert(v[1], v[2], v[3])
segments(from$x, from$y, to$x, to$y,
col = colors_axes[i], lwd = 2)
}
# Note: For true 3D ellipsoid, use rgl::plot3d() and rgl::ellipse3d()
```
Purpose: Visualize the 2σ boundary that should contain ~74% of points.
Step 9. Count points in 3D ellipsoids
```r
# Mahalanobis distances
distances_sq <- mahalanobis(data, mean_vec, cov_mat)
for(k in 1:3) {
inside <- sum(distances_sq <= k^2)
percent <- inside / nrow(data) * 100
cat(sprintf("%dσ ellipsoid contains %.1f%% of points\n", k, percent))
}
```
Purpose: Verify coverage: 1σ ~21%, 2σ ~74%, 3σ ~97% in 3D.
## MATLAB/Octave
Step 1. Load the 3D data
```matlab
data = csvread('S26_Jan_16_IW_Data.csv');
size(data) % expect 7543 3
```
Purpose: Load numerical data from headerless CSV.
Step 2. Mean-center the data
```matlab
mean_vec = mean(data);
X_centered = data - mean_vec;
round(mean_vec) % expect near [70, 72, 100]
```
Purpose: Center by subtracting mean from each row.
Step 3. Compute the covariance matrix
```matlab
cov_mat = cov(data); # MATLAB uses n-1 divisor
size(cov_mat) % expect 3 3
```
Purpose: Get 3×3 covariance matrix.
Step 4. Find eigenvalues and eigenvectors
```matlab
[V, D] = eig(cov_mat);
eigenvalues = diag(D);
[eigenvalues, idx] = sort(eigenvalues, 'descend');
eigenvectors = V(:, idx);
round(eigenvalues) # expect approximately [101; 72; 2]
```
Purpose: Extract and sort eigenvalues/vectors.
Step 5. Calculate variance explained
```matlab
var_explained = eigenvalues / sum(eigenvalues) * 100;
cumulative = cumsum(var_explained);
for i = 1:3
fprintf('PC%d: %.1f%% (total: %.1f%%)\n', ...
i, var_explained(i), cumulative(i));
end
```
Purpose: Quantify variance captured by each PC.
Step 6. Simple 3D scatter plot
```matlab
figure;
scatter3(X_centered(:,1), X_centered(:,2), X_centered(:,3), 1, 'b', 'filled');
xlabel('X1-mean'); ylabel('X2-mean'); zlabel('X3-mean');
title('Centered Data Cloud');
grid on;
view(40, 30);
```
Purpose: Visualize the centered cloud of points before adding structure.
Step 7. Add principal axes (eigenvector lines)
```matlab
figure;
scatter3(X_centered(:,1), X_centered(:,2), X_centered(:,3), ...
1, 'b', 'filled', 'MarkerFaceAlpha', 0.2);
hold on;
% Add eigenvector lines scaled by ±sqrt(eigenvalue)
colors = {'r', 'g', [1 0.5 0]}; % red, green, orange
for i = 1:3
v = eigenvectors(:, i) * sqrt(eigenvalues(i));
plot3([-v(1) v(1)], [-v(2) v(2)], [-v(3) v(3)], ...
'Color', colors{i}, 'LineWidth', 3, ...
'DisplayName', sprintf('PC%d', i));
end
legend('Location', 'best');
xlabel('X1-mean'); ylabel('X2-mean'); zlabel('X3-mean');
grid on; hold off;
```
Purpose: Visualize principal directions scaled by their standard deviations.
Step 8. Add 2σ ellipsoid (transparent)
```matlab
figure;
scatter3(X_centered(:,1), X_centered(:,2), X_centered(:,3), ...
1, 'b', 'filled', 'MarkerFaceAlpha', 0.2);
hold on;
% Create ellipsoid using parametric form
[u,v] = meshgrid(linspace(0,2*pi,30), linspace(0,pi,20));
x = cos(u).*sin(v);
y = sin(u).*sin(v);
z = cos(v);
% Scale by 2σ along principal axes
k = 2; % 2-sigma
ellipsoid_pts = zeros(3, numel(x));
for idx = 1:numel(x)
pt = [x(idx); y(idx); z(idx)] * k;
ellipsoid_pts(:,idx) = eigenvectors * (pt .* sqrt(eigenvalues));
end
% Reshape for surface plot
X = reshape(ellipsoid_pts(1,:), size(x));
Y = reshape(ellipsoid_pts(2,:), size(x));
Z = reshape(ellipsoid_pts(3,:), size(x));
% Plot ellipsoid surface
surf(X, Y, Z, 'FaceAlpha', 0.2, 'FaceColor', 'y', 'EdgeColor', 'none');
% Add principal axes at 2σ scale
for i = 1:3
v = eigenvectors(:, i) * sqrt(eigenvalues(i)) * k;
plot3([-v(1) v(1)], [-v(2) v(2)], [-v(3) v(3)], ...
'Color', colors{i}, 'LineWidth', 2);
end
xlabel('X1-mean'); ylabel('X2-mean'); zlabel('X3-mean');
title('2σ Ellipsoid');
grid on; hold off;
```
Purpose: Visualize the 2σ boundary that should contain ~74% of points.*
**Step 9. Count points in 3D ellipsoids**
```matlab
% Mahalanobis distances
cov_inv = inv(cov_mat);
distances_sq = sum((X_centered * cov_inv) .* X_centered, 2);
for k = 1:3
inside = sum(distances_sq <= k^2);
percent = inside / size(data, 1) * 100;
fprintf('%dσ ellipsoid contains %.1f%% of points\n', k, percent);
end
```
*Purpose: Verify coverage: 1σ ~21%, 2σ ~74%, 3σ ~97% in 3D.*
# Key insights
## Variance and dimensional reduction
- [PCA](https://en.wikipedia.org/wiki/Principal_component_analysis) PC1 captures ~58% of [variance](https://en.wikipedia.org/wiki/Variance)
- PC2 captures ~41% of the variance
- PC3 captures <1% of the variance
- First two PCs capture >99% of total variance
- This mathematically justifies reducing from 3D to 2D (see [dimensionality reduction](https://en.wikipedia.org/wiki/Dimensionality_reduction))
## Ellipsoid coverage in 3D
Unlike 1D where 1σ contains 68% of data, in 3D (under a [multivariate normal distribution](https://en.wikipedia.org/wiki/Multivariate_normal_distribution) and measured via [Mahalanobis distance](https://en.wikipedia.org/wiki/Mahalanobis_distance)):
- 1σ ellipsoid: ~21% of points
- 2σ ellipsoid: ~74% of points
- 3σ ellipsoid: ~97% of points
The coverage depends on dimension! Higher dimensions need larger multipliers (in units of [standard deviation](https://en.wikipedia.org/wiki/Standard_deviation)) to capture the same percentage of data.
## Connection to Day 3
- Day 3: We projected 3D→2D by simply using columns 1 and 2 (see [[MATH310S26-Day3-WorkdayMaterials (Two-by-two matrices, normal equations, and ordinary least squares)]])
- Day 6: PCA finds the optimal 2D projection that preserves maximum [variance](https://en.wikipedia.org/wiki/Variance)
- The arbitrary projection (Day 3) kept coordinate meanings
- The PCA projection (Day 6) creates new coordinates that maximize information (see eigen/covariance discussion in [[MATH310S26-Day5-Notes]])