๐ง PCA Fundamentals
๐ฏ What is PCA?
Principal Component Analysis (PCA) is a dimensionality reduction technique that transforms data into a new coordinate system where the greatest variance lies on the first axis (first principal component), the second greatest on the second axis, and so on. It's one of the most widely used techniques in data science and machine learning.
Why PCA Matters
Dimensionality reduction: Reduce features while preserving information. Visualization: Project high-dimensional data to 2D/3D. Noise reduction: Keep signal, discard noise. Speed: Faster training for ML models. Multicollinearity: Removes correlated features. Compression: Store data more efficiently.
๐ Dimensionality Reduction๐ Variance๐ Orthogonal๐ Linear
# PCA fundamentals (conceptual)
import numpy as np
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
# Load data
iris = load_iris()
X = iris.data # (150, 4)
y = iris.target
print(f"Original shape: {X.shape}")
# Apply PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
print(f"Reduced shape: {X_pca.shape}")
print(f"Explained variance ratio: {pca.explained_variance_ratio_}")
print(f"Total variance explained: {sum(pca.explained_variance_ratio_):.4f}")
# PCA in one line:
# 1. Center the data
# 2. Compute covariance matrix
# 3. Find eigenvectors/eigenvalues
# 4. Sort by eigenvalue
# 5. Project data onto top K eigenvectors
# Manual implementation
X_centered = X - X.mean(axis=0)
cov_matrix = np.cov(X_centered.T)
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
# Sort by eigenvalue (descending)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# Project
X_manual = X_centered @ eigenvectors[:, :2]
print(f"\nManual PCA shape: {X_manual.shape}")๐ Variance & Covariance
๐ฏ The Heart of PCA
Variance measures how spread out data is along a direction. Covariance measures how two variables vary together. PCA finds directions (principal components) that maximize variance and minimize covariance.
Variance & Covariance Formulas
Variance: Var(X) = (1/n) ฮฃ (x_i - ฮผ)ยฒ. Covariance: Cov(X,Y) = (1/n) ฮฃ (x_i - ฮผ_x)(y_i - ฮผ_y). Covariance matrix: nรn symmetric matrix. Diagonal: Variances. Off-diagonal: Covariances. Goal: Diagonalize the covariance matrix.
๐ Variance๐ Covariance๐ Matrix๐ฏ Diagonal
# Variance and covariance in PCA
import numpy as np
# Sample data
np.random.seed(42)
X = np.random.randn(100, 3) @ np.array([
[1.0, 0.7, 0.3],
[0.0, 0.5, 0.2],
[0.0, 0.0, 0.1]
])
print("=== Data Statistics ===")
print(f"Shape: {X.shape}")
print(f"Mean: {X.mean(axis=0)}")
print(f"Standard deviation: {X.std(axis=0)}")
# 1. Variance of each feature
print("\n=== Variance ===")
variance = X.var(axis=0)
print(f"Variance per feature: {variance}")
print(f"Total variance: {variance.sum()}")
# 2. Covariance matrix
print("\n=== Covariance Matrix ===")
cov_matrix = np.cov(X.T)
print(cov_matrix)
# Interpretation:
# - Diagonal: variance of each feature
# - Off-diagonal: covariance between features
# - Symmetric matrix
# 3. Correlation matrix (normalized covariance)
print("\n=== Correlation Matrix ===")
corr_matrix = np.corrcoef(X.T)
print(corr_matrix)
# 4. Eigendecomposition of covariance
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# Sort descending
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print("\n=== Eigendecomposition ===")
print(f"Eigenvalues: {eigenvalues}")
print(f"Sum of eigenvalues: {eigenvalues.sum()}")
print(f"Variance explained: {eigenvalues / eigenvalues.sum()}")
print(f"\nEigenvectors (columns):")
print(eigenvectors)
# 5. Verify: eigenvectors are orthogonal
print("\n=== Orthogonality Check ===")
for i in range(len(eigenvectors)):
for j in range(i+1, len(eigenvectors)):
dot = np.dot(eigenvectors[:, i], eigenvectors[:, j])
print(f"Dot(eig{i}, eig{j}) = {dot:.6f}")
# 6. Project data
print("\n=== Projected Data ===")
X_centered = X - X.mean(axis=0)
X_pca = X_centered @ eigenvectors
print(f"Original shape: {X.shape}")
print(f"Projected shape: {X_pca.shape}")
print(f"Variance of projected: {X_pca.var(axis=0)}")
print(f"Matches eigenvalues: {np.allclose(X_pca.var(axis=0), eigenvalues)}")๐งฎ Eigenvalues & Eigenvectors
๐ฏ The Mathematical Foundation
An eigenvector of a matrix is a vector that doesn't change direction when the matrix is applied. The eigenvalue is how much it's scaled. In PCA, eigenvectors of the covariance matrix are the principal components; eigenvalues are their variances.
Eigendecomposition
Definition: A v = ฮป v. A: Covariance matrix. v: Eigenvector (principal direction). ฮป: Eigenvalue (variance). Ordering: Sort by ฮป descending. Interpretation: Largest ฮป = most important direction.
๐งฎ Eigenvalue๐งญ Eigenvector๐ Variance๐ฏ Principal
# Eigenvalues and eigenvectors in PCA
import numpy as np
import matplotlib.pyplot as plt
# Generate correlated 2D data
np.random.seed(42)
n = 500
mean = [0, 0]
cov = [[3, 2], [2, 2]]
X = np.random.multivariate_normal(mean, cov, n)
# 1. Center data
X_centered = X - X.mean(axis=0)
# 2. Compute covariance matrix
cov_matrix = np.cov(X_centered.T)
# 3. Eigen decomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print("=== Covariance Matrix ===")
print(cov_matrix)
print("\n=== Eigenvalues (variances) ===")
print(f"ฮป1 = {eigenvalues[0]:.4f} (PC1 - largest variance)")
print(f"ฮป2 = {eigenvalues[1]:.4f} (PC2 - smallest variance)")
print(f"Total variance: {eigenvalues.sum():.4f}")
print(f"PC1 explains: {eigenvalues[0]/eigenvalues.sum()*100:.2f}%")
print(f"PC2 explains: {eigenvalues[1]/eigenvalues.sum()*100:.2f}%")
print("\n=== Eigenvectors (principal directions) ===")
print(f"PC1: {eigenvectors[:, 0]}")
print(f"PC2: {eigenvectors[:, 1]}")
# 4. Verify: A @ v = ฮป @ v
for i in range(2):
v = eigenvectors[:, i]
lam = eigenvalues[i]
A_v = cov_matrix @ v
lambda_v = lam * v
print(f"\nVerify PC{i+1}: A@v = {A_v.round(4)}, ฮป*v = {lambda_v.round(4)}")
# 5. Project data
X_pca = X_centered @ eigenvectors
print("\n=== Projected Data ===")
print(f"Original variance: {X_centered.var(axis=0)}")
print(f"PC1 variance: {X_pca[:, 0].var():.4f} (should โ ฮป1)")
print(f"PC2 variance: {X_pca[:, 1].var():.4f} (should โ ฮป2)")
# 6. Dimensionality reduction
X_reduced = X_pca[:, :1] # Keep only PC1
X_reconstructed = X_reduced @ eigenvectors[:, :1].T + X.mean(axis=0)
print("\n=== Dimensionality Reduction ===")
print(f"Original shape: {X.shape}")
print(f"Reduced shape: {X_reduced.shape}")
print(f"Reconstruction error: {np.mean((X - X_reconstructed)**2):.4f}")
# How much info did we lose?
loss = eigenvalues[1] / eigenvalues.sum()
print(f"\nInformation lost: {loss*100:.2f}%")
print(f"Information kept: {(1-loss)*100:.2f}%")
# Visualization
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.5)
for i, (eigval, eigvec) in enumerate(zip(eigenvalues, eigenvectors.T)):
plt.arrow(0, 0, eigvec[0]*np.sqrt(eigval)*2, eigvec[1]*np.sqrt(eigval)*2,
color=['red', 'blue'][i], width=0.05, head_width=0.3)
plt.text(eigvec[0]*np.sqrt(eigval)*2.2, eigvec[1]*np.sqrt(eigval)*2.2,
f'PC{i+1}', color=['red', 'blue'][i], fontsize=12)
plt.title('Original Data with Principal Components')
plt.axis('equal')
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
plt.scatter(X_pca[:, 0], X_pca[:, 1], alpha=0.5, c='green')
plt.title('Data in PC Space (Decorrelated)')
plt.axis('equal')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Eigenvalue decomposition properties
print("\n=== Properties ===")
print("1. Eigenvectors are orthogonal (perpendicular)")
print("2. Eigenvalues = variance along each PC")
print("3. Sum of eigenvalues = total variance")
print("4. Covariance matrix in PC space is diagonal")
print(f"5. Diagonal covariance: {np.diag(eigenvalues)}")๐ข SVD Approach
๐ฏ Singular Value Decomposition
SVD is a more numerically stable way to compute PCA. Instead of computing the covariance matrix (which squares the condition number), SVD decomposes the data matrix directly.
SVD vs Eigendecomposition
SVD: X = U ฮฃ V^T. Eigendecomposition: Cov = V ฮ V^T. Relationship: V = principal components, ฮฃยฒ / (n-1) = eigenvalues. Advantage: More numerically stable. Speed: Faster for tall matrices. Used by: scikit-learn's PCA.
๐ข SVD๐ฏ Stableโก Fast๐ง Scikit-learn
# SVD-based PCA
import numpy as np
# Generate data
np.random.seed(42)
X = np.random.randn(100, 5) @ np.diag([5, 3, 1, 0.5, 0.1])
# Center data
X_centered = X - X.mean(axis=0)
# 1. SVD approach
print("=== SVD Approach ===")
U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)
print(f"U shape: {U.shape}")
print(f"S shape: {S.shape}")
print(f"Vt shape: {Vt.shape}")
print(f"Singular values: {S}")
# Principal components are rows of Vt
# Variances are Sยฒ / (n-1)
n = X_centered.shape[0]
explained_variance = (S ** 2) / (n - 1)
print(f"Explained variance: {explained_variance}")
print(f"Explained variance ratio: {explained_variance / explained_variance.sum()}")
# Project data
X_pca_svd = U * S # or X_centered @ Vt.T
# 2. Covariance approach (comparison)
print("\n=== Covariance Approach ===")
cov_matrix = np.cov(X_centered.T)
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print(f"Eigenvalues: {eigenvalues}")
print(f"Match SVD: {np.allclose(eigenvalues, explained_variance)}")
# 3. Why SVD is better
print("\n=== SVD vs Eigendecomposition ===")
print("SVD:")
print(" - More numerically stable")
print(" - Works directly on data matrix")
print(" - Avoids squaring condition number")
print(" - Handles rank-deficient matrices")
print(" - Faster for tall matrices (n >> p)")
print("\nEigendecomposition:")
print(" - Requires computing covariance matrix")
print(" - Squares condition number (numerical issue)")
print(" - Slower for tall matrices")
# 4. Demonstrate numerical stability
print("\n=== Numerical Stability Test ===")
# Create ill-conditioned matrix
np.random.seed(42)
X_ill = np.random.randn(50, 10)
X_ill[:, 9] = X_ill[:, 0] + 1e-10 * np.random.randn(50) # Nearly dependent
# SVD
U, S, Vt = np.linalg.svd(X_ill - X_ill.mean(axis=0), full_matrices=False)
print(f"SVD singular values: {S[:5]}")
# Covariance
cov = np.cov((X_ill - X_ill.mean(axis=0)).T)
eigs = np.linalg.eigvalsh(cov)[::-1]
print(f"Cov eigenvalues: {eigs[:5]}")
# 5. Truncated SVD (for sparse matrices)
from scipy.sparse.linalg import svds
from scipy.sparse import random as sparse_random
# Sparse matrix
X_sparse = sparse_random(1000, 100, density=0.05, format='csr', random_state=42)
# Truncated SVD
U_k, S_k, Vt_k = svds(X_sparse, k=10)
# Sort descending
idx = np.argsort(S_k)[::-1]
S_k = S_k[idx]
print(f"\nTruncated SVD:")
print(f"Original shape: {X_sparse.shape}")
print(f"U shape: {U_k.shape}")
print(f"S shape: {S_k.shape}")
print(f"Vt shape: {Vt_k.shape}")
# 6. Randomized SVD (for very large matrices)
from sklearn.utils.extmath import randomized_svd
X_large = np.random.randn(10000, 100)
U_r, S_r, Vt_r = randomized_svd(X_large, n_components=10, random_state=42)
print(f"\nRandomized SVD:")
print(f"Original: {X_large.shape}")
print(f"Components: {len(S_r)}")
print(f"Speed: Much faster than exact SVD for large matrices")
# scikit-learn PCA uses SVD by default
from sklearn.decomposition import PCA
pca = PCA(n_components=3, svd_solver='full') # 'full', 'randomized', 'arpack'
pca.fit(X_centered)
print(f"\nSklearn PCA explained variance: {pca.explained_variance_}")
print(f"Matches SVD: {np.allclose(pca.explained_variance_, explained_variance[:3])}")๐ PCA Step-by-Step
๐ฏ The Complete Algorithm
PCA is a linear transformation with clearly defined steps. Understanding each step helps you apply PCA correctly and interpret the results.
The 6 Steps of PCA
1. Standardize: Center (and scale) the data. 2. Covariance: Compute covariance matrix. 3. Eigen: Find eigenvalues and eigenvectors. 4. Sort: Order by eigenvalue descending. 5. Select: Choose top K components. 6. Transform: Project data onto components.
๐ Standardize๐ Covariance๐งฎ Eigen๐ฏ Project
# PCA from scratch โ complete implementation
import numpy as np
import matplotlib.pyplot as plt
class PCAFromScratch:
def __init__(self, n_components=None):
self.n_components = n_components
self.components_ = None
self.explained_variance_ = None
self.explained_variance_ratio_ = None
self.mean_ = None
self.scale_ = None
def fit(self, X):
"""Fit PCA to data."""
n_samples, n_features = X.shape
# Step 1: Standardize (center and optionally scale)
self.mean_ = X.mean(axis=0)
X_centered = X - self.mean_
# Optional: Scale to unit variance
# self.scale_ = X.std(axis=0)
# X_centered = X_centered / self.scale_
# Step 2: Compute covariance matrix
cov_matrix = np.cov(X_centered, rowvar=False)
# Step 3: Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
# Step 4: Sort by eigenvalue descending
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# Step 5: Select top K components
if self.n_components is None:
self.n_components = n_features
self.components_ = eigenvectors[:, :self.n_components].T
self.explained_variance_ = eigenvalues[:self.n_components]
self.explained_variance_ratio_ = (
self.explained_variance_ / eigenvalues.sum()
)
return self
def transform(self, X):
"""Project data onto principal components."""
X_centered = X - self.mean_
return X_centered @ self.components_.T
def fit_transform(self, X):
"""Fit and transform in one step."""
self.fit(X)
return self.transform(X)
def inverse_transform(self, X_transformed):
"""Reconstruct original data from PCA."""
X_reconstructed = X_transformed @ self.components_ + self.mean_
return X_reconstructed
# Example with iris data
from sklearn.datasets import load_iris
iris = load_iris()
X = iris.data
# Apply PCA
pca = PCAFromScratch(n_components=2)
X_pca = pca.fit_transform(X)
print("=== PCA Results ===")
print(f"Original shape: {X.shape}")
print(f"Reduced shape: {X_pca.shape}")
print(f"\nExplained variance: {pca.explained_variance_}")
print(f"Explained variance ratio: {pca.explained_variance_ratio_}")
print(f"Total variance explained: {pca.explained_variance_ratio_.sum():.4f}")
print(f"\nPrincipal components (top 2):")
print(pca.components_)
# Step-by-step demonstration
print("\n\n=== Step-by-Step Demonstration ===")
# Step 1: Center
X_centered = X - X.mean(axis=0)
print(f"Step 1 - Centered mean: {X_centered.mean(axis=0).round(10)}")
# Step 2: Covariance
cov_matrix = np.cov(X_centered, rowvar=False)
print(f"\nStep 2 - Covariance matrix shape: {cov_matrix.shape}")
print(f"Covariance matrix:\n{cov_matrix}")
# Step 3: Eigen
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
print(f"\nStep 3 - Eigenvalues (unsorted): {eigenvalues}")
print(f"Eigenvectors shape: {eigenvectors.shape}")
# Step 4: Sort
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print(f"\nStep 4 - Sorted eigenvalues: {eigenvalues}")
# Step 5: Select top K
k = 2
top_k_components = eigenvectors[:, :k]
print(f"\nStep 5 - Top {k} components shape: {top_k_components.shape}")
# Step 6: Project
X_pca_step = X_centered @ top_k_components
print(f"\nStep 6 - Projected data shape: {X_pca_step.shape}")
# Verify match
print(f"\nMatches our implementation: {np.allclose(X_pca, X_pca_step)}")
# Reconstruction
X_reconstructed = X_pca @ pca.components_ + pca.mean_
reconstruction_error = np.mean((X - X_reconstructed) ** 2)
print(f"\nReconstruction error: {reconstruction_error:.4f}")
# Visualize
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(X[:, 0], X[:, 1], c=iris.target, cmap='viridis', alpha=0.7)
plt.xlabel(iris.feature_names[0])
plt.ylabel(iris.feature_names[1])
plt.title('Original Data (first 2 features)')
plt.colorbar(label='Class')
plt.subplot(1, 2, 2)
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=iris.target, cmap='viridis', alpha=0.7)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
plt.title('PCA Projection')
plt.colorbar(label='Class')
plt.tight_layout()
plt.show()๐ฏ Principal Components
๐ฏ What Are They?
Principal components are the new axes of the transformed space. Each PC is a linear combination of original features, and they are ordered by the amount of variance they explain. They are orthogonal (uncorrelated) to each other.
Properties of Principal Components
Ordered: PC1 has max variance, PC2 next, etc. Orthogonal: Perpendicular to each other. Uncorrelated: Zero correlation between PCs. Linear: Combinations of original features. Loadings: Coefficients of original features. Interpretation: Can analyze what each PC means.
๐ Ordered๐ Orthogonal๐ฏ Uncorrelated๐ Loadings
# Understanding principal components
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
# Load iris
iris = load_iris()
X = iris.data
y = iris.target
feature_names = iris.feature_names
# Apply PCA
pca = PCA()
X_pca = pca.fit_transform(X)
print("=== Explained Variance ===")
for i, (var, ratio) in enumerate(zip(pca.explained_variance_,
pca.explained_variance_ratio_)):
print(f"PC{i+1}: variance={var:.4f}, ratio={ratio:.4f} ({ratio*100:.2f}%)")
print(f"\nCumulative variance:")
cumsum = np.cumsum(pca.explained_variance_ratio_)
for i, cs in enumerate(cumsum):
print(f" PC1-{i+1}: {cs:.4f} ({cs*100:.2f}%)")
# 1. Components (loadings)
print("\n=== Components (Loadings) ===")
loadings_df = pd.DataFrame(
pca.components_.T,
columns=[f'PC{i+1}' for i in range(len(pca.components_))],
index=feature_names
)
print(loadings_df.round(4))
# 2. Interpretation
print("\n=== Interpretation ===")
for i in range(min(4, len(pca.components_))):
component = pca.components_[i]
# Find most important features
importance = np.abs(component)
top_features = np.argsort(importance)[::-1][:3]
print(f"\nPC{i+1}:")
for feat_idx in top_features:
direction = "+" if component[feat_idx] > 0 else "-"
print(f" {direction} {feature_names[feat_idx]}: {component[feat_idx]:+.4f}")
# 3. Loadings visualization
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Scree plot
axes[0].bar(range(1, len(pca.explained_variance_) + 1),
pca.explained_variance_, alpha=0.7, color='steelblue')
axes[0].plot(range(1, len(pca.explained_variance_) + 1),
pca.explained_variance_, 'ro-')
axes[0].set_xlabel('Principal Component')
axes[0].set_ylabel('Explained Variance')
axes[0].set_title('Scree Plot')
axes[0].grid(alpha=0.3)
# Cumulative variance
axes[1].plot(range(1, len(cumsum) + 1), cumsum, 'bo-')
axes[1].axhline(y=0.95, color='r', linestyle='--', label='95% threshold')
axes[1].axhline(y=0.99, color='g', linestyle='--', label='99% threshold')
axes[1].set_xlabel('Number of Components')
axes[1].set_ylabel('Cumulative Explained Variance')
axes[1].set_title('Cumulative Variance')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
# 4. Biplot โ how features relate to PCs
fig, ax = plt.subplots(figsize=(10, 8))
# Project data
X_2d = X_pca[:, :2]
# Plot data points
scatter = ax.scatter(X_2d[:, 0], X_2d[:, 1], c=y, cmap='viridis', alpha=0.5, s=50)
# Plot loading vectors
for i, feature in enumerate(feature_names):
ax.arrow(0, 0,
pca.components_[0, i] * 3,
pca.components_[1, i] * 3,
color='red', alpha=0.7, head_width=0.15)
ax.text(pca.components_[0, i] * 3.3,
pca.components_[1, i] * 3.3,
feature, color='red', fontsize=10, fontweight='bold')
ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
ax.set_title('PCA Biplot')
ax.axhline(0, color='gray', linestyle='--', alpha=0.5)
ax.axvline(0, color='gray', linestyle='--', alpha=0.5)
ax.grid(alpha=0.3)
plt.colorbar(scatter, label='Class')
plt.tight_layout()
plt.show()
# 5. Feature contributions to PC1
print("\n=== Feature Contributions to PC1 ===")
contributions = pd.DataFrame({
'feature': feature_names,
'loading': pca.components_[0],
'abs_loading': np.abs(pca.components_[0])
}).sort_values('abs_loading', ascending=False)
print(contributions.to_string(index=False))
# 6. Reconstructed data with different numbers of PCs
print("\n=== Reconstruction Quality ===")
for k in [1, 2, 3, 4]:
pca_k = PCA(n_components=k)
X_reduced = pca_k.fit_transform(X)
X_reconstructed = pca_k.inverse_transform(X_reduced)
error = np.mean((X - X_reconstructed) ** 2)
variance_kept = pca_k.explained_variance_ratio_.sum()
print(f"k={k}: MSE={error:.6f}, Variance kept={variance_kept:.4f} ({variance_kept*100:.1f}%)")โ๏ธ Feature Scaling
๐ฏ Why Standardization Matters
PCA is sensitive to feature scales. If one feature ranges 0-1 and another 0-1000, PCA will be dominated by the larger-scale feature. Standardization (z-score normalization) is usually essential before PCA.
When to Scale
Always scale: When features have different units. Don't scale: When features are already comparable. StandardScaler: Mean=0, Std=1. MinMaxScaler: Range [0, 1]. RobustScaler: Uses median and IQR. Impact: Changes which directions are "principal".
๐ Different Units๐ฏ Z-Scoreโ๏ธ Balance๐ Impact
# Feature scaling before PCA
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler
# Create data with different scales
np.random.seed(42)
X = np.column_stack([
np.random.randn(100) * 0.1, # Small scale
np.random.randn(100) * 100, # Large scale
np.random.randn(100) * 10 + 50, # Different mean and scale
np.random.randn(100) * 1 # Normal scale
])
print("=== Original Data Statistics ===")
print(f"Mean: {X.mean(axis=0)}")
print(f"Std: {X.std(axis=0)}")
print(f"Range: {X.max(axis=0) - X.min(axis=0)}")
# 1. PCA without scaling
print("\n=== PCA WITHOUT Scaling ===")
pca_no_scale = PCA()
X_pca_no_scale = pca_no_scale.fit_transform(X)
print(f"Explained variance ratios: {pca_no_scale.explained_variance_ratio_}")
print(f"PC1 variance: {pca_no_scale.explained_variance_[0]:.2f}")
print(f"Feature with largest variance dominates PC1!")
# Find which feature dominates
dominant_feature = np.argmax(np.abs(pca_no_scale.components_[0]))
print(f"Dominant feature in PC1: Feature {dominant_feature} (scale: {X[:, dominant_feature].std():.2f})")
# 2. PCA with StandardScaler
print("\n=== PCA WITH StandardScaler ===")
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
print(f"Scaled mean: {X_scaled.mean(axis=0).round(10)}")
print(f"Scaled std: {X_scaled.std(axis=0)}")
pca_scaled = PCA()
X_pca_scaled = pca_scaled.fit_transform(X_scaled)
print(f"\nExplained variance ratios: {pca_scaled.explained_variance_ratio_}")
print(f"PC1 variance: {pca_scaled.explained_variance_[0]:.2f}")
print(f"Now variance is more balanced!")
# 3. Compare components
print("\n=== Component Comparison ===")
print("Without scaling - PC1 loadings:")
print(pca_no_scale.components_[0].round(4))
print("\nWith scaling - PC1 loadings:")
print(pca_scaled.components_[0].round(4))
# 4. Different scalers
scalers = {
'StandardScaler': StandardScaler(),
'MinMaxScaler': MinMaxScaler(),
'RobustScaler': RobustScaler(),
}
print("\n=== Different Scalers ===")
for name, scaler in scalers.items():
X_scaled = scaler.fit_transform(X)
pca = PCA(n_components=2)
pca.fit(X_scaled)
print(f"\n{name}:")
print(f" Explained variance: {pca.explained_variance_ratio_}")
print(f" Total: {pca.explained_variance_ratio_.sum():.4f}")
# 5. When NOT to scale
print("\n=== When NOT to Scale ===")
print("""
1. All features already have same units/scale
Example: pixel intensities (0-255), word counts
2. You want to preserve relative importance
Example: financial data where larger values = more important
3. Data is already standardized
Example: z-scored data from previous step
4. PCA is used for compression, not interpretation
Example: image compression where all pixels are comparable
""")
# 6. Real-world example
print("\n=== Real-World Example ===")
from sklearn.datasets import load_wine
wine = load_wine()
X_wine = wine.data
feature_names = wine.feature_names
print(f"Feature scales in wine dataset:")
for i, name in enumerate(feature_names):
print(f" {name}: mean={X_wine[:, i].mean():.2f}, std={X_wine[:, i].std():.2f}")
# Without scaling
pca_no = PCA(n_components=2)
X_pca_no = pca_no.fit_transform(X_wine)
print(f"\nWithout scaling: PC1 explains {pca_no.explained_variance_ratio_[0]*100:.1f}%")
# With scaling
scaler = StandardScaler()
X_wine_scaled = scaler.fit_transform(X_wine)
pca_yes = PCA(n_components=2)
X_pca_yes = pca_yes.fit_transform(X_wine_scaled)
print(f"With scaling: PC1 explains {pca_yes.explained_variance_ratio_[0]*100:.1f}%")
# 7. Best practices
print("\n=== Best Practices ===")
print("""
1. Always visualize feature distributions first
2. Scale if features have different units or ranges
3. StandardScaler is the most common choice
4. Use RobustScaler for data with outliers
5. Fit scaler on training data only (not test data)
6. Store scaler parameters for production
7. Consider if scaling is meaningful for your problem
""")
# 8. Pipeline
print("\n=== Recommended Pipeline ===")
from sklearn.pipeline import Pipeline
pipeline = Pipeline([
('scaler', StandardScaler()),
('pca', PCA(n_components=0.95)) # Keep 95% variance
])
X_pipeline = pipeline.fit_transform(X_wine)
print(f"Pipeline result shape: {X_pipeline.shape}")
print(f"Components kept: {pipeline.named_steps['pca'].n_components_}")
print(f"Variance explained: {pipeline.named_steps['pca'].explained_variance_ratio_.sum():.4f}")๐ Choosing K Components
๐ฏ How Many Components to Keep?
Choosing the right number of components (K) is critical. Too few and you lose information; too many and you don't reduce dimensionality. Several methods help: variance threshold, scree plot, Kaiser rule, and cross-validation.
Methods for Choosing K
Variance threshold: Keep 95% or 99% of variance. Scree plot: Look for the "elbow". Kaiser rule: Keep eigenvalues > 1. Cross-validation: Optimize downstream task. Domain knowledge: Interpretability concerns. Reconstruction error: Minimize error.
๐ Variance๐ Scree๐ข Kaiser๐ฏ CV
# Choosing K components
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_digits, load_iris
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
# Load a dataset with more features
digits = load_digits()
X = digits.data # 64 features
y = digits.target
print(f"Original features: {X.shape[1]}")
# Fit PCA with all components
pca = PCA().fit(X)
explained_variance = pca.explained_variance_
explained_variance_ratio = pca.explained_variance_ratio_
cumulative_variance = np.cumsum(explained_variance_ratio)
# 1. Variance threshold
print("\n=== Variance Threshold ===")
for threshold in [0.90, 0.95, 0.99, 0.999]:
n_components = np.argmax(cumulative_variance >= threshold) + 1
print(f" {threshold*100:.1f}% variance โ {n_components} components")
# 2. Scree plot (elbow method)
print("\n=== Scree Plot ===")
# Find elbow using second derivative
diffs = np.diff(explained_variance)
diffs2 = np.diff(diffs)
elbow = np.argmax(diffs2) + 2
print(f"Elbow at component: {elbow}")
print(f"Variance at elbow: {cumulative_variance[elbow-1]:.4f}")
# 3. Kaiser rule
print("\n=== Kaiser Rule ===")
kaiser_components = np.sum(explained_variance > 1)
print(f"Components with eigenvalue > 1: {kaiser_components}")
print(f"Variance explained: {cumulative_variance[kaiser_components-1]:.4f}")
# 4. Reconstruction error
print("\n=== Reconstruction Error ===")
errors = []
for k in range(1, 65):
pca_k = PCA(n_components=k)
X_reduced = pca_k.fit_transform(X)
X_reconstructed = pca_k.inverse_transform(X_reduced)
mse = np.mean((X - X_reconstructed) ** 2)
errors.append(mse)
# Find k where error stops decreasing rapidly
errors = np.array(errors)
error_diffs = np.diff(errors)
elbow_error = np.argmax(np.abs(np.diff(error_diffs))) + 2
print(f"Elbow in reconstruction error: {elbow_error}")
# 5. Cross-validation
print("\n=== Cross-Validation ===")
cv_scores = []
k_range = range(1, 65, 5)
for k in k_range:
pipeline = Pipeline([
('scaler', StandardScaler()),
('pca', PCA(n_components=k)),
('clf', LogisticRegression(max_iter=1000))
])
scores = cross_val_score(pipeline, X, y, cv=5, scoring='accuracy')
cv_scores.append(scores.mean())
print(f" k={k}: accuracy={scores.mean():.4f} (ยฑ{scores.std():.4f})")
best_k_cv = list(k_range)[np.argmax(cv_scores)]
print(f"\nBest k by CV: {best_k_cv}")
# 6. Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Scree plot
axes[0, 0].bar(range(1, len(explained_variance) + 1),
explained_variance, alpha=0.7, color='steelblue')
axes[0, 0].plot(range(1, len(explained_variance) + 1),
explained_variance, 'ro-', markersize=4)
axes[0, 0].axhline(y=1, color='red', linestyle='--', label='Kaiser (ฮป=1)')
axes[0, 0].axvline(x=elbow, color='green', linestyle='--', label=f'Elbow ({elbow})')
axes[0, 0].set_xlabel('Principal Component')
axes[0, 0].set_ylabel('Eigenvalue')
axes[0, 0].set_title('Scree Plot')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)
# Cumulative variance
axes[0, 1].plot(range(1, len(cumulative_variance) + 1),
cumulative_variance, 'bo-', markersize=4)
axes[0, 1].axhline(y=0.95, color='r', linestyle='--', label='95%')
axes[0, 1].axhline(y=0.99, color='g', linestyle='--', label='99%')
for thresh, color in [(0.95, 'r'), (0.99, 'g')]:
k = np.argmax(cumulative_variance >= thresh) + 1
axes[0, 1].axvline(x=k, color=color, linestyle=':', alpha=0.5)
axes[0, 1].text(k+1, thresh-0.05, f'k={k}', color=color)
axes[0, 1].set_xlabel('Number of Components')
axes[0, 1].set_ylabel('Cumulative Variance')
axes[0, 1].set_title('Cumulative Explained Variance')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)
# Reconstruction error
axes[1, 0].plot(range(1, 65), errors, 'go-', markersize=4)
axes[1, 0].axvline(x=elbow_error, color='red', linestyle='--',
label=f'Elbow ({elbow_error})')
axes[1, 0].set_xlabel('Number of Components')
axes[1, 0].set_ylabel('Reconstruction MSE')
axes[1, 0].set_title('Reconstruction Error')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)
# Cross-validation
axes[1, 1].plot(k_range, cv_scores, 'mo-', markersize=6)
axes[1, 1].axvline(x=best_k_cv, color='red', linestyle='--',
label=f'Best k={best_k_cv}')
axes[1, 1].set_xlabel('Number of Components')
axes[1, 1].set_ylabel('CV Accuracy')
axes[1, 1].set_title('Cross-Validation Score')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
# 7. Recommendations
print("\n=== Recommendations ===")
print("""
Rule of thumb:
- 95% variance: Most common choice
- 99% variance: When accuracy is critical
- Kaiser rule: Good default for exploration
- Elbow method: Visual, subjective
- CV: Best for supervised learning
For different use cases:
- Visualization: 2-3 components
- Noise reduction: 90-95% variance
- Feature engineering: 95-99% variance
- Compression: As few as possible
- Downstream ML: Optimize via CV
""")
# 8. Practical function
def choose_k(X, method='variance', threshold=0.95):
"""Choose K components based on method."""
pca = PCA().fit(X)
if method == 'variance':
cumsum = np.cumsum(pca.explained_variance_ratio_)
return np.argmax(cumsum >= threshold) + 1
elif method == 'kaiser':
return np.sum(pca.explained_variance_ > 1)
elif method == 'elbow':
diffs2 = np.diff(np.diff(pca.explained_variance_))
return np.argmax(diffs2) + 2
else:
raise ValueError(f"Unknown method: {method}")
print("\n=== Practical Usage ===")
for method in ['variance', 'kaiser', 'elbow']:
k = choose_k(X, method=method)
print(f"{method}: k={k}")๐ Kernel PCA
๐ฏ Non-Linear Dimensionality Reduction
Kernel PCA extends PCA to non-linear data by using the kernel trick. It implicitly maps data to a higher-dimensional space where linear PCA can separate non-linear structures.
Kernel PCA Features
Kernel trick: Compute dot products without explicit mapping. RBF kernel: Gaussian โ most common. Polynomial: For polynomial relationships. Sigmoid: Neural network-like. Use cases: Non-linear manifolds, clustering. Cost: Slower than linear PCA.
๐ Non-linear๐ฏ Kernel Trick๐ RBF๐ฌ Manifolds
# Kernel PCA for non-linear dimensionality reduction
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA, KernelPCA
from sklearn.datasets import make_circles, make_moons, make_swiss_roll
# 1. Concentric circles (non-linear)
X_circles, y_circles = make_circles(n_samples=400, factor=0.3, noise=0.05, random_state=42)
# Linear PCA fails
pca_linear = PCA(n_components=2)
X_linear = pca_linear.fit_transform(X_circles)
# Kernel PCA works
kpca_rbf = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_kpca_rbf = kpca_rbf.fit_transform(X_circles)
# 2. Moons (non-linear)
X_moons, y_moons = make_moons(n_samples=400, noise=0.05, random_state=42)
pca_moons = PCA(n_components=2)
X_pca_moons = pca_moons.fit_transform(X_moons)
kpca_moons = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_kpca_moons = kpca_moons.fit_transform(X_moons)
# 3. Swiss roll (3D โ 2D)
X_swiss, color_swiss = make_swiss_roll(n_samples=1000, noise=0.05, random_state=42)
pca_swiss = PCA(n_components=2)
X_pca_swiss = pca_swiss.fit_transform(X_swiss)
kpca_swiss = KernelPCA(n_components=2, kernel='rbf', gamma=0.1)
X_kpca_swiss = kpca_swiss.fit_transform(X_swiss)
# 4. Comparison plot
fig, axes = plt.subplots(3, 4, figsize=(16, 12))
# Circles
axes[0, 0].scatter(X_circles[:, 0], X_circles[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[0, 0].set_title('Circles - Original')
axes[0, 1].scatter(X_linear[:, 0], X_linear[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[0, 1].set_title('Circles - Linear PCA (fails)')
axes[0, 2].scatter(X_kpca_rbf[:, 0], X_kpca_rbf[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[0, 2].set_title('Circles - Kernel PCA RBF')
axes[0, 3].axis('off')
# Moons
axes[1, 0].scatter(X_moons[:, 0], X_moons[:, 1], c=y_moons, cmap='viridis', alpha=0.6)
axes[1, 0].set_title('Moons - Original')
axes[1, 1].scatter(X_pca_moons[:, 0], X_pca_moons[:, 1], c=y_moons, cmap='viridis', alpha=0.6)
axes[1, 1].set_title('Moons - Linear PCA')
axes[1, 2].scatter(X_kpca_moons[:, 0], X_kpca_moons[:, 1], c=y_moons, cmap='viridis', alpha=0.6)
axes[1, 2].set_title('Moons - Kernel PCA RBF')
axes[1, 3].axis('off')
# Swiss roll
axes[2, 0].scatter(X_swiss[:, 0], X_swiss[:, 1], c=color_swiss, cmap='viridis', alpha=0.6)
axes[2, 0].set_title('Swiss Roll - Original (3D)')
axes[2, 1].scatter(X_pca_swiss[:, 0], X_pca_swiss[:, 1], c=color_swiss, cmap='viridis', alpha=0.6)
axes[2, 1].set_title('Swiss Roll - Linear PCA')
axes[2, 2].scatter(X_kpca_swiss[:, 0], X_kpca_swiss[:, 1], c=color_swiss, cmap='viridis', alpha=0.6)
axes[2, 2].set_title('Swiss Roll - Kernel PCA RBF')
axes[2, 3].axis('off')
plt.tight_layout()
plt.show()
# 5. Different kernels
print("=== Different Kernels ===")
kernels = ['linear', 'poly', 'rbf', 'sigmoid', 'cosine']
for kernel in kernels:
try:
kpca = KernelPCA(n_components=2, kernel=kernel)
X_kpca = kpca.fit_transform(X_circles)
print(f"Kernel: {kernel}")
print(f" Explained variance: {kpca.lambdas_[:3]}")
except Exception as e:
print(f"Kernel: {kernel} โ {e}")
# 6. Gamma parameter (for RBF)
print("\n=== RBF Gamma Parameter ===")
gammas = [0.1, 1, 5, 15, 50, 100]
fig, axes = plt.subplots(1, len(gammas), figsize=(20, 4))
for i, gamma in enumerate(gammas):
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=gamma)
X_kpca = kpca.fit_transform(X_circles)
axes[i].scatter(X_kpca[:, 0], X_kpca[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[i].set_title(f'gamma={gamma}')
axes[i].axis('off')
plt.tight_layout()
plt.show()
print("""
Gamma interpretation:
- Small gamma: Smooth, underfitting
- Medium gamma: Good balance
- Large gamma: Complex, overfitting
- Rule of thumb: gamma โ 1 / (n_features * X.var())
""")
# 7. Kernel PCA vs other non-linear methods
print("\n=== Kernel PCA vs Other Methods ===")
from sklearn.manifold import TSNE, Isomap, LocallyLinearEmbedding
# Apply different methods
tsne = TSNE(n_components=2, random_state=42, perplexity=30)
X_tsne = tsne.fit_transform(X_circles)
isomap = Isomap(n_components=2, n_neighbors=10)
X_isomap = isomap.fit_transform(X_circles)
lle = LocallyLinearEmbedding(n_components=2, n_neighbors=10)
X_lle = lle.fit_transform(X_circles)
fig, axes = plt.subplots(1, 4, figsize=(20, 5))
axes[0].scatter(X_kpca_rbf[:, 0], X_kpca_rbf[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[0].set_title('Kernel PCA (RBF)')
axes[1].scatter(X_tsne[:, 0], X_tsne[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[1].set_title('t-SNE')
axes[2].scatter(X_isomap[:, 0], X_isomap[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[2].set_title('Isomap')
axes[3].scatter(X_lle[:, 0], X_lle[:, 1], c=y_circles, cmap='viridis', alpha=0.6)
axes[3].set_title('LLE')
plt.tight_layout()
plt.show()
# 8. Practical example
print("\n=== Practical Usage ===")
from sklearn.datasets import load_digits
from sklearn.model_selection import cross_val_score
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
digits = load_digits()
X_digits = digits.data
y_digits = digits.target
# Compare linear vs kernel PCA
for name, reducer in [
('Linear PCA', PCA(n_components=30)),
('Kernel PCA (RBF)', KernelPCA(n_components=30, kernel='rbf', gamma=0.001)),
('Kernel PCA (poly)', KernelPCA(n_components=30, kernel='poly', degree=3)),
]:
pipeline = Pipeline([
('reduce', reducer),
('clf', SVC(kernel='linear'))
])
scores = cross_val_score(pipeline, X_digits, y_digits, cv=3, scoring='accuracy')
print(f"{name}: {scores.mean():.4f} (ยฑ{scores.std():.4f})")๐ Incremental PCA
๐ฏ PCA for Large Datasets
Incremental PCA processes data in batches, allowing PCA on datasets that don't fit in memory. It's essential for big data applications.
When to Use Incremental PCA
Large datasets: Data doesn't fit in RAM. Streaming data: Continuous incoming data. Online learning: Model updates incrementally. Memory constraints: Limited RAM. Trade-off: Slightly less accurate than full PCA. Batch size: Balance memory and accuracy.
๐ Batches๐พ Memory๐ Streamingโก Online
# Incremental PCA for large datasets
import numpy as np
import time
from sklearn.decomposition import PCA, IncrementalPCA
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
# Generate a large dataset
n_samples = 100000
n_features = 100
X, y = make_classification(
n_samples=n_samples, n_features=n_features,
n_informative=50, random_state=42
)
print(f"Dataset: {X.shape}, size: {X.nbytes / 1e6:.1f} MB")
# 1. Regular PCA (needs all data in memory)
print("\n=== Regular PCA ===")
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
t0 = time.time()
pca_full = PCA(n_components=50)
X_pca_full = pca_full.fit_transform(X_scaled)
full_time = time.time() - t0
print(f"Time: {full_time:.2f}s")
print(f"Explained variance: {pca_full.explained_variance_ratio_.sum():.4f}")
# 2. Incremental PCA (processes in batches)
print("\n=== Incremental PCA ===")
batch_size = 1000
n_batches = n_samples // batch_size
t0 = time.time()
ipca = IncrementalPCA(n_components=50, batch_size=batch_size)
for i in range(n_batches):
start = i * batch_size
end = start + batch_size
batch = X_scaled[start:end]
ipca.partial_fit(batch)
X_pca_inc = ipca.transform(X_scaled)
inc_time = time.time() - t0
print(f"Time: {inc_time:.2f}s")
print(f"Explained variance: {ipca.explained_variance_ratio_.sum():.4f}")
# 3. Compare results
print("\n=== Comparison ===")
print(f"Full PCA variance: {pca_full.explained_variance_ratio_.sum():.4f}")
print(f"Incremental PCA variance: {ipca.explained_variance_ratio_.sum():.4f}")
print(f"Difference: {abs(pca_full.explained_variance_ratio_.sum() - ipca.explained_variance_ratio_.sum()):.6f}")
# Compare components
component_diff = np.abs(pca_full.components_ - ipca.components_).mean()
print(f"Component difference: {component_diff:.6f}")
# 4. Different batch sizes
print("\n=== Effect of Batch Size ===")
batch_sizes = [100, 500, 1000, 5000, 10000]
for bs in batch_sizes:
ipca = IncrementalPCA(n_components=50, batch_size=bs)
n_batches = max(1, n_samples // bs)
for i in range(n_batches):
batch = X_scaled[i*bs:(i+1)*bs]
ipca.partial_fit(batch)
variance = ipca.explained_variance_ratio_.sum()
print(f"Batch size {bs}: variance={variance:.4f}")
# 5. Streaming data example
print("\n=== Streaming Data ===")
ipca_stream = IncrementalPCA(n_components=50, batch_size=500)
# Simulate streaming data
all_variances = []
for i in range(0, n_samples, 1000):
batch = X_scaled[i:i+1000]
ipca_stream.partial_fit(batch)
# Track cumulative variance
if i % 10000 == 0:
total_var = ipca_stream.explained_variance_ratio_.sum()
all_variances.append(total_var)
print(f"After {i+1000} samples: variance={total_var:.4f}")
# 6. Memory comparison
print("\n=== Memory Usage ===")
print(f"Full dataset: {X.nbytes / 1e6:.1f} MB")
print(f"Per batch ({batch_size} samples): {X[:batch_size].nbytes / 1e6:.2f} MB")
print(f"Memory savings: {X.nbytes / X[:batch_size].nbytes:.0f}x")
# 7. Real-world example: Image data
print("\n=== Real-World: Image Data ===")
from sklearn.datasets import load_digits
digits = load_digits()
X_img = digits.data
print(f"Image data: {X_img.shape}")
# Incremental PCA on images
ipca_img = IncrementalPCA(n_components=32, batch_size=100)
for i in range(0, len(X_img), 100):
ipca_img.partial_fit(X_img[i:i+100])
X_img_pca = ipca_img.transform(X_img)
print(f"Reduced shape: {X_img_pca.shape}")
print(f"Compression ratio: {X_img.shape[1] / X_img_pca.shape[1]:.1f}x")
# Reconstruct
X_img_reconstructed = ipca_img.inverse_transform(X_img_pca)
mse = np.mean((X_img - X_img_reconstructed) ** 2)
print(f"Reconstruction MSE: {mse:.4f}")
# 8. Best practices
print("\n=== Best Practices ===")
print("""
1. Choose batch size carefully:
- Too small: slow convergence
- Too large: memory issues
- Rule: 1000-10000 samples per batch
2. Always fit scaler on full data first
- Or use partial_fit on scaler too
3. Watch for numerical stability:
- Use float64 for accumulation
- Normalize data
4. Validate with holdout set
- Compare to full PCA when possible
5. Use for:
- Datasets > RAM
- Streaming data
- Online learning
- Memory-constrained environments
6. Consider alternatives:
- Randomized PCA (faster for many features)
- Truncated SVD (for sparse data)
""")
# 9. Complete pipeline example
print("\n=== Complete Pipeline ===")
from sklearn.pipeline import Pipeline
# Simulate loading data in chunks
chunk_sizes = [1000] * 10
pipeline_ipca = IncrementalPCA(n_components=20)
print("Processing chunks:")
for i, chunk_size in enumerate(chunk_sizes):
chunk = X_scaled[i*chunk_size:(i+1)*chunk_size]
pipeline_ipca.partial_fit(chunk)
print(f" Chunk {i+1}: {chunk_size} samples, "
f"variance={pipeline_ipca.explained_variance_ratio_.sum():.4f}")
print(f"\nFinal variance: {pipeline_ipca.explained_variance_ratio_.sum():.4f}")
print(f"Components: {pipeline_ipca.n_components_}")๐ค PCA Whitening
๐ฏ Decorrelation + Normalization
Whitening transforms data so that it has identity covariance matrix โ features are uncorrelated and have unit variance. It's used in image processing and as preprocessing for ML algorithms.
Two Types of Whitening
PCA whitening: PCA then divide by eigenvalues. ZCA whitening: PCA whitening then rotate back to original space. Benefit: Removes redundancy. Use case: Image preprocessing, ICA, deep learning. Effect: Identity covariance matrix.
๐ค Whitening๐ ZCA๐ Unit Variance๐ผ๏ธ Images
# PCA whitening and ZCA whitening
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
# Generate correlated data
np.random.seed(42)
n = 500
X = np.random.randn(n, 2) @ np.array([[3, 2], [1, 1]])
# 1. Standard PCA (no whitening)
pca = PCA(n_components=2, whiten=False)
X_pca = pca.fit_transform(X)
# 2. PCA whitening
pca_white = PCA(n_components=2, whiten=True)
X_pca_white = pca_white.fit_transform(X)
# 3. Manual whitening
# Step 1: Center
X_centered = X - X.mean(axis=0)
# Step 2: Compute covariance
cov = np.cov(X_centered.T)
# Step 3: Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# Step 4: PCA whitening (rotate + scale)
X_pca_manual = X_centered @ eigenvectors
X_pca_manual_white = X_pca_manual / np.sqrt(eigenvalues)
# Step 5: ZCA whitening (scale + rotate back)
X_zca_white = X_pca_manual_white @ eigenvectors.T
# Verify
print("=== Covariance Matrices ===")
print(f"Original covariance:\n{np.cov(X.T)}")
print(f"\nAfter PCA (no whitening):\n{np.cov(X_pca.T)}")
print(f"\nAfter PCA whitening:\n{np.cov(X_pca_white.T)}")
print(f"\nAfter ZCA whitening:\n{np.cov(X_zca_white.T)}")
# 4. Compare transformations
fig, axes = plt.subplots(1, 4, figsize=(20, 5))
datasets = [
(X, 'Original Data'),
(X_pca, 'PCA (no whitening)'),
(X_pca_white, 'PCA Whitened'),
(X_zca_white, 'ZCA Whitened'),
]
for ax, (data, title) in zip(axes, datasets):
ax.scatter(data[:, 0], data[:, 1], alpha=0.5, s=20)
ax.set_title(title)
ax.axis('equal')
ax.grid(alpha=0.3)
# Show covariance
cov = np.cov(data.T)
text = f"cov:\n[{cov[0,0]:.2f}, {cov[0,1]:.2f}]\n[{cov[1,0]:.2f}, {cov[1,1]:.2f}]"
ax.text(0.05, 0.95, text, transform=ax.transAxes,
verticalalignment='top', fontsize=8,
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
plt.tight_layout()
plt.show()
# 5. Eigenvalue analysis
print("\n=== Eigenvalue Analysis ===")
print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")
# PCA whitening: divide by sqrt(eigenvalue)
print(f"\nPCA whitened variance: {X_pca_manual_white.var(axis=0)}")
print("Expected: [1, 1]")
# ZCA whitening
print(f"ZCA whitened variance: {X_zca_white.var(axis=0)}")
# 6. Application: Image whitening
print("\n=== Image Whitening ===")
from sklearn.datasets import load_digits
digits = load_digits()
X_img = digits.data
# Normalize to [0, 1]
X_img_norm = X_img / 16.0
# PCA whitening
pca_img = PCA(n_components=64, whiten=True)
X_img_white = pca_img.fit_transform(X_img_norm)
# ZCA whitening (manual)
X_img_centered = X_img_norm - X_img_norm.mean(axis=0)
cov_img = np.cov(X_img_centered.T)
eigenvalues_img, eigenvectors_img = np.linalg.eigh(cov_img)
idx = np.argsort(eigenvalues_img)[::-1]
eigenvalues_img = eigenvalues_img[idx]
eigenvectors_img = eigenvectors_img[:, idx]
# PCA whitening
X_img_pca_white = X_img_centered @ eigenvectors_img
X_img_pca_white = X_img_pca_white / np.sqrt(eigenvalues_img + 1e-5)
# ZCA whitening
X_img_zca_white = X_img_pca_white @ eigenvectors_img.T
# Visualize a sample
sample_idx = 0
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
axes[0].imshow(X_img_norm[sample_idx].reshape(8, 8), cmap='gray')
axes[0].set_title('Original')
axes[1].imshow(X_img_centered[sample_idx].reshape(8, 8), cmap='gray')
axes[1].set_title('Centered')
axes[2].imshow(X_img_pca_white[sample_idx][:64].reshape(8, 8), cmap='gray')
axes[2].set_title('PCA Whitened')
axes[3].imshow(X_img_zca_white[sample_idx].reshape(8, 8), cmap='gray')
axes[3].set_title('ZCA Whitened')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
# 7. When to use which
print("\n=== When to Use Which ===")
print("""
PCA Whitening:
+ Removes correlations
+ Unit variance per component
+ Good for ICA, clustering
- Doesn't preserve original space
ZCA Whitening:
+ Removes correlations
+ Unit variance
+ Preserves original orientation
+ Better for CNNs, image processing
- More computation
Use PCA whitening when:
- You want maximally decorrelated features
- Downstream algorithm assumes independent features
- Doing ICA or independent component analysis
Use ZCA whitening when:
- Preserving visual structure matters
- Preprocessing images for CNNs
- You want decorrelated but recognizable features
- Working with spatial data
Don't use whitening when:
- Feature scale is meaningful
- Sparse data (whitening destroys sparsity)
- Interpretability is important
""")
# 8. Whitening in deep learning
print("\n=== Whitening in Deep Learning ===")
print("""
Used in:
- CNN input preprocessing (ZCA)
- Batch normalization (a form of whitening)
- Style transfer (whitening and coloring)
- Domain adaptation
- Self-supervised learning
Modern alternatives:
- Batch Normalization (per-layer)
- Layer Normalization
- Instance Normalization
- Group Normalization
Batch Norm is essentially:
1. Compute mean and variance per batch
2. Normalize
3. Learn scale and shift parameters
""")
# 9. Complete whitening function
def whiten(X, method='zca', epsilon=1e-5):
"""Whiten data using PCA or ZCA."""
# Center
X_centered = X - X.mean(axis=0)
# Covariance
cov = np.cov(X_centered.T)
# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# Whitening matrix
D = np.diag(1.0 / np.sqrt(eigenvalues + epsilon))
W_pca = D @ eigenvectors.T
W_zca = eigenvectors @ D @ eigenvectors.T
if method == 'pca':
return X_centered @ W_pca.T
elif method == 'zca':
return X_centered @ W_zca.T
else:
raise ValueError(f"Unknown method: {method}")
# Test
X_test = np.random.randn(100, 3) @ np.array([[2, 1, 0], [0, 2, 1], [0, 0, 1]])
X_pca_w = whiten(X_test, method='pca')
X_zca_w = whiten(X_test, method='zca')
print("\n=== Whitening Function Test ===")
print(f"Original covariance:\n{np.cov(X_test.T).round(3)}")
print(f"\nPCA whitened covariance:\n{np.cov(X_pca_w.T).round(3)}")
print(f"\nZCA whitened covariance:\n{np.cov(X_zca_w.T).round(3)}")๐ผ Applications (Vision, NLP, Finance)
๐ฏ Where PCA Shines
PCA is used across many domains: computer vision (Eigenfaces, image compression), NLP (word embeddings, topic modeling), finance (risk factors, portfolio optimization), and more.
Domain Applications
Vision: Eigenfaces, image compression, denoising. NLP: Latent Semantic Analysis, word embeddings. Finance: Risk factors, portfolio optimization. Biology: Gene expression, population genetics. Robotics: Sensor fusion, state estimation. Time series: Trend extraction, anomaly detection.
๐๏ธ Vision๐ NLP๐ฐ Finance๐งฌ Biology
# Real-world PCA applications
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import fetch_olivetti_faces, load_digits
# 1. EIGENFACES (Face Recognition)
print("=== Eigenfaces ===")
faces = fetch_olivetti_faces(shuffle=True, random_state=42)
X_faces = faces.data
y_faces = faces.target
print(f"Face data: {X_faces.shape}")
# Apply PCA
n_components = 50
pca_faces = PCA(n_components=n_components, whiten=True)
X_faces_pca = pca_faces.fit_transform(X_faces)
print(f"Reduced shape: {X_faces_pca.shape}")
print(f"Variance explained: {pca_faces.explained_variance_ratio_.sum():.4f}")
# Visualize eigenfaces
fig, axes = plt.subplots(4, 5, figsize=(12, 10))
for i, ax in enumerate(axes.flat):
eigenface = pca_faces.components_[i].reshape(64, 64)
ax.imshow(eigenface, cmap='gray')
ax.set_title(f'Eigenface {i+1}')
ax.axis('off')
plt.suptitle('Top 20 Eigenfaces', fontsize=16)
plt.tight_layout()
plt.show()
# Face recognition with PCA
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score
X_train, X_test, y_train, y_test = train_test_split(
X_faces_pca, y_faces, test_size=0.3, random_state=42
)
clf = SVC(kernel='rbf', C=10, gamma=0.001)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
print(f"Face recognition accuracy: {accuracy_score(y_test, y_pred):.4f}")
# 2. IMAGE COMPRESSION
print("\n=== Image Compression ===")
digits = load_digits()
X_digits = digits.data / 16.0 # Normalize
# Try different compression levels
k_values = [5, 10, 20, 30, 50, 64]
fig, axes = plt.subplots(2, len(k_values), figsize=(20, 6))
sample_idx = 0
original = X_digits[sample_idx].reshape(8, 8)
# Original
axes[0, 0].imshow(original, cmap='gray')
axes[0, 0].set_title('Original')
axes[0, 0].axis('off')
for i, k in enumerate(k_values):
if k == 5:
continue
pca_k = PCA(n_components=k)
X_compressed = pca_k.fit_transform(X_digits)
X_reconstructed = pca_k.inverse_transform(X_compressed)
reconstructed = X_reconstructed[sample_idx].reshape(8, 8)
ax = axes[0, i]
ax.imshow(reconstructed, cmap='gray')
ax.set_title(f'k={k}')
ax.axis('off')
# Compression ratio
original_size = 64
compressed_size = k * 2 + k # components + mean
ratio = original_size / compressed_size
ax = axes[1, i]
ax.bar(['Original', 'Compressed'], [original_size, compressed_size])
ax.set_title(f'Ratio: {ratio:.1f}x')
plt.tight_layout()
plt.show()
# 3. LATENT SEMANTIC ANALYSIS (NLP)
print("\n=== Latent Semantic Analysis ===")
from sklearn.feature_extraction.text import TfidfVectorizer
documents = [
"The cat sat on the mat",
"The dog played in the park",
"Machine learning is fascinating",
"Deep learning uses neural networks",
"Cats and dogs are pets",
"AI and ML are related fields",
"The cat chased the mouse",
"Neural networks learn patterns",
"Natural language processing is hard",
"Pets need love and care"
]
vectorizer = TfidfVectorizer(stop_words='english')
X_tfidf = vectorizer.fit_transform(documents).toarray()
print(f"TF-IDF matrix: {X_tfidf.shape}")
# Apply LSA (PCA on TF-IDF)
pca_lsa = PCA(n_components=3)
X_lsa = pca_lsa.fit_transform(X_tfidf)
print(f"LSA shape: {X_lsa.shape}")
print(f"Variance explained: {pca_lsa.explained_variance_ratio_.sum():.4f}")
# Visualize documents in 2D
fig, ax = plt.subplots(figsize=(10, 8))
ax.scatter(X_lsa[:, 0], X_lsa[:, 1], s=100, alpha=0.6)
for i, doc in enumerate(documents):
ax.annotate(f"Doc {i+1}", (X_lsa[i, 0], X_lsa[i, 1]), fontsize=9)
ax.set_xlabel(f'PC1 ({pca_lsa.explained_variance_ratio_[0]*100:.1f}%)')
ax.set_ylabel(f'PC2 ({pca_lsa.explained_variance_ratio_[1]*100:.1f}%)')
ax.set_title('Document Clustering with LSA')
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# 4. FINANCE: PORTFOLIO RISK FACTORS
print("\n=== Finance: Portfolio Risk ===")
np.random.seed(42)
n_days = 1000
n_stocks = 50
# Simulate stock returns
# Market factor + sector factors + idiosyncratic
market_factor = np.random.randn(n_days) * 0.01
sector_factors = np.random.randn(n_days, 5) * 0.005
idiosyncratic = np.random.randn(n_days, n_stocks) * 0.02
# Stock returns = exposure to factors + noise
X_stocks = (market_factor[:, None] * np.random.randn(n_stocks) +
sector_factors @ np.random.randn(5, n_stocks) * 0.5 +
idiosyncratic)
# Apply PCA
pca_finance = PCA()
X_finance = pca_finance.fit_transform(X_stocks)
print(f"Stock returns: {X_stocks.shape}")
print(f"Top 5 PC variances: {pca_finance.explained_variance_[:5]}")
# Variance explained by first few factors
cumulative = np.cumsum(pca_finance.explained_variance_ratio_)
print(f"\nVariance explained by top factors:")
for i in [1, 3, 5, 10, 20]:
print(f" Top {i:2d} factors: {cumulative[i-1]*100:.2f}%")
# 5. TIME SERIES: TREND EXTRACTION
print("\n=== Time Series: Trend Extraction ===")
np.random.seed(42)
n_points = 500
# Generate multiple time series with common trend
common_trend = np.sin(np.linspace(0, 4*np.pi, n_points)) + np.linspace(0, 1, n_points)
noise = np.random.randn(n_points, 10) * 0.3
X_timeseries = noise + common_trend[:, None]
# Apply PCA
pca_ts = PCA(n_components=3)
X_ts_pca = pca_ts.fit_transform(X_timeseries)
print(f"Time series: {X_timeseries.shape}")
print(f"Variance explained: {pca_ts.explained_variance_ratio_}")
# Plot
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
axes[0].plot(X_timeseries[:, :5], alpha=0.5)
axes[0].set_title('Original Time Series (5 shown)')
axes[0].set_xlabel('Time')
axes[0].set_ylabel('Value')
axes[0].grid(alpha=0.3)
axes[1].plot(X_ts_pca[:, :3], linewidth=2)
axes[1].set_title('First 3 Principal Components')
axes[1].set_xlabel('Time')
axes[1].set_ylabel('Component Value')
axes[1].legend([f'PC{i+1}' for i in range(3)])
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
# 6. BIOLOGY: GENE EXPRESSION
print("\n=== Biology: Gene Expression ===")
# Simulate gene expression data
n_samples = 200
n_genes = 1000
# Simulate: few biological processes + noise
n_processes = 5
process_activity = np.random.randn(n_samples, n_processes)
gene_loadings = np.random.randn(n_processes, n_genes)
noise = np.random.randn(n_samples, n_genes) * 0.5
X_genes = process_activity @ gene_loadings + noise
# PCA
pca_genes = PCA(n_components=20)
X_genes_pca = pca_genes.fit_transform(X_genes)
print(f"Gene expression: {X_genes.shape}")
print(f"Top PCs explain:")
cumsum_genes = np.cumsum(pca_genes.explained_variance_ratio_)
for k in [1, 3, 5, 10]:
print(f" PC1-{k}: {cumsum_genes[k-1]*100:.2f}%")
# 7. ROBOTICS: SENSOR FUSION
print("\n=== Robotics: Sensor Fusion ===")
# Simulate redundant sensor readings
n_timesteps = 1000
n_sensors = 20
# True state + sensor noise
true_state = np.cumsum(np.random.randn(n_timesteps) * 0.1)
sensor_readings = true_state[:, None] + np.random.randn(n_timesteps, n_sensors) * 0.5
# PCA to extract common signal
pca_sensors = PCA(n_components=1)
X_sensors_pca = pca_sensors.fit_transform(sensor_readings)
# The first PC should recover the true state
correlation = np.corrcoef(X_sensors_pca[:, 0], true_state)[0, 1]
print(f"Sensor data: {sensor_readings.shape}")
print(f"PC1 correlation with true state: {correlation:.4f}")
print(f"Variance explained: {pca_sensors.explained_variance_ratio_[0]*100:.2f}%")
# 8. Applications summary
print("\n\n=== PCA Applications Summary ===")
applications = {
'Computer Vision': [
'Eigenfaces for face recognition',
'Image compression',
'Image denoising',
'Feature extraction'
],
'NLP': [
'Latent Semantic Analysis',
'Word embeddings',
'Topic modeling',
'Document similarity'
],
'Finance': [
'Risk factor modeling',
'Portfolio optimization',
'Interest rate curves',
'Credit scoring'
],
'Biology': [
'Gene expression analysis',
'Population genetics',
'Protein structure',
'Single-cell RNA-seq'
],
'Robotics': [
'Sensor fusion',
'State estimation',
'SLAM',
'Motion planning'
],
'Time Series': [
'Trend extraction',
'Seasonality decomposition',
'Anomaly detection',
'Factor models'
]
}
for domain, uses in applications.items():
print(f"\n{domain}:")
for use in uses:
print(f" โข {use}")โ๏ธ Advantages & Limitations
๐ฏ When to Use PCA
PCA is a powerful, interpretable technique โ but it's not always the right choice. Understanding its strengths and weaknesses helps you decide when to use it and when to reach for alternatives.
Trade-offs
Pros: Fast, interpretable, removes correlation, reduces noise. Cons: Linear only, assumes Gaussian, may lose information, hard to interpret components. Alternatives: t-SNE, UMAP, autoencoders, kernel PCA, ICA. Best for: Linear data, preprocessing, compression, exploration.
โก Fast๐ Interpretable๐ Linear๐ Alternatives
# PCA vs alternatives โ comparison
import numpy as np
import time
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA, KernelPCA, FastICA, NMF
from sklearn.manifold import TSNE, Isomap, LocallyLinearEmbedding
from sklearn.datasets import load_digits, make_swiss_roll
from sklearn.preprocessing import StandardScaler
# Load data
digits = load_digits()
X = digits.data
y = digits.target
# Standardize
X_scaled = StandardScaler().fit_transform(X)
print(f"Data: {X.shape}")
# Compare methods
methods = {
'PCA': PCA(n_components=2),
'Kernel PCA (RBF)': KernelPCA(n_components=2, kernel='rbf', gamma=0.001),
'FastICA': FastICA(n_components=2, random_state=42),
'NMF': NMF(n_components=2, random_state=42, init='nndsvd', max_iter=500),
't-SNE': TSNE(n_components=2, random_state=42, perplexity=30),
'Isomap': Isomap(n_components=2, n_neighbors=10),
'LLE': LocallyLinearEmbedding(n_components=2, n_neighbors=10),
}
results = {}
timings = {}
for name, method in methods.items():
try:
t0 = time.time()
X_transformed = method.fit_transform(X_scaled)
elapsed = time.time() - t0
results[name] = X_transformed
timings[name] = elapsed
print(f"{name}: {elapsed:.2f}s")
except Exception as e:
print(f"{name}: failed โ {e}")
# Visualize
n_methods = len(results)
fig, axes = plt.subplots(2, 4, figsize=(20, 10))
axes = axes.flatten()
for i, (name, X_transformed) in enumerate(results.items()):
if i >= len(axes):
break
axes[i].scatter(X_transformed[:, 0], X_transformed[:, 1],
c=y, cmap='tab10', alpha=0.6, s=10)
axes[i].set_title(f'{name}\n({timings[name]:.2f}s)')
axes[i].axis('off')
# Hide extra axes
for j in range(len(results), len(axes)):
axes[j].axis('off')
plt.tight_layout()
plt.show()
# Comparison table
print("\n=== Comparison Table ===")
print(f"{'Method':<20} {'Time':<10} {'Linear':<10} {'Preserves':<15}")
print("-" * 60)
method_props = {
'PCA': ('Yes', 'Variance'),
'Kernel PCA (RBF)': ('No', 'Non-linear structure'),
'FastICA': ('No', 'Independence'),
'NMF': ('No', 'Parts (non-negative)'),
't-SNE': ('No', 'Local structure'),
'Isomap': ('No', 'Geodesic distances'),
'LLE': ('No', 'Local linear structure'),
}
for name, (linear, preserves) in method_props.items():
if name in timings:
print(f"{name:<20} {timings[name]:<10.2f} {linear:<10} {preserves:<15}")
# 1. PCA ADVANTAGES
print("\n=== PCA Advantages ===")
print("""
1. FAST
- O(min(n,p)ยฒ ร max(n,p))
- Faster than t-SNE, UMAP
- Scales to large datasets
2. INTERPRETABLE
- Components are linear combinations
- Loadings show feature contributions
- Explained variance is meaningful
3. DETERMINISTIC
- Same input โ same output
- No random initialization
- Reproducible
4. DECORRELATION
- Removes multicollinearity
- Uncorrelated features
- Diagonal covariance
5. NOISE REDUCTION
- Removes directions with low variance
- Effective for noisy data
- Acts as regularizer
6. COMPRESSION
- Reduces storage
- Speeds up downstream ML
- Preserves most information
""")
# 2. PCA LIMITATIONS
print("\n=== PCA Limitations ===")
print("""
1. LINEAR ONLY
- Cannot capture non-linear relationships
- Fails on curved manifolds
- Use Kernel PCA for non-linear
2. ASSUMES GAUSSIAN
- Assumes data from normal distribution
- May not work for skewed data
- Consider transformations first
3. SENSITIVE TO SCALING
- Must standardize features
- Results depend on scaling
- Different scales โ different PCs
4. LOSES INTERPRETABILITY
- PCs don't have physical meaning
- Complex feature combinations
- Hard to explain to stakeholders
5. LOSES INFORMATION
- Reducing dimensions loses info
- Choose K carefully
- Reconstruction error grows
6. OUTLIER SENSITIVE
- Outliers affect covariance
- Use RobustScaler
- Consider robust PCA
7. ASSUMES LINEAR CORRELATION
- Only captures linear relationships
- Misses complex dependencies
- Consider other methods
""")
# 3. WHEN TO USE PCA
print("\n=== When to Use PCA ===")
print("""
USE PCA WHEN:
โ Data is linearly correlated
โ You need fast dimensionality reduction
โ You want to visualize high-dim data
โ You need to remove multicollinearity
โ You want to reduce noise
โ You have large sparse data
โ Interpretability is not critical
โ You need deterministic results
โ Preprocessing for supervised learning
โ Feature extraction for clustering
DON'T USE PCA WHEN:
โ Data has non-linear structure (use t-SNE, UMAP)
โ You need to preserve local structure (use LLE, Isomap)
โ Features have different meaning (you lose interpretability)
โ You need to preserve the original features
โ Data is categorical (use MCA)
โ Outliers are present (use robust PCA)
โ You need probabilistic interpretation (use PPCA)
โ Scale matters (don't standardize)
ALTERNATIVES TO CONSIDER:
- t-SNE: Visualization, local structure
- UMAP: Visualization + preservation
- Autoencoders: Non-linear, deep learning
- ICA: Independent components
- NMF: Non-negative, parts-based
- LDA: Supervised dimensionality reduction
""")
# 4. Practical decision tree
print("\n=== Decision Tree ===")
print("""
Data type:
โโโ Linear + Gaussian โ Use PCA
โโโ Non-linear โ Consider Kernel PCA, t-SNE, or UMAP
โโโ Sparse โ Consider TruncatedSVD
โโโ Non-negative โ Consider NMF
โโโ Time series โ Consider Dynamic PCA
โโโ Supervised โ Consider LDA or PLS
Goal:
โโโ Visualization โ t-SNE, UMAP, PCA
โโโ Preprocessing โ PCA (fast)
โโโ Compression โ PCA, autoencoders
โโโ Feature extraction โ PCA, autoencoders
โโโ Interpretability โ PCA (with limitations)
โโโ Non-linear embedding โ Kernel PCA, autoencoders
""")
# 5. Hybrid approaches
print("\n=== Hybrid Approaches ===")
print("""
Often, the best solution is a combination:
1. PCA โ t-SNE
- PCA first to reduce dimensions (e.g., 50)
- t-SNE on the reduced data
- Faster than t-SNE on original
2. PCA โ Clustering
- Reduce dimensions with PCA
- Cluster in PC space
- Faster and often better
3. PCA โ Neural Network
- Preprocess with PCA
- Feed to neural network
- Faster training
4. PCA โ UMAP
- Similar to PCA + t-SNE
- Preserves more global structure
5. Randomized PCA for large data
- Approximate PCA
- Much faster
- Minimal accuracy loss
""")
# 6. Summary recommendation
print("\n\n=== RECOMMENDATION ===")
print("""
For most cases:
1. START with PCA
- Fast, interpretable, well-understood
- Good baseline for dimensionality reduction
2. Standardize your features
- Unless units are identical
3. Choose K components carefully
- 95% variance is a good default
- Use cross-validation when possible
4. If PCA doesn't work well:
- Check for non-linearity
- Try Kernel PCA
- Consider t-SNE for visualization
- Consider autoencoders for deep learning
5. Combine methods
- PCA + t-SNE for visualization
- PCA + clustering for analysis
- PCA + ML for predictive modeling
Remember:
- PCA is a tool, not a solution
- Understand your data first
- Try multiple methods
- Validate with downstream tasks
""")๐ Future: Deep PCA & Autoencoders
๐ฏ The Next Generation
PCA is evolving: deep autoencoders, variational autoencoders, self-supervised learning, and neural dimensionality reduction are pushing the boundaries of what's possible.
Modern Alternatives
Autoencoders: Neural networks that learn compressed representations. VAEs: Probabilistic autoencoders with smooth latent space. Diffusion models: For generation and dimensionality reduction. Contrastive learning: SimCLR, MoCo for representation learning. Non-linear PCA: Kernel methods, deep learning.
๐ง Autoencoders๐ฒ VAE๐ Deep๐ค Neural
# Deep dimensionality reduction โ autoencoders vs PCA
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.decomposition import PCA
from sklearn.datasets import load_digits
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
# Load data
digits = load_digits()
X = digits.data
y = digits.target
# Standardize
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# Convert to tensor
X_tensor = torch.FloatTensor(X_scaled)
print(f"Data: {X.shape}")
# 1. LINEAR AUTOENCODER (equivalent to PCA)
class LinearAutoencoder(nn.Module):
def __init__(self, input_dim, latent_dim):
super().__init__()
self.encoder = nn.Linear(input_dim, latent_dim)
self.decoder = nn.Linear(latent_dim, input_dim)
def forward(self, x):
z = self.encoder(x)
return self.decoder(z), z
# Train linear autoencoder
linear_ae = LinearAutoencoder(X.shape[1], 2)
optimizer = optim.Adam(linear_ae.parameters(), lr=0.001)
criterion = nn.MSELoss()
dataset = TensorDataset(X_tensor)
loader = DataLoader(dataset, batch_size=64, shuffle=True)
losses = []
for epoch in range(100):
epoch_loss = 0
for batch in loader:
x = batch[0]
optimizer.zero_grad()
x_reconstructed, _ = linear_ae(x)
loss = criterion(x_reconstructed, x)
loss.backward()
optimizer.step()
epoch_loss += loss.item()
losses.append(epoch_loss / len(loader))
# Get latent representation
with torch.no_grad():
_, X_linear_ae = linear_ae(X_tensor)
X_linear_ae = X_linear_ae.numpy()
# 2. DEEP AUTOENCODER (non-linear)
class DeepAutoencoder(nn.Module):
def __init__(self, input_dim, latent_dim=2):
super().__init__()
self.encoder = nn.Sequential(
nn.Linear(input_dim, 128),
nn.ReLU(),
nn.Linear(128, 64),
nn.ReLU(),
nn.Linear(64, 32),
nn.ReLU(),
nn.Linear(32, latent_dim)
)
self.decoder = nn.Sequential(
nn.Linear(latent_dim, 32),
nn.ReLU(),
nn.Linear(32, 64),
nn.ReLU(),
nn.Linear(64, 128),
nn.ReLU(),
nn.Linear(128, input_dim)
)
def forward(self, x):
z = self.encoder(x)
return self.decoder(z), z
# Train deep autoencoder
deep_ae = DeepAutoencoder(X.shape[1], 2)
optimizer = optim.Adam(deep_ae.parameters(), lr=0.001)
for epoch in range(200):
epoch_loss = 0
for batch in loader:
x = batch[0]
optimizer.zero_grad()
x_reconstructed, _ = deep_ae(x)
loss = criterion(x_reconstructed, x)
loss.backward()
optimizer.step()
epoch_loss += loss.item()
with torch.no_grad():
_, X_deep_ae = deep_ae(X_tensor)
X_deep_ae = X_deep_ae.numpy()
# 3. VARIATIONAL AUTOENCODER (VAE)
class VAE(nn.Module):
def __init__(self, input_dim, latent_dim=2):
super().__init__()
self.encoder = nn.Sequential(
nn.Linear(input_dim, 128),
nn.ReLU(),
nn.Linear(128, 64),
nn.ReLU()
)
self.fc_mu = nn.Linear(64, latent_dim)
self.fc_logvar = nn.Linear(64, latent_dim)
self.decoder = nn.Sequential(
nn.Linear(latent_dim, 64),
nn.ReLU(),
nn.Linear(64, 128),
nn.ReLU(),
nn.Linear(128, input_dim)
)
def encode(self, x):
h = self.encoder(x)
return self.fc_mu(h), self.fc_logvar(h)
def reparameterize(self, mu, logvar):
std = torch.exp(0.5 * logvar)
eps = torch.randn_like(std)
return mu + eps * std
def forward(self, x):
mu, logvar = self.encode(x)
z = self.reparameterize(mu, logvar)
return self.decoder(z), z, mu, logvar
# Train VAE
vae = VAE(X.shape[1], 2)
optimizer = optim.Adam(vae.parameters(), lr=0.001)
for epoch in range(200):
for batch in loader:
x = batch[0]
optimizer.zero_grad()
x_reconstructed, z, mu, logvar = vae(x)
# VAE loss: reconstruction + KL divergence
recon_loss = criterion(x_reconstructed, x)
kl_loss = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
kl_loss = kl_loss / x.size(0)
loss = recon_loss + 0.001 * kl_loss
loss.backward()
optimizer.step()
with torch.no_grad():
_, X_vae, _, _ = vae(X_tensor)
X_vae = X_vae.numpy()
# 4. Compare all methods
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)
fig, axes = plt.subplots(2, 4, figsize=(20, 10))
methods = {
'PCA (Linear)': X_pca,
'Linear AE': X_linear_ae,
'Deep AE': X_deep_ae,
'VAE': X_vae,
}
for ax, (name, X_transformed) in zip(axes.flatten(), methods.items()):
scatter = ax.scatter(X_transformed[:, 0], X_transformed[:, 1],
c=y, cmap='tab10', alpha=0.6, s=10)
ax.set_title(name)
ax.axis('off')
plt.tight_layout()
plt.show()
# 5. Reconstruction error comparison
def reconstruction_error(original, reconstructed):
return np.mean((original - reconstructed) ** 2)
# For PCA
X_pca_reconstructed = pca.inverse_transform(X_pca)
pca_error = reconstruction_error(X_scaled, X_pca_reconstructed)
# For autoencoders
with torch.no_grad():
X_linear_reconstructed, _ = linear_ae(X_tensor)
X_deep_reconstructed, _ = deep_ae(X_tensor)
X_vae_reconstructed, _, _, _ = vae(X_tensor)
linear_ae_error = reconstruction_error(X_scaled, X_linear_reconstructed.numpy())
deep_ae_error = reconstruction_error(X_scaled, X_deep_reconstructed.numpy())
vae_error = reconstruction_error(X_scaled, X_vae_reconstructed.numpy())
print("\n=== Reconstruction Error (MSE) ===")
print(f"PCA: {pca_error:.4f}")
print(f"Linear AE: {linear_ae_error:.4f}")
print(f"Deep AE: {deep_ae_error:.4f}")
print(f"VAE: {vae_error:.4f}")
# 6. Future of dimensionality reduction
print("\n=== Future of Dimensionality Reduction ===")
print("""
1. DEEP LEARNING APPROACHES
- Deep autoencoders
- Variational autoencoders (VAE)
- Generative models (GANs, diffusion)
- Self-supervised learning
2. NEURAL DIMENSIONALITY REDUCTION
- UMAP with neural networks
- Neural PCA
- Deep kernel methods
- Contrastive learning
3. HYBRID METHODS
- PCA + Autoencoders
- Linear + Non-linear
- Probabilistic methods
4. FOUNDATION MODELS
- Pre-trained representations
- CLIP, DINOv2 embeddings
- Transfer learning
- Zero-shot dimensionality reduction
5. ADVANCED METHODS
- Topological data analysis
- Manifold learning
- Sparse coding
- Dictionary learning
6. PRACTICAL TRENDS
- On-device dimensionality reduction
- Federated learning with PCA
- Differential privacy
- Streaming algorithms
""")
# 7. Comparison: PCA vs Deep Learning
print("\n=== PCA vs Deep Learning ===")
comparison = {
'Speed': ('Very fast', 'Slower (needs GPU)'),
'Interpretability': ('High', 'Low'),
'Non-linearity': ('No', 'Yes'),
'Data requirements': ('Any size', 'Large data preferred'),
'Training': ('No training', 'Requires training'),
'Deterministic': ('Yes', 'No (random init)'),
'Hyperparameters': ('Few', 'Many'),
'Deployment': ('Simple', 'Complex'),
'Robustness': ('Good', 'Variable'),
'Best for': ('Linear, fast, interpretable', 'Complex, non-linear, large'),
}
print(f"{'Aspect':<25} {'PCA':<25} {'Deep Learning':<25}")
print("-" * 75)
for aspect, (pca_val, dl_val) in comparison.items():
print(f"{aspect:<25} {pca_val:<25} {dl_val:<25}")
# 8. Recommendations
print("\n=== Recommendations ===")
print("""
Use PCA when:
โ Data is linearly correlated
โ You need fast, interpretable results
โ Dataset is small or medium
โ You want deterministic output
โ You can't afford GPU training
โ Interpretability matters
โ You need a baseline
Use Autoencoders when:
โ Data has non-linear structure
โ You have large datasets
โ You need the best reconstruction
โ You can afford training time
โ Interpretability is less important
โ You have GPU resources
Use VAE when:
โ You need a probabilistic model
โ You want to generate new samples
โ You need a smooth latent space
โ You're doing representation learning
โ You want to interpolate between samples
โ You need uncertainty estimates
Hybrid approach (recommended):
1. Start with PCA as baseline
2. If results are insufficient, try autoencoders
3. Compare reconstruction error
4. Consider task-specific metrics
5. Combine PCA + deep learning when possible
""")