Dimensionality Reduction

import numpy as np

rng = np.random.default_rng(0)
for d in (2, 10, 100, 1000):
    X = rng.standard_normal((2000, d))
    q = rng.standard_normal(d)
    dist = np.linalg.norm(X - q, axis=1)
    contrast = (dist.max() - dist.min()) / dist.min()
    print(f"d={d:>5}  nearest={dist.min():7.2f}  "
          f"farthest={dist.max():7.2f}  contrast={contrast:5.2f}")

# Volume: the fraction of a unit ball within 10% of its surface is 1 - 0.9**d.
[(d, round(1 - 0.9 ** d, 3)) for d in (2, 10, 50, 100)]

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

pca = make_pipeline(StandardScaler(), PCA(n_components=2, random_state=0))
Z = pca.fit_transform(X_train)          # shape (n, 2): the projected coordinates

p = pca[-1]                              # the PCA step inside the pipeline
p.components_                            # shape (2, n_features): the directions themselves
p.explained_variance_                    # the variance along each, in original units
Z.var(axis=0)                            # the same numbers, measured on the output

# Every component is orthogonal to every other, by construction.
np.round(p.components_ @ p.components_.T, 10)   # the 2x2 identity

import numpy as np
from sklearn.decomposition import PCA

full = PCA().fit(X_train_scaled)          # every component, to see the whole curve
evr = full.explained_variance_ratio_
cum = np.cumsum(evr)
k95 = int(np.searchsorted(cum, 0.95) + 1)
print(f"{k95} components reach 95% of the variance ({cum[k95 - 1]:.3f})")

# Or state the budget and let PCA choose k for you:
pca = PCA(n_components=0.95, svd_solver="full", random_state=0).fit(X_train_scaled)
pca.n_components_

import numpy as np
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(0)
wide = rng.normal(0, 2.0, (600, 1))            # huge variance, no signal
y = rng.integers(0, 2, 600)
thin = (y * 1.6 - 0.8 + rng.normal(0, 0.3, 600)).reshape(-1, 1)   # tiny variance, all signal
X = np.hstack([wide, thin])

cross_val_score(LogisticRegression(), X, y, cv=5).mean()           # about 0.99
Z = PCA(n_components=1).fit_transform(X)
cross_val_score(LogisticRegression(), Z, y, cv=5).mean()           # about 0.50

# LDA maximises class separation instead of variance, and keeps the signal.
LinearDiscriminantAnalysis(n_components=1).fit_transform(X, y)

from sklearn.manifold import TSNE
import umap                                # pip install umap-learn

# Reduce with PCA first: it denoises and makes the neighbour search tractable.
from sklearn.decomposition import PCA
X50 = PCA(n_components=50, random_state=0).fit_transform(X_scaled)

Z_tsne = TSNE(n_components=2, perplexity=30, init="pca", random_state=0).fit_transform(X50)
Z_umap = umap.UMAP(n_neighbors=15, min_dist=0.1, random_state=0).fit_transform(X50)

# t-SNE in scikit-learn has no .transform(): there is no way to place a new point.
mapper = umap.UMAP(n_neighbors=15, random_state=0).fit(X50_train)
mapper.transform(X50_test)                 # UMAP can, which is why it is usable in a pipeline

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

# income in currency units, age in years, score in 0-1.
raw = PCA(n_components=2).fit(X_raw)
np.round(raw.components_[0], 4)          # about [1.0, 0.0002, 0.0] — income, alone
raw.explained_variance_ratio_            # about [0.999, 0.001]

scaled = make_pipeline(StandardScaler(), PCA(n_components=2)).fit(X_raw)
np.round(scaled[-1].components_[0], 3)   # a real blend of all three

import numpy as np
from sklearn.decomposition import PCA

pca = PCA(n_components=7, random_state=0).fit(X_train_scaled)
Z = pca.transform(X_test_scaled)
X_hat = pca.inverse_transform(Z)

mse = np.mean((X_test_scaled - X_hat) ** 2)
print(f"reconstruction MSE {mse:.4f}  ~= 1 - {pca.explained_variance_ratio_.sum():.4f}")

# Per-row error: large means 'this row does not live in the training subspace'.
row_error = np.mean((X_test_scaled - X_hat) ** 2, axis=1)
flagged = row_error > np.quantile(row_error, 0.99)

import numpy as np
import pandas as pd

loadings = pd.DataFrame(
    pca.components_.T,
    index=feature_names,
    columns=[f"PC{i + 1}" for i in range(pca.n_components_)],
)
# The features that define each component, largest absolute loading first.
for pc in loadings.columns[:2]:
    top = loadings[pc].abs().sort_values(ascending=False).head(4).index
    print(pc, [(f, round(loadings.loc[f, pc], 2)) for f in top])

# SparsePCA forces most loadings to zero, buying back some readability.
from sklearn.decomposition import SparsePCA
SparsePCA(n_components=3, alpha=1.0, random_state=0).fit(X_train_scaled)

from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

pipe = Pipeline([
    ("scale", StandardScaler()),
    ("pca", PCA(random_state=0)),
    ("clf", LogisticRegression(max_iter=2000)),
])
grid = {
    "pca__n_components": [2, 5, 10, 20, None],   # None = no reduction, the baseline
    "clf__C": [0.1, 1.0, 10.0],
}
search = GridSearchCV(pipe, grid, cv=5, scoring="roc_auc").fit(X_train, y_train)
search.best_params_        # if it picks None, the reduction was not helping