PCA — Dimensionality Reduction

You are photographing a 3D sculpture to put on a website. You can only take one photo. From most angles you capture the full shape — height, width, and some sense of depth. From a few bad angles the sculpture looks like a flat line. The best angle is the one that shows the most variation — where the sculpture looks most different from one end to the other.

PCA finds the "best angle" to project high-dimensional data onto a lower-dimensional space — the projection that preserves the most variation. The first principal component is the direction of maximum variance. The second is the direction of maximum remaining variance perpendicular to the first. And so on.

import numpy as np
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
n = 1000

# Amazon customer features — correlated by design
order_freq   = np.random.normal(8, 3, n)
avg_spend    = 200 + 80 * order_freq + np.random.normal(0, 150, n)   # correlated
n_categories = 1 + 0.3 * order_freq + np.random.normal(0, 1, n)
sessions     = 5 + 1.2 * order_freq + np.random.normal(0, 3, n)
return_rate  = np.random.uniform(0, 0.4, n)   # independent
coupon_rate  = np.random.uniform(0, 1, n)      # independent

X = np.column_stack([order_freq, avg_spend, n_categories,
                     sessions, return_rate, coupon_rate])
feat_names = ['order_freq', 'avg_spend', 'n_categories',
              'sessions', 'return_rate', 'coupon_rate']

# ── PCA from scratch — every step visible ─────────────────────────────

# Step 1: Standardise
sc    = StandardScaler()
X_std = sc.fit_transform(X)
print(f"Step 1 — After standardisation:")
print(f"  Mean (should be ~0): {X_std.mean(axis=0).round(3)}")
print(f"  Std  (should be ~1): {X_std.std(axis=0).round(3)}")

# Step 2: Covariance matrix
C = np.cov(X_std.T)   # (n_features, n_features)
print(f"
Step 2 — Covariance matrix shape: {C.shape}")
print("  High covariance between correlated features:")
for i, fi in enumerate(feat_names):
    for j, fj in enumerate(feat_names):
        if i < j and abs(C[i, j]) > 0.3:
            print(f"    {fi} ↔ {fj}: {C[i,j]:.3f}")

# Step 3: Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(C)   # eigh for symmetric matrices
# Sort descending
idx          = np.argsort(eigenvalues)[::-1]
eigenvalues  = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

print(f"
Step 3 — Eigenvalues (variance per component):")
total_var = eigenvalues.sum()
for i, (ev, vec) in enumerate(zip(eigenvalues, eigenvectors.T)):
    pct = ev / total_var * 100
    bar = '█' * int(pct / 3)
    print(f"  PC{i+1}: {bar} {pct:.1f}%  (λ={ev:.4f})")

# Step 4: Project onto top 2 components
k = 2
V         = eigenvectors[:, :k]        # (n_features, k)
X_reduced = X_std @ V                  # (n_samples, k)
print(f"
Step 4 — Projection:")
print(f"  Original shape: {X.shape}   →   Reduced shape: {X_reduced.shape}")
print(f"  Variance retained: {eigenvalues[:k].sum()/total_var*100:.1f}%")

# ── Verify against sklearn PCA ─────────────────────────────────────────
from sklearn.decomposition import PCA
pca_sk = PCA(n_components=2)
X_sk   = pca_sk.fit_transform(X_std)
print(f"
sklearn PCA variance retained: {pca_sk.explained_variance_ratio_.sum()*100:.1f}%")
print(f"(Matches scratch implementation — same result)")

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
n = 2000

# Amazon customer: 20 features, many correlated
base1 = np.random.normal(0, 1, n)   # purchase behaviour axis
base2 = np.random.normal(0, 1, n)   # browsing behaviour axis
base3 = np.random.normal(0, 1, n)   # price sensitivity axis

X = np.column_stack([
    base1 + np.random.normal(0, 0.3, n),           # order_frequency
    base1 * 80 + 800 + np.random.normal(0, 100, n), # avg_spend
    base1 * 3 + 8 + np.random.normal(0, 1, n),     # orders_last_90d
    base1 * 0.5 + np.random.normal(0, 0.2, n),     # repeat_rate
    base2 + np.random.normal(0, 0.3, n),            # app_sessions
    base2 * 5 + 15 + np.random.normal(0, 3, n),    # pages_per_session
    base2 * 0.3 + np.random.normal(0, 0.2, n),     # search_rate
    base2 * 2 + np.random.normal(0, 1, n),         # wishlist_adds
    base3 + np.random.normal(0, 0.3, n),            # coupon_usage
    base3 * 0.4 + np.random.normal(0, 0.2, n),     # discount_sensitivity
    base3 * (-0.3) + np.random.normal(0, 0.2, n),  # premium_brand_pref
    np.random.normal(0, 1, n),                      # return_rate (independent)
    np.random.normal(0, 1, n),                      # review_rate
    np.random.normal(0, 1, n),                      # rating_avg
    np.random.normal(0, 1, n),                      # device_switches
    np.random.normal(0, 1, n),                      # time_of_day_pref
    np.random.normal(0, 1, n),                      # n_payment_methods
    np.random.normal(0, 1, n),                      # referral_count
    np.random.normal(0, 1, n),                      # support_contacts
    np.random.normal(0, 1, n),                      # account_age_days
])

sc   = StandardScaler()
X_sc = sc.fit_transform(X)

# ── Full PCA — examine all components ─────────────────────────────────
pca_full = PCA()
pca_full.fit(X_sc)

evr  = pca_full.explained_variance_ratio_
cumr = np.cumsum(evr)

print("Explained variance per component:")
print(f"{'PC':<6} {'Variance':>10} {'Cumulative':>12} {'Bar'}")
print("─" * 50)
for i, (var, cum) in enumerate(zip(evr, cumr)):
    bar  = '█' * int(var * 100)
    flag = ' ← 95% threshold' if cum >= 0.95 and (i == 0 or cumr[i-1] < 0.95) else ''
    flag = flag or (' ← 99% threshold' if cum >= 0.99 and (i == 0 or cumr[i-1] < 0.99) else '')
    print(f"  PC{i+1:<3}  {var*100:>8.2f}%  {cum*100:>10.2f}%  {bar}{flag}")

# ── Choose k automatically ────────────────────────────────────────────
k_95  = np.searchsorted(cumr, 0.95) + 1
k_99  = np.searchsorted(cumr, 0.99) + 1
k_90  = np.searchsorted(cumr, 0.90) + 1

print(f"
Components needed to retain:")
print(f"  90% variance: k={k_90}  (from {X.shape[1]} features)")
print(f"  95% variance: k={k_95}  (from {X.shape[1]} features)")
print(f"  99% variance: k={k_99}  (from {X.shape[1]} features)")

# ── Apply PCA with chosen k ────────────────────────────────────────────
pca = PCA(n_components=k_95)
X_reduced = pca.fit_transform(X_sc)

print(f"
Reduced dataset shape: {X.shape} → {X_reduced.shape}")
print(f"Variance retained: {pca.explained_variance_ratio_.sum()*100:.2f}%")

# ── n_components as variance fraction ─────────────────────────────────
# Instead of specifying k, specify minimum variance to retain
pca_auto = PCA(n_components=0.95)   # keep enough components for 95%
X_auto   = pca_auto.fit_transform(X_sc)
print(f"
PCA(n_components=0.95): kept {pca_auto.n_components_} components")

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
n = 2000

feat_names = [
    'order_frequency', 'avg_spend', 'orders_90d', 'repeat_rate',
    'app_sessions', 'pages_per_session', 'search_rate', 'wishlist_adds',
    'coupon_usage', 'discount_sensitivity', 'premium_brand_pref',
    'return_rate', 'review_rate', 'rating_avg', 'device_switches',
    'time_of_day_pref', 'n_payment_methods', 'referral_count',
    'support_contacts', 'account_age_days',
]

base1 = np.random.normal(0, 1, n)
base2 = np.random.normal(0, 1, n)
base3 = np.random.normal(0, 1, n)

X = np.column_stack([
    base1 + np.random.normal(0,0.3,n),
    base1*80+800+np.random.normal(0,100,n),
    base1*3+8+np.random.normal(0,1,n),
    base1*0.5+np.random.normal(0,0.2,n),
    base2+np.random.normal(0,0.3,n),
    base2*5+15+np.random.normal(0,3,n),
    base2*0.3+np.random.normal(0,0.2,n),
    base2*2+np.random.normal(0,1,n),
    base3+np.random.normal(0,0.3,n),
    base3*0.4+np.random.normal(0,0.2,n),
    base3*(-0.3)+np.random.normal(0,0.2,n),
    *[np.random.normal(0,1,n) for _ in range(9)],
])

sc    = StandardScaler()
X_sc  = sc.fit_transform(X)
pca   = PCA(n_components=5)
X_pca = pca.fit_transform(X_sc)

# ── Loadings matrix ────────────────────────────────────────────────────
# pca.components_: shape (n_components, n_features)
# components_[i, j] = loading of feature j on component i
loadings = pd.DataFrame(
    pca.components_.T,     # transpose: (n_features, n_components)
    index   = feat_names,
    columns = [f'PC{i+1}' for i in range(5)],
)

print("Top contributing features per principal component:")
for pc in ['PC1', 'PC2', 'PC3']:
    top3_pos = loadings[pc].nlargest(3)
    top3_neg = loadings[pc].nsmallest(3)
    print(f"
  {pc} ({pca.explained_variance_ratio_[int(pc[2])-1]*100:.1f}% variance):")
    for feat, val in top3_pos.items():
        bar = '▸' * int(abs(val) * 20)
        print(f"    +{bar} {val:+.3f}  {feat}")
    for feat, val in top3_neg.items():
        bar = '◂' * int(abs(val) * 20)
        print(f"    -{bar} {val:+.3f}  {feat}")

# ── Visualise loadings as heatmap values ───────────────────────────────
print("
Loadings heatmap (top 8 features, first 3 PCs):")
top_feats = loadings.abs().max(axis=1).nlargest(8).index
subset    = loadings.loc[top_feats, ['PC1','PC2','PC3']]

print(f"
{'Feature':<25} {'PC1':>8} {'PC2':>8} {'PC3':>8}")
print("─" * 52)
for feat, row in subset.iterrows():
    vals = ""
    for v in row:
        # Colorize output: Green for positive, Red for negative
        color = "32" if v > 0.1 else "31" if v < -0.1 else "0"
        vals += f"  [{color}m{v:+.3f}"
    print(f"  {feat:<23} {vals}")

Fitting PCA on the full dataset before splitting is data leakage — the same as fitting a scaler on all data before splitting. The principal components are computed using test set statistics. Always use Pipeline so PCA is refit on the training fold in each cross-validation split, just like every other transformer.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.metrics import accuracy_score
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 2000

# Simulated high-dimensional Amazon classification:
# Predict segment (0=budget, 1=mid, 2=premium) from 50 features
base1 = np.random.normal(0, 1, n)
base2 = np.random.normal(0, 1, n)
y     = np.where(base1 > 1, 2, np.where(base1 < -1, 0, 1))

X_core = np.column_stack([base1, base2, base1+base2])
X_noise = np.random.normal(0, 1, (n, 47))   # 47 noise dimensions
X = np.hstack([X_core, X_noise])   # 50 total features

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2,
                                            stratify=y, random_state=42)

# ── Compare: no PCA vs PCA before KNN ─────────────────────────────────
print("KNN: effect of PCA preprocessing on 50-dimensional data:")
print(f"{'Pipeline':<35} {'CV Accuracy'}")
print("─" * 48)

for n_comp in [None, 30, 20, 10, 5, 3]:
    if n_comp is None:
        pipe = Pipeline([
            ('scaler', StandardScaler()),
            ('model',  KNeighborsClassifier(n_neighbors=15)),
        ])
        label = 'No PCA (50 features)'
    else:
        pipe = Pipeline([
            ('scaler', StandardScaler()),
            ('pca',    PCA(n_components=n_comp)),
            ('model',  KNeighborsClassifier(n_neighbors=15)),
        ])
        label = f'PCA → {n_comp} components'

    cv = cross_val_score(pipe, X_tr, y_tr, cv=5, scoring='accuracy').mean()
    print(f"  {label:<33}  {cv:.4f}")

# ── LogisticRegression — PCA removes multicollinearity ────────────────
print("
LogisticRegression: effect of PCA on correlated features:")
for n_comp in [None, 10, 5, 3]:
    if n_comp is None:
        pipe = Pipeline([
            ('scaler', StandardScaler()),
            ('model',  LogisticRegression(max_iter=1000, random_state=42)),
        ])
        label = 'No PCA'
    else:
        pipe = Pipeline([
            ('scaler', StandardScaler()),
            ('pca',    PCA(n_components=n_comp)),
            ('model',  LogisticRegression(max_iter=1000, random_state=42)),
        ])
        label = f'PCA({n_comp})'

    cv = cross_val_score(pipe, X_tr, y_tr, cv=5, scoring='accuracy').mean()
    print(f"  {label:<15}: CV accuracy = {cv:.4f}")

# ── PCA for visualisation — project to 2D ─────────────────────────────
sc_viz   = StandardScaler()
X_viz_sc = sc_viz.fit_transform(X_tr)
pca_viz  = PCA(n_components=2)
X_2d     = pca_viz.fit_transform(X_viz_sc)

print(f"
2D visualisation:")
print(f"  Original: {X_tr.shape}  →  2D: {X_2d.shape}")
print(f"  Variance retained: {pca_viz.explained_variance_ratio_.sum()*100:.1f}%")
print("  Can now plot clusters in a scatter plot — 2 axes instead of 50")

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
n = 1000

base1 = np.random.normal(0, 1, n)
base2 = np.random.normal(0, 1, n)
X = np.column_stack([
    base1+np.random.normal(0,0.3,n),
    base1*2+np.random.normal(0,0.5,n),
    base2+np.random.normal(0,0.3,n),
    base2*1.5+np.random.normal(0,0.4,n),
    *[np.random.normal(0,1,n) for _ in range(16)],
])   # 20 features

sc   = StandardScaler()
X_sc = sc.fit_transform(X)

# ── Reconstruct from k components and measure error ───────────────────
print("Reconstruction error vs components kept:")
print(f"{'k components':<16} {'Variance kept':>14} {'MSE (recon)':>13} {'RMSE':>8}")
print("─" * 56)

for k in [1, 2, 3, 5, 8, 10, 15, 20]:
    pca   = PCA(n_components=k)
    X_red = pca.fit_transform(X_sc)         # compress: (n, 20) → (n, k)
    X_rec = pca.inverse_transform(X_red)    # reconstruct: (n, k) → (n, 20)

    mse  = np.mean((X_sc - X_rec) ** 2)
    rmse = np.sqrt(mse)
    var  = pca.explained_variance_ratio_.sum()

    bar  = '█' * int(var * 30)
    print(f"  k={k:<3}  {k:>3}/{X.shape[1]}    {bar:<30} {var*100:>5.1f}%    {mse:>8.4f}   {rmse:>6.4f}")

# ── Anomaly detection using reconstruction error ───────────────────────
# Normal customers: low reconstruction error (fit the PCA model well)
# Anomalous customers: high error (don't fit the patterns PCA found)

X_normal  = X_sc[:950]
X_anomaly = np.random.normal(5, 3, (50, 20))   # outliers — very different pattern
X_test    = np.vstack([X_normal, X_anomaly])
true_labels = np.array([0]*950 + [1]*50)        # 0=normal, 1=anomaly

# Fit PCA on "normal" training data
pca_anom = PCA(n_components=5)
pca_anom.fit(X_sc)   # fit on all — in practice, fit on known normal data only

# Reconstruction error per sample
X_rec_test = pca_anom.inverse_transform(pca_anom.transform(X_test))
rec_errors = np.mean((X_test - X_rec_test) ** 2, axis=1)

# High reconstruction error = anomaly
threshold = np.percentile(rec_errors[:950], 95)   # 95th percentile of normal errors
predicted_anomaly = (rec_errors > threshold).astype(int)

tp = ((predicted_anomaly == 1) & (true_labels == 1)).sum()
fp = ((predicted_anomaly == 1) & (true_labels == 0)).sum()
print(f"
Anomaly detection via reconstruction error:")
print(f"  Threshold (95th pct of normal): {threshold:.4f}")
print(f"  True positives (anomalies caught): {tp}/50")
print(f"  False positives (normal flagged):  {fp}/950")

import numpy as np
from sklearn.decomposition import PCA, KernelPCA, TruncatedSVD
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 2000

base = np.random.normal(0, 1, (n, 3))
y    = (base[:, 0] + base[:, 1] > 0).astype(int)
X_corr = np.hstack([
    base + np.random.normal(0, 0.2, (n, 3)),
    base * 2 + np.random.normal(0, 0.3, (n, 3)),
    base * 0.5 + np.random.normal(0, 0.1, (n, 3)),
    np.random.normal(0, 1, (n, 11)),
])   # 20 features, 9 correlated

# ── PCA helps logistic regression but not Random Forest ───────────────
print("Effect of PCA on different algorithms:")
print(f"{'Algorithm':<35} {'No PCA':>10} {'With PCA':>10} {'Change'}")
print("─" * 60)

sc = StandardScaler()

for name, model in [
    ('LogisticRegression', LogisticRegression(max_iter=1000, random_state=42)),
    ('RandomForest',       RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)),
]:
    # Without PCA
    pipe_no_pca = Pipeline([('sc', StandardScaler()), ('m', model)])
    cv_no  = cross_val_score(pipe_no_pca, X_corr, y, cv=5).mean()

    # With PCA (keep 95% variance)
    pipe_pca = Pipeline([('sc', StandardScaler()), ('pca', PCA(n_components=0.95)), ('m', model)])
    cv_pca = cross_val_score(pipe_pca, X_corr, y, cv=5).mean()

    change = cv_pca - cv_no
    arrow  = '↑' if change > 0.005 else '↓' if change < -0.005 else '→'
    print(f"  {name:<33}  {cv_no:.4f}    {cv_pca:.4f}  {arrow} {abs(change):.4f}")

# ── Sparse data: TruncatedSVD instead of PCA ──────────────────────────
# PCA requires centering which destroys sparsity of sparse matrices
# TruncatedSVD = PCA without centering — works on sparse matrices
from scipy.sparse import random as sparse_random
X_sparse = sparse_random(n, 500, density=0.05, random_state=42)
y_sparse = np.random.randint(0, 2, n)

svd = TruncatedSVD(n_components=20, random_state=42)
X_dense = svd.fit_transform(X_sparse)
print(f"
TruncatedSVD on sparse matrix:")
print(f"  {X_sparse.shape} sparse → {X_dense.shape} dense")
print(f"  Variance retained: {svd.explained_variance_ratio_.sum()*100:.1f}%")

# ── Kernel PCA — non-linear dimensionality reduction ──────────────────
from sklearn.datasets import make_circles
X_rings, y_rings = make_circles(n_samples=500, noise=0.05, factor=0.4, random_state=42)

# Standard PCA cannot separate rings (linear)
pca_linear = PCA(n_components=2)
X_lin      = pca_linear.fit_transform(X_rings)

# KernelPCA with RBF — projects to space where rings are separable
kpca   = KernelPCA(n_components=2, kernel='rbf', gamma=10)
X_kpca = kpca.fit_transform(X_rings)

lr_lin  = LogisticRegression().fit(X_lin, y_rings)
lr_kpca = LogisticRegression().fit(X_kpca, y_rings)
print(f"
Ring-shaped clusters — PCA vs KernelPCA:")
print(f"  Standard PCA accuracy:  {lr_lin.score(X_lin, y_rings):.4f}")
print(f"  KernelPCA accuracy:     {lr_kpca.score(X_kpca, y_rings):.4f}")

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.pipeline import Pipeline
import joblib, warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 10_000

# ── High-dimensional Amazon customer dataset: 50 features ───────────
purchase_base  = np.random.normal(0, 1, n)
browsing_base  = np.random.normal(0, 1, n)
price_base     = np.random.normal(0, 1, n)

feature_matrix = np.column_stack([
    # Purchase behaviour (15 features — all correlated with purchase_base)
    *[purchase_base + np.random.normal(0, 0.3+i*0.05, n) for i in range(15)],
    # Browsing behaviour (15 features — all correlated with browsing_base)
    *[browsing_base + np.random.normal(0, 0.3+i*0.05, n) for i in range(15)],
    # Price sensitivity (10 features)
    *[price_base + np.random.normal(0, 0.3+i*0.05, n) for i in range(10)],
    # Truly independent features (10 features)
    *[np.random.normal(0, 1, n) for _ in range(10)],
])

col_names = (
    [f'purchase_{i:02d}' for i in range(15)]
    + [f'browsing_{i:02d}' for i in range(15)]
    + [f'price_sens_{i:02d}' for i in range(10)]
    + [f'independent_{i:02d}' for i in range(10)]
)
df_customers = pd.DataFrame(feature_matrix, columns=col_names)
print(f"Original dataset: {df_customers.shape}")

# ── Step 1: Standardise ────────────────────────────────────────────────
sc   = StandardScaler()
X_sc = sc.fit_transform(df_customers)

# ── Step 2: Find the right number of PCA components ───────────────────
pca_full = PCA().fit(X_sc)
cumvar   = np.cumsum(pca_full.explained_variance_ratio_)
k_95     = np.searchsorted(cumvar, 0.95) + 1
k_90     = np.searchsorted(cumvar, 0.90) + 1

print(f"
Variance analysis:")
print(f"  Original features:  {df_customers.shape[1]}")
print(f"  Components for 90%: {k_90}")
print(f"  Components for 95%: {k_95}")
print(f"  Compression ratio:  {df_customers.shape[1]/k_95:.1f}× reduction")

# ── Step 3: Apply PCA ─────────────────────────────────────────────────
pca      = PCA(n_components=k_95)
X_pca    = pca.fit_transform(X_sc)
print(f"
Reduced dataset: {X_pca.shape}")

# ── Step 4: Cluster in PCA space ──────────────────────────────────────
# Clustering is faster and more meaningful in reduced space
print("
K-Means clustering in PCA space:")
for k in [3, 4, 5, 6]:
    km  = KMeans(n_clusters=k, random_state=42, n_init=10)
    lbl = km.fit_predict(X_pca)
    sil = silhouette_score(X_pca, lbl, sample_size=2000)
    print(f"  k={k}: silhouette={sil:.4f}")

# ── Step 5: Final model ────────────────────────────────────────────────
best_k   = 4
km_final = KMeans(n_clusters=best_k, random_state=42, n_init=20)
df_customers['segment'] = km_final.fit_predict(X_pca)

# Profile segments using ORIGINAL features (more interpretable than PCA features)
profile = df_customers.groupby('segment').agg(
    n_customers  = ('purchase_00', 'count'),
    purchase_avg = ('purchase_00', 'mean'),
    browsing_avg = ('browsing_00', 'mean'),
    price_avg    = ('price_sens_00', 'mean'),
).round(3)
print("
Segment profiles (original feature means):")
print(profile.to_string())

# ── Step 6: Save everything needed for production ─────────────────────
# New customers: standardise → PCA transform → predict segment
bundle = {
    'scaler':    sc,
    'pca':       pca,
    'kmeans':    km_final,
    'n_features_original': df_customers.shape[1] - 1,  # exclude 'segment'
    'n_components_pca':    k_95,
    'version':   'v1.0',
}
joblib.dump(bundle, '/tmp/amazon_pca_segments.pkl')

# Assign new customer to a segment
new_customer = np.random.normal(0, 1, (1, 50))
bundle_loaded = joblib.load('/tmp/amazon_pca_segments.pkl')
X_new_sc  = bundle_loaded['scaler'].transform(new_customer)
X_new_pca = bundle_loaded['pca'].transform(X_new_sc)
segment   = bundle_loaded['kmeans'].predict(X_new_pca)[0]
print(f"
New customer → Segment {segment}")

print("
Production bundle saved: /tmp/amazon_pca_segments.pkl")

A principal component is a weighted combination of every original feature, not a stand-in for any single one. PC1 might be 0.6×avg_spend + 0.4×order_freq − 0.3×return_rate + small contributions from every other column. It is not "avg_spend, renamed." This is exactly why the loadings matrix exists — you have to look at which features load heavily on a component to give it a human label like "purchasing intensity," and even then it is an approximation, not an identity. Feed PCA output into a model and you get strong predictions but zero built-in interpretability at the individual-feature level — that trade is the price of compression, not a bug.

PCA does not detect whether your data has linear or nonlinear structure — it simply assumes linear structure and projects onto the best linear approximation it can find, whether or not that assumption is true. Feed it data shaped like two concentric rings or a Swiss roll and PCA will still confidently return two components — it just will not separate the classes, because no straight-line projection can unroll a curved manifold. That failure mode is precisely why t-SNE, UMAP, and Kernel PCA exist: they model nonlinear neighborhood structure that a linear projection structurally cannot represent. A low reconstruction error or a clean scree plot only ever tells you that *linear* redundancy was compressed well, never that the underlying structure was linear.

They share an algorithm but not a goal. For dimensionality reduction before a model, you choose k to retain a target amount of variance (typically 90–95%) because the objective is preserving predictive signal — k is often 10, 20, or more. For visualisation, you are constrained to k=2 or k=3 regardless of how much variance that retains, because a human can only look at 2 or 3 axes — that projection might capture only 35% of total variance and still be the only option. A 2D PCA scatter plot that shows overlapping clusters does not necessarily mean the classes are inseparable; it may just mean the separating signal lives in components 3 through 10, which the plot discarded. Never conclude "no structure" from a 2D PCA plot alone.

PCA maximises variance, and variance is not scale-invariant — a feature measured in dollars (variance in the millions) will mathematically dominate a feature measured as a proportion (variance under 1) regardless of which one is more informative. Without standardisation, PC1 is not "the direction of maximum meaningful variation" — it is frequently just "the feature with the largest units," which has nothing to do with correlation structure. This is not a minor accuracy loss; it can make the entire eigendecomposition reflect measurement units instead of genuine relationships in the data, silently invalidating everything downstream — the scree plot, the loadings, the reduced dataset.

Explained variance measures how much of the *original data's spread* you kept, not how much *useful signal* you kept. The components with the smallest eigenvalues are often dominated by measurement noise rather than meaningful variation — this is exactly why reconstruction-error anomaly detection works: those low-variance directions are largely noise, so a point that violates them looks anomalous. Chasing 99% instead of 95% explained variance can mean keeping components that are mostly noise, which can add variance to your inputs without adding predictive value, and in some cases hurts downstream model performance by re-introducing the very redundancy and noise PCA was supposed to remove. More retained variance is not synonymous with a better dataset.

"You have many correlated variables and want to summarise them with fewer numbers without throwing away information. PCA finds new axes — linear combinations of your original variables — ordered so the first axis captures as much of the total variance as possible, the second captures the most variance left over after removing the first (and is uncorrelated with it), and so on. Mathematically these axes are the eigenvectors of the covariance matrix, and their eigenvalues tell you how much variance each one captures. You keep the first k axes and discard the rest, trading a small, controlled amount of information loss for a much lower-dimensional dataset."

PCA finds the directions of maximum variance, and variance scales with the square of a feature's units — so a feature measured in dollars (thousands of units of spread) will swamp a feature measured as a 0–1 ratio in the covariance matrix, regardless of which one actually correlates more meaningfully with everything else. Skip standardisation and your first principal component often just tracks whichever raw feature happens to have the largest numeric range, not the direction with the richest correlation structure. StandardScaler before PCA (mean 0, variance 1 per feature) puts every feature on equal footing so the eigendecomposition reflects genuine relationships instead of arbitrary units of measurement.

Start with the cumulative explained-variance curve: fit PCA with no component limit, look at the scree plot, and find where the curve elbows or crosses a threshold like 90–95% — sklearn's PCA(n_components=0.95) does this automatically. But I would not stop there for a supervised task — I'd wrap PCA and the downstream model in a Pipeline and cross-validate the actual metric I care about (AUC, RMSE) across a few candidate values of k, because variance retained and predictive performance retained are correlated but not identical. For a visualisation use case, k is fixed at 2 or 3 regardless of the variance curve, since the deliverable is a plot a human can look at, not a model input.

Three cases: tree-based models (Random Forest, XGBoost) already handle correlated features natively via their split logic and don't benefit from PCA's rotation — it can even hurt them by destroying the original, interpretable split boundaries. Situations requiring feature-level interpretability or regulatory explainability — a compliance team cannot act on "component 3 increased," they need "credit score decreased." And genuinely nonlinear data, where I'd reach for Kernel PCA, t-SNE, or UMAP instead, since standard PCA can only ever find linear projections. I'd also skip it for naturally sparse data like one-hot encodings or text — TruncatedSVD avoids PCA's mean-centering step, which would otherwise destroy the sparsity and blow up memory.

PCA is a linear, deterministic, globally-variance-maximising projection — it preserves large-scale distances but has no way to represent nonlinear structure like clusters that curve around each other. t-SNE and UMAP are nonlinear methods that instead optimise to preserve local neighborhood relationships — points close together in the original space stay close in the projection — which is exactly what makes them good at revealing cluster structure that PCA's straight-line projection would flatten or overlap. The trade-off: t-SNE/UMAP distances between distant clusters are not meaningful and results vary between runs (stochastic optimisation, hyperparameter sensitive), whereas PCA is reproducible and its axes have an actual mathematical meaning (variance explained). In practice I'd use PCA first as a fast, cheap sanity check, and reach for UMAP when PCA's 2D plot shows a blob with no visible separation.